dat <- ilri.sheep
dat <- transform(dat, lamb=factor(lamb), ewe=factor(ewe), ram=factor(ram),
year=factor(year))
# dl is linear covariate, same as damage, but truncated to [2,8]
dat <- within(dat, {
dl <- damage
dl <- ifelse(dl < 3, 2, dl)
dl <- ifelse(dl > 7, 8, dl)
dq <- dl^2
})
dat <- subset(dat, !is.na(weanage))
# EDA
require("lattice")
bwplot(weanwt ~ year, dat, main="ilri.sheep") # Year effect
bwplot(weanwt ~ factor(dl), dat) # Dam age effect
bwplot(weanwt ~ gen, dat) # Genotype differences
xyplot(weanwt ~ weanage, dat, type=c('p','smooth')) # Age covariate
# case study page 4.18
lm1 <- lm(weanwt ~ year + sex + weanage + dl + dq + ewegen + ramgen, data=dat)
summary(lm1)
anova(lm1)
require("lme4")
lme1 <- lmer(weanwt ~ year + sex + weanage + dl + dq + ewegen + ramgen +
(1|ewe) + (1|ram), data=dat)
print(lme1, corr=FALSE)
lme2 <- lmer(weanwt ~ year + sex + weanage + dl + dq + ewegen + ramgen +
(1|ewe), data=dat)
lme3 <- lmer(weanwt ~ year + sex + weanage + dl + dq + ewegen + ramgen +
(1|ram), data=dat)
anova(lme1, lme2, lme3)
require("asreml")
# case study page 4.20
m1 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen,
data=dat)
anova(m1)
# case study page 4.26
m2 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen,
random = ~ ram + ewe, data=dat)
anova(m2)
# case study page 4.37, year means
predict(m2, classify="year")$predictionsRun the code above in your browser using DataLab