[R] xyplot: key inside the plot region / lme: confidence bands for predicted
Kingsford Jones
kingsfordjones at gmail.com
Thu Jul 8 08:02:58 CEST 2010
On Wed, Jul 7, 2010 at 8:59 AM, Michael Friendly <friendly at yorku.ca> wrote:
> No one replied to my second question: how to get standard errors or
> confidence intervals for the
> estimated fixed effects from lme(). AFAICS, intervals() only gives CIs
> for coefficients.
Assuming you mean 'predictions based on the estimated fixed effects'
rather than 'estimated fixed effects', for one solution have a look at
Ben Bolker's r-sig-mixed-models FAQ, about 2/3 of the way down:
http://glmm.wikidot.com/faq
Googling "lme prediction intervals" turns up quite a bit. Here's a
recent conversation:
<http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/3291>
hope it helps,
Kingsford Jones
>
> My working example is:
>
> library(nlme)
> library(lattice)
>
> Ortho <- Orthodont
> Ortho$year <- Ortho$age - 8 # make intercept = initial status
>
> Ortho.mix1 <- lme(distance ~ year * Sex, data=Ortho,
> random = ~ 1 + year | Subject, method="ML")
> anova(Ortho.mix1)
>
> # get predicted values
> grid <- expand.grid(year=0:6, Sex=c("Male", "Female"))
> grid$age <- grid$year+8 # plot vs. age
> fm.mix1 <-cbind(grid, distance = predict(Ortho.mix1, newdata = grid,
> level=0))
>
> xyplot(distance ~ age, data=fm.mix1, groups=Sex, type="b",
> par.settings = list(superpose.symbol = list(cex = 1.2, pch=c(15,16))),
> auto.key=list(text=levels(fm.mix1$Sex), points = TRUE, x=0.05, y=0.9,
> corner=c(0,1)),
> main="Linear mixed model: predicted growth")
>
>> intervals(Ortho.mix1)
> Approximate 95% confidence intervals
>
> Fixed effects:
> lower est. upper
> (Intercept) 21.607212 22.615625 23.6240383
> year 0.619660 0.784375 0.9490904
> SexFemale -3.041252 -1.406534 0.2281834
> year:SexFemale -0.562889 -0.304830 -0.0467701
> attr(,"label")
> [1] "Fixed effects:"
>
> Random Effects:
> Level: Subject
> lower est. upper
> sd((Intercept)) 1.0710615 1.7045122 2.712600
> sd(year) 0.0281228 0.1541351 0.844783
> cor((Intercept),year) -0.9358505 -0.0311195 0.927652
>
> Within-group standard error:
> lower est. upper
> 1.07081 1.31004 1.60272
>>
>
>
> Peter Ehlers wrote:
>>
>> On 2010-07-02 9:37, Michael Friendly wrote:
>>>
>>> I have two questions related to plotting predicted values for a linear
>>> mixed model using xyplot:
>>>
>>> 1: With a groups= argument, I can't seem to get the key to appear inside
>>> the xyplot. (I have the Lattice book,
>>> but don't find an example that actually does this.)
>>> 2: With lme(), how can I generate confidence bands or prediction
>>> intervals around the fitted values? Once
>>> I get them, I'd like to add them to the plot.
>>> Example:
>>>
>>> library(nlme)
>>> library(lattice)
>>> Ortho <- Orthodont
>>> Ortho$year <- Ortho$year - 8 # make intercept = initial status
>>>
>>> Ortho.mix1 <- lme(distance ~ year * Sex, data=Ortho,
>>> random = ~ 1 + year | Subject, method="ML",
>>> # correlation = corAR1 (form = ~ 1 | Subject)
>>> )
>>> Ortho.mix1
>>>
>>> # predicted values
>>> grid <- expand.grid(year=0:6, Sex=c("Male", "Female"))
>>> grid$age <- grid$year+8 # plot vs. age
>>> fm.mix1 <-cbind(grid, distance = predict(Ortho.mix1, newdata = grid,
>>> level=0))
>>>
>>> xyplot(distance ~ age, data=fm.mix1, groups=Sex, type="b", pch=c(15,16),
>>> cex=1.2,
>>> auto.key=list(text=c("Male", "Female"), points = TRUE, x=8, y=26),
>>> main="Linear mixed model: predicted growth")
>>>
>>> Can someone help?
>>>
>>> -Michael
>>
>>
>> Michael,
>>
>> the x,y location should be specified in the unit square
>> (and you might want to set the 'corner' component).
>>
>> Ortho$year <- Ortho$age - 8 ##fix typo
>> xyplot(distance ~ age, data=fm.mix1, groups=Sex,
>> type="b", pch=c(15,16), cex=1.2,
>> auto.key=list(text=c("Male", "Female"),
>> points = TRUE,
>> x=0.1, y=0.8, corner=c(0,1)))
>>
>> See description of 'key' in ?xyplot.
>>
>> -Peter Ehlers
>
>
> --
> Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology
> Dept.
> York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Street Web: http://www.datavis.ca
> Toronto, ONT M3J 1P3 CANADA
>
> ______________________________________________
> 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.
>
More information about the R-help
mailing list