Skip to content

effects package --- add abline to plot

7 messages · Prew, Paul, David Winsemius, John Fox

#
Dear Paul,

I can't think of a way to do what you want using the plot() methods provided
by the package. (You don't say what kind of model you're fitting, but this
is true for all of the plot() methods.) On the other hand, as described in
?effect, the object returned by effect() contains all of the information
needed to construct a custom graph, and you should be able to easily get
what you want.

I hope this helps,
 John

------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox
On
package,
#
On Apr 28, 2009, at 12:00 AM, Prew, Paul wrote:

            
I have concerns that you are considering these probabilities. They are  
not going to be probabilities. They are effects.
I do not know why that is happening and you have not provided a  
minimal executable example. The vectorized use of abline does succeed  
in a simper example:

 > plot(.5,.5)
 > abline(h=c(0.1,0.2,0.3,0.4,0.5,0.6))

  .... so the problem may lie in how the effects package completes its  
plot function for this particular object. You ought to provide at a  
minimum the results of str on that object. Perhaps it executes a  
device call and then turns off the device? However I loaded the  
effects package and ran that abline call after the example:

 > mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion,
+     data=Cowles, family=binomial)
 > eff.cowles <- allEffects(mod.cowles, xlevels=list(neuroticism=0:24,
+     extraversion=seq(0, 24, 6)))
 > eff.cowles


I did not get what I expected, which would have been a single  
horizontal line at 0.4 but rather got four lines roughly at 0.351,  
0.378, 0.408, 0.439. Even then, I would have expected one more line  
before the upper limits of that plot, which makes me think these four  
lines were the results of arguments 0.3 ,0.4, 0.5, 0.6.   Most R  
plotting is done in the coordinate system rather than with absolute  
coordinates, but perhaps the mixture of base graphics with lattice  
graphis is ht eproblem
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
#
Dear David,
The plot() methods in the effects package make lattice graphs which in most instances will have more than one panel. For a binomial GLM, the default is to plot on the scale of the linear predictor (e.g., the logit scale) but to label the response axis on the scale of the response (i.e., the probability scale). To draw a line on the graph, even if you could do it, would require that you translate to the scale of the linear predictor [e.g., for a logit model, log(p/(1 - p))].
"ticks: a two-item list controlling the placement of tick marks on the vertical axis, with elements at and n. If at=NULL (the default), the program attempts to find `nice' locations for the ticks, and the value of n (default, 5) gives the approximate number of tick marks desired; if at is non-NULL, then the value of n is ignored."
You can't just modify a lattice graph on the screen like that. plot() invisibly returns the lattice object; I suppose that you could try to modify that, but I think that my original suggestion -- to make a custom plot from the object returned by effect() -- is likely simpler.

John
#
Dear John and David,  thank you for your help.  I apologize for not defining the analysis as an ordinal regression, or including a structure --- could have taken some of the guesswork out for you.

John --- for the ticks, I would still like to make this work for future analyses, but still not sure what specifically needs changing.  Before initially posting, I did read ?effect, and several other searches around "tick" and "at", but couldn't find a workable description or example for how to use "at".  The one example I did find I thought I had copied pretty closely with my command below.

plot(..., ticks=c(0.1,0.2,0.3,0.4,0.5,0.6))

I tried other combinations that probably look pretty silly...
plot(..., ticks(at=c(0.1,0.2,0.3,0.4,0.5,0.6)))
plot(..., ticks=c(0.1:0.6/0.1))

Just don't know how to properly populate this list ---  
"ticks: a two-item list controlling the placement of tick marks on the vertical axis, with elements at and n. If at=NULL (the default), the program attempts to find `nice' locations for the ticks, and the value of n (default, 5) gives the approximate number of tick marks desired; if at is non-NULL, then the value of n is ignored."

Thanks again.  Effects is a terrific package.
Paul

**************  Model Specification ****************
Clean.label <- polr(Clean.lbl ~ City.Abbr, method="logistic", data=Safeway, 
  Hess=TRUE)
List of 18
 $ coefficients : Named num [1:8] -1.0887 -1.1449 0.6923 0.0894 -0.8229 ...
  ..- attr(*, "names")= chr [1:8] "City.Abbr[T.Dublin]" "City.Abbr[T.Englwd]" "City.Abbr[T.Fairfax]" "City.Abbr[T.Falls Ch]" ...
 $ zeta         : Named num [1:4] -2.529 0.894 3.447 6.406
  ..- attr(*, "names")= chr [1:4] "Excellent|V.Good" "V.Good|Good" "Good|Fair" "Fair|Poor"
 $ deviance     : num 2327
 $ fitted.values: num [1:1248, 1:5] 0.0568 0.0568 0.0568 0.0568 0.0568 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:1248] "1" "2" "3" "4" ...
  .. ..$ : chr [1:5] "Excellent" "V.Good" "Good" "Fair" ...
 $ lev          : chr [1:5] "Excellent" "V.Good" "Good" "Fair" ...
 $ terms        :Classes 'terms', 'formula' length 3 Clean.lbl ~ City.Abbr
  .. ..- attr(*, "variables")= language list(Clean.lbl, City.Abbr)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "Clean.lbl" "City.Abbr"
  .. .. .. ..$ : chr "City.Abbr"
  .. ..- attr(*, "term.labels")= chr "City.Abbr"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(Clean.lbl, City.Abbr)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "ordered" "factor"
  .. .. ..- attr(*, "names")= chr [1:2] "Clean.lbl" "City.Abbr"
 $ df.residual  : num 1236
 $ edf          : num 12
 $ n            : num 1248
 $ nobs         : num 1248
 $ call         : language polr(formula = Clean.lbl ~ City.Abbr, data = Safeway, Hess = TRUE,      method = "logistic")
 $ method       : chr "logistic"
 $ convergence  : int 0
 $ niter        : Named int [1:2] 62 22
  ..- attr(*, "names")= chr [1:2] "f.evals.function" "g.evals.gradient"
 $ Hessian      : num [1:12, 1:12] 31.9 0 0 0 0 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:12] "City.Abbr[T.Dublin]" "City.Abbr[T.Englwd]" "City.Abbr[T.Fairfax]" "City.Abbr[T.Falls Ch]" ...
  .. ..$ : chr [1:12] "City.Abbr[T.Dublin]" "City.Abbr[T.Englwd]" "City.Abbr[T.Fairfax]" "City.Abbr[T.Falls Ch]" ...
 $ model        :'data.frame':	1248 obs. of  2 variables:
  ..$ Clean.lbl: Ord.factor w/ 5 levels "Excellent"<"V.Good"<..: 1 2 3 3 3 3 3 3 2 2 ...
  ..$ City.Abbr: Factor w/ 9 levels "DC","Dublin",..: 9 9 9 9 9 9 9 9 9 9 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Clean.lbl ~ City.Abbr
  .. .. ..- attr(*, "variables")= language list(Clean.lbl, City.Abbr)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "Clean.lbl" "City.Abbr"
  .. .. .. .. ..$ : chr "City.Abbr"
  .. .. ..- attr(*, "term.labels")= chr "City.Abbr"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(Clean.lbl, City.Abbr)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "ordered" "factor"
  .. .. .. ..- attr(*, "names")= chr [1:2] "Clean.lbl" "City.Abbr"
 $ contrasts    :List of 1
  ..$ City.Abbr: chr "contr.Treatment"
 $ xlevels      :List of 1
  ..$ City.Abbr: chr [1:9] "DC" "Dublin" "Englwd" "Fairfax" ...
 - attr(*, "class")= chr "polr"

Paul Prew  |  Statistician
651-795-5942?? |?? fax 651-204-7504 
Ecolab Research Center  | Mail Stop ESC-F4412-A 
655 Lone Oak Drive  |  Eagan, MN 55121-1560 


-----Original Message-----
From: John Fox [mailto:jfox at mcmaster.ca] 
Sent: Tuesday, April 28, 2009 10:34 AM
To: 'David Winsemius'
Cc: r-help at r-project.org; Prew, Paul
Subject: RE: [R] effects package --- add abline to plot

Dear David,
The plot() methods in the effects package make lattice graphs which in most instances will have more than one panel. For a binomial GLM, the default is to plot on the scale of the linear predictor (e.g., the logit scale) but to label the response axis on the scale of the response (i.e., the probability scale). To draw a line on the graph, even if you could do it, would require that you translate to the scale of the linear predictor [e.g., for a logit model, log(p/(1 - p))].
#
Dear Paul,
You need to specify ticks as a *list*; in your case, you just need the first element, so ticks=list(at=c(0.1,0.2,0.3,0.4,0.5,0.6)) or more compactly ticks=list(at=0.1*1:6) should do the trick. 

Can you tell me a bit more about why you want to draw horizontal lines on the plot? If this seems like a generally useful idea, I can put it my to-do list.

I hope this helps,
 John