[R] How do I overlay two trellis plots of lme fitted lines produced by plot.augPred?

[Dennis Murphy]

Hi:

Here's one way, using the latticeExtra package:

library(lattice)
library(latticeExtra)
p1 <- plot(augPred(       fit.lme               ), ylim=c(20,65))
p2 <- plot(augPred(update(fit.lme, y ~ log(age))), ylim=c(20,65),
col.line = 'red')
p1 + p2

HTH,
Dennis

On Thu, Jul 7, 2011 at 10:09 AM, Ramzi Nahhas <ramzi.nahhas at wright.edu> wrote:
> Hello,
>
> I want to use lme to fit two (or more) models, and then compare the fits on
> each individual. I know how to write my own code to do this (for each
> individual, plot the raw data, followed by lines() to plot each fitted
> curve) but I would like to use plot(augPred(... as it produces a nice
> trellis plot.  I thought I could do this with par(new=T) but it does not
> seem to work. trellis.par.set() does not have a component called "new" so I
> cannot use that to imitate what I would do with par.
>
> Does anyone know how to overlay two or more plots produced by plot.augPred?
>
> Thank you!
>
> Ramzi
>
> # Here is code that reproduces my problem:
> tmpdat =
> structure(list(ptno = c(4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 10L, 10L, 10L, 10L, 10L, 10L,
> 10L, 10L, 10L, 10L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,
> 20L, 20L, 20L, 20L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
> 23L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L), age =
> c(6.9897330595,
> 9.0075290897, 10.009582478, 12.005475702, 13.013004791, 13.990417522,
> 16.101300479, 17.062286105, 17.995893224, 20.996577687, 8.0027378508,
> 8.9965776865, 10.006844627, 15.085557837, 15.994524298, 16.993839836,
> 17.984941821, 22.431211499, 7.0034223135, 8, 8.9965776865, 11.047227926,
> 11.997262149, 14.124572211, 15.006160164, 15.986310746, 17.059548255,
> 18.001368925, 6.9952087611, 8.0027378508, 9.0075290897, 10.001368925,
> 11.000684463, 11.997262149, 12.988364134, 14.009582478, 14.986995209,
> 16.002737851, 16.963723477, 17.995893224, 22.702258727, 7.0034223135,
> 7.9972621492, 9.0047912389, 12.027378508, 12.982888433, 13.995893224,
> 15.282683094, 16.041067762, 17.029431896, 19.498973306, 7.0061601643,
> 8.9993155373, 10.015058179, 11.011635866, 12.999315537, 14.006844627,
> 15.030800821, 16.065708419, 17.015742642, 18.015058179), y = c(40.632301331,
> 40.996398926, 44.707000732, 39.288200378, 42.174301147, 50.073101044,
> 45.052200317, 44.862499237, 53.557498932, 54.667499542, 42.447601318,
> 43.082500458, 45.340499878, 47.740501404, 50.521499634, 52.097400665,
> 53.834999084, 53.002498627, 31.815000534, 41.812698364, 45.349998474,
> 42.350700378, 43.663898468, 43.941699982, 47.007400513, 46.071899414,
> 48.099998474, 50.319999695, 33.178501129, 38.366100311, 33.740398407,
> 40.996498108, 45.250499725, 40.964000702, 49.529201508, 46.635799408,
> 43.848800659, 49.038299561, 59.698799133, 57.349998474, 56.979999542,
> 37.90530014, 39.726600647, 34.828800201, 37.519298553, 36.681400299,
> 41.897899628, 48.111301422, 52.746299744, 45.144901276, 53.650001526,
> 42.632099152, 42.175498962, 46.064498901, 48.142799377, 51.856700897,
> 55.925800323, 51.002101898, 56.268901825, 61.367401123, 57.844799042
> )), .Names = c("ptno", "age", "y"), row.names = c(1L, 3L, 4L,
> 5L, 6L, 7L, 8L, 9L, 10L, 11L, 25L, 26L, 27L, 28L, 29L, 30L, 31L,
> 32L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 77L, 78L,
> 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L,
> 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L,
> 104L, 105L, 106L, 107L, 108L, 109L), class = "data.frame")
>
> library(nlme)
> grpdat = groupedData(y      ~ age | ptno,
>                     data   = tmpdat,
>                     FUN    = mean,
>                     labels = list(x="Age", y="Y"))
> fit.lme = lme(grpdat)
> plot(augPred(       fit.lme               ), ylim=c(20,65))
> plot(augPred(update(fit.lme, y ~ log(age))), ylim=c(20,65))
>
> # I want these two plots to be in the same window so I can visually compare
> the two fits
>
> # The following does not work:
> plot(augPred(       fit.lme               ), ylim=c(20,65),
> col.line="green")
> par(new=T)
> plot(augPred(update(fit.lme, y ~ log(age))), ylim=c(20,65), col.line="red")
> par(new=F)
>
> # trellis.par.set() does not have a component called "new" so I cannot use
> that to imitate what I would do with par
>
> --
> ---------------------------------------------------------------------
> - Ramzi W. Nahhas, Ph.D.
> - Assistant Professor
> - Lifespan Health Research Center, Department of Community Health
> - Boonshoft School of Medicine, Wright State University
> - 3171 Research Boulevard, Dayton OH 45420
> - Phone: (937) 775-1421   Cell: (937) 269-6797   Fax: (937) 775-1456
> - www.med.wright.edu/lhrc/fac/rn.html
>
> ______________________________________________
> 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.
>