[R] effects package --- add abline to plot
John Fox
jfox at mcmaster.ca
Tue Apr 28 21:48:09 CEST 2009
Dear Paul,
> -----Original Message-----
> From: Prew, Paul [mailto:Paul.Prew at ecolab.com]
> Sent: April-28-09 2:19 PM
> To: John Fox; David Winsemius
> Cc: r-help at r-project.org
> Subject: RE: [R] effects package --- add abline to plot
>
> 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."
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
> Thanks again. Effects is a terrific package.
> Paul
>
> ************** Model Specification ****************
> Clean.label <- polr(Clean.lbl ~ City.Abbr, method="logistic", data=Safeway,
> Hess=TRUE)
>
> > str(Clean.label)
> 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,
>
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
> > Behalf Of David Winsemius
> > Sent: April-28-09 10:12 AM
> > To: Prew, Paul
> > Cc: r-help at r-project.org
> > Subject: Re: [R] effects package --- add abline to plot
> >
> >
> > On Apr 28, 2009, at 12:00 AM, Prew, Paul wrote:
> >
> > > Hello, I am not having success in a simple task. Using the effects
> > > package, I would like to add reference lines at probability values
> > > of 0.1 – 0.6 on a plot of the effects.
> >
> > I have concerns that you are considering these probabilities. They are
> > not going to be probabilities. They are effects.
> >
> > > The plot command works, but following up with an abline command
> > > produces the message “plot .new has not been called yet”, and of
> > > course the reference lines were not added.
> > >
> > > Looking through past R help lists, there was a similar request for
> > > help --- trying to add an abline but “got the error "plot.new has
> > > not been called yet".
> > >
> > > The help list reply was
> > >
> > > “ ?abline: "This function adds one or more straight lines through
> > > the
> > > current plot.", i.e. the already existing *current plot*.
> > >
> > > So plot your data (e.g. with plot(x, y)) before adding a regression
> > > line.”
> > >
> > > I interpreted the above to suggest the following ---
> > >
> > > plot(allEffects(Clean.label),ask=FALSE, alternating = TRUE,
> > > ylab="Probability of Rating", xlab="City",main="Cleanliness Ratings
> > > by City",
> > > factor.names=FALSE, ticks=c(0.1,0.2,0.3,0.4,0.5,0.6))
> > > abline(h=c(0.1,0.2,0.3,0.4,0.5,0.6))
> >
> > 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
>
> 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))].
>
> > >
> > >
> > > Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) :
> > > plot.new has not been called yet
> > >
> > > Less bothersome is the fact that the tick marks weren’t modified to
> > > 0.1, 0.2, 0.3, etc.
>
> From ?plot.eff
>
> "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."
>
> > >
> > > Further searching brought the panel.abline command to light, but
> > > that didn’t produce any results, not even an error message.
> > >
> > >> plot(allEffects(Clean.label),ask=FALSE, alternating = TRUE,
> > > + ylab="Probability of Rating", xlab="City",main="Cleanliness
> > > Ratings by City",
> > > + factor.names=FALSE, ticks=c(0.1,0.2,0.3,0.4,0.5,0.6))
> > >
> > >> panel.abline(h=c(0.1,0.2,0.3,0.4,0.5,0.6))
>
> 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
>
> >
> > >>
> > >
> > >
> >
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> > ;;;;;;;;;;;;;;
> > >> sessionInfo()
> > > R version 2.9.0 RC (2009-04-10 r48318)
> > > i386-pc-mingw32
> > >
> > > locale:
> > > LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.
> > > 1252;LC_MONETARY=English_United States.
> > > 1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
> > >
> > > attached base packages:
> > > [1] tcltk grid stats graphics grDevices utils
> > > datasets methods
> > > [9] base
> > >
> > > other attached packages:
> > > [1] relimp_1.0-1 Rcmdr_1.4-9 car_1.2-13 effects_2.0-4
> > > [5] colorspace_1.0-0 nnet_7.2-46 MASS_7.2-46 lattice_0.17-22
> > >
> > > loaded via a namespace (and not attached):
> > > [1] tools_2.9.0
> > > Thank you for any advice.
> > > Paul
> > >
> > > 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
> > >
> > >
> >
> > David Winsemius, MD
> > Heritage Laboratories
> > West Hartford, CT
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>
>
> CONFIDENTIALITY NOTICE:
> This e-mail communication and any attachments may contain proprietary and
> privileged information for the use of the designated recipients named above.
> Any unauthorized review, use, disclosure or distribution is prohibited.
> If you are not the intended recipient, please contact the sender by reply e-
> mail and destroy all copies of the original message.
More information about the R-help
mailing list