[R] How do I overlay two trellis plots of lme fitted lines produced by plot.augPred?
Ramzi Nahhas
ramzi.nahhas at wright.edu
Thu Jul 7 19:09:19 CEST 2011
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
More information about the R-help
mailing list