[R] xyplot: adding pooled regression lines to a paneled type="r" plot

[Deepayan Sarkar]

On Wed, Jun 23, 2010 at 11:35 PM, Michael Friendly <friendly at yorku.ca> wrote:
> Thanks, Deepayan
>
> I read your presentation and understand how this works for the case you
> presented, but I can't
> get it to work for my case, where I want to superimpose model fitted lines
> over individual
> subject regression lines.  Here's what I tried
>
> library(nlme)
> library(lattice)
>
> ## ------------------
> ## pooled OLS model
> ## ------------------
>
> Ortho.OLS <- lm(distance ~ age * Sex, data=Orthodont)
> #coef(Ortho.OLS)
>
>   # plot individual lines
> xyplot(distance ~ age|Sex, data=Orthodont, type='r', groups=Subject,
> col=gray(.50),
>   main="Individual linear regressions ~ age")
>
> grid <- expand.grid(age=8:14, Sex=c("Male", "Female"))
>
> fm.OLS <-cbind(grid, distance = predict(Ortho.OLS, newdata = grid))
>
> Ortho <-Orthodont[c("age", "Sex", "distance")]
> combined <- make.groups(original = Ortho,
>                           OLS = fm.OLS)
> str(combined)
> xyplot(distance ~ age|Sex, data=combined, groups=which, col="black", lwd=2,
>  type = c("r", "l"), distribute.type = TRUE
>   )
>
> This last just gives me the pooled within-Sex regression lines, which is
> what I want to overlay
> on the first plot.

Your 'combined' dataset no longer has the "Subject" id-s, which you
need to be able to separate the per-subject lines. Try this instead:

grid <- expand.grid(age=8:14, Sex=c("Male", "Female"))

fm.OLS <-cbind(grid, distance = predict(Ortho.OLS, newdata = grid))
fm.OLS$Subject <- "Pooled"

combined <- rbind(Orthodont, fm.OLS[names(Orthodont)])

str(combined$Subject) ## last level is the "Pooled" `subject'

xyplot(distance ~ age | Sex, data=combined, groups=Subject,
       col=rep(c("grey80", "black"), c(nlevels(combined$Subject)-1, 1)),
       lwd=2, type = "r")


> Further, if I try a mixed model, I get errors trying to get the predicted
> values in a similar form
>
> Ortho.MLM <- lme(distance ~ age * Sex, data=Orthodont,
>       random = ~ 1 + age | Subject,
>       correlation = corAR1 (form = ~ 1 | Subject))
>
> fm.MLM <-cbind(grid, distance = predict(Ortho.MLM, newdata = grid))
>
>> fm.MLM <-cbind(grid, distance = predict(Ortho.MLM, newdata = grid))
> Error in predict.lme(Ortho.MLM, newdata = grid) :
>  Cannot evaluate groups for desired levels on "newdata"

Again, the subject id-s are missing. I think you want

predict(Ortho.MLM, newdata = grid, level = 0)

See ?predict.lme

-Deepayan

>
> Deepayan Sarkar wrote:
>>
>> On Tue, Jun 22, 2010 at 9:30 AM, Michael Friendly <friendly at yorku.ca>
>> wrote:
>>
>>>
>>> Consider the following plot that shows separate regression lines ~ age
>>> for
>>> each subject in the Pothoff-Roy Orthodont data,
>>> with separate panels by Sex:
>>>
>>> library(nlme)
>>> #plot(Orthodont)
>>> xyplot(distance ~ age|Sex, data=Orthodont, type='r', groups=Subject,
>>> col=gray(.50),
>>>  main="Individual linear regressions ~ age")
>>>
>>> I'd like to also show in each panel the pooled OLS regression line for
>>> each
>>> Sex in the corresponding panel,
>>> generated by the following model:
>>>
>>> Ortho.OLS <- lm(distance ~ age * Sex, data=Orthodont)
>>>
>>> Sex is a factor, with Male=0, so the coefficients are:
>>>
>>>>
>>>> coef(Ortho.OLS)
>>>>
>>>
>>>  (Intercept)           age     SexFemale age:SexFemale
>>>  16.3406250     0.7843750     1.0321023    -0.3048295
>>>
>>> I anticipate wanting to fit other models to these data, and also
>>> displaying
>>> the model-predicted
>>> regression lines in the same or similar plot, e.g., for a simple linear
>>> mixed model:
>>>
>>> Ortho.MLM <- lme(distance ~ age * Sex, data=Orthodont,
>>>      random = ~ 1 + age | Subject,
>>>      correlation = corAR1 (form = ~ 1 | Subject))
>>>
>>
>> Have a look at
>>
>> http://user2007.org/program/presentations/sarkar.pdf
>>
>> -Deepayan
>>
>
>
> --
> 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
>
>