if (FALSE) {
library(agridat)
data(diggle.cow)
dat <- diggle.cow
# Figure 1 of Verbyla 1999
libs(latticeExtra)
useOuterStrips(xyplot(weight ~ day|iron*infect, dat, group=animal,
type='b', cex=.5,
main="diggle.cow"))
# Scaling
dat <- transform(dat, time = (day-122)/10)
if(require("asreml", quietly=TRUE)) {
libs(asreml, latticeExtra)
## # Smooth for each animal. No treatment effects. Similar to SAS Output 38.6.9
m1 <- asreml(weight ~ 1 + lin(time) + animal + animal:lin(time), data=dat,
random = ~ animal:spl(time))
p1 <- predict(m1, data=dat, classify="animal:time",
design.points=list(time=seq(0,65.9, length=50)))
p1 <- p1$pvals
p1 <- merge(dat, p1, all=TRUE) # to get iron/infect merged in
foo1 <- xyplot(weight ~ day|iron*infect, dat, group=animal,
main="diggle.cow")
foo2 <- xyplot(predicted.value ~ day|iron*infect, p1, type='l', group=animal)
print(foo1+foo2)
}
}
Run the code above in your browser using DataLab