data(diggle.cow)
dat <- diggle.cow
# Figure 1 of Verbyla 1999
require("lattice")
xyplot(weight ~ day|iron*infect, dat, type='b',group=animal)
# Scaling
dat <- transform(dat, time = (day-122)/10)
# Smooth for each animal. No treatment effects. Similar to Output 38.6.9 here:
# http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect018.htm
require("asreml")
require("latticeExtra")
m1 <- asreml(weight ~ 1 + lin(time) + animal + animal:lin(time), data=dat,
random = ~ animal:spl(time))
p1 <- predict(m1, classify="animal:time",
predictpoints=list(time=seq(0,65.9, length=50)))
p1 <- p1$pred$pval
p1 <- merge(dat, p1, all=TRUE) # to get iron/infect merged in
foo1 <- xyplot(weight ~ day|iron*infect, dat, group=animal)
foo2 <- xyplot(predicted.value ~ day|iron*infect, p1, type='l', group=animal)
print(foo1+foo2)Run the code above in your browser using DataLab