dat <- diggle.cow
# Figure 1 of Verbyla 1999
require(lattice)
require(latticeExtra)
useOuterStrips(xyplot(weight ~ day|iron*infect, dat, type='b',group=animal,
main="diggle.cow"))
# Scaling
dat <- transform(dat, time = (day-122)/10)
# Smooth for each animal. No treatment effects. Similar to SAS Output 38.6.9
require(asreml)
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