John,
 
Thank you for the code, it looks like just what's needed.  I'll give it a try tomorrow.  
 
The desire for reference lines is to aid comparisons of different factor levels, to look at whether confidence intervals overlap.  I've attached a graphic file of the particular allEffects plot, if you're interested in the types of comparisons.  
 
A prediction: the client is going to lay down a ruler and examine which levels don't overlap.  I'd like to save them some trouble.
 
Thanks again, Paul

________________________________

From: John Fox [mailto:jfox@mcmaster.ca]
Sent: Tue 4/28/2009 2:48 PM
To: Prew, Paul
Cc: r-help@r-project.org; 'David Winsemius'
Subject: RE: [R] effects package --- add abline to plot



Dear Paul,

> -----Original Message-----
> From: Prew, Paul [mailto:Paul.Prew@ecolab.com]
> Sent: April-28-09 2:19 PM
> To: John Fox; David Winsemius
> Cc: r-help@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@mcmaster.ca]
> Sent: Tuesday, April 28, 2009 10:34 AM
> To: 'David Winsemius'
> Cc: r-help@r-project.org; Prew, Paul
> Subject: RE: [R] effects package --- add abline to plot
>
> Dear David,
>
> > -----Original Message-----
> > From: r-help-bounces@r-project.org [mailto:r-help-bounces@r-project.org] On
> > Behalf Of David Winsemius
> > Sent: April-28-09 10:12 AM
> > To: Prew, Paul
> > Cc: r-help@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@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.






CONFIDENTIALITY NOTICE: \ This e-mail communication an...{{dropped:11}}


