[R] xyplot: key inside the plot region / lme: confidence bands for predicted

[Michael Friendly]

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.

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