Dear Paul, > -----Original Message----- > From: Prew, Paul [mailto:paul.p...@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:j...@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-boun...@r-project.org [mailto:r-help-boun...@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. ______________________________________________ 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.