[R] xyplot: adding pooled regression lines to a paneled type="r" plot
Deepayan Sarkar
deepayan.sarkar at gmail.com
Thu Jun 24 10:37:11 CEST 2010
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
>
>
More information about the R-help
mailing list