if (FALSE) {
library(agridat)
data(theobald.covariate)
dat <- theobald.covariate
libs(lattice)
xyplot(yield ~ chu|gen, dat, type=c('p','smooth'),
xlab = "chu = corn heat units",
main="theobald.covariate - yield vs heat")
# REML estimates (Means) in table 3 of Theobald 2002
libs(lme4)
dat <- transform(dat, year=factor(year))
m0 <- lmer(yield ~ -1 + gen + (1|year/env) + (1|gen:year), data=dat)
round(fixef(m0),2)
# Use JAGS to fit Theobald (2002) model 3.2 with 'Expert' prior
# Requires JAGS to be installed
if(0) {
libs(reshape2)
ymat <- acast(dat, year+env~gen, value.var='yield')
chu <- acast(dat, year+env~., mean, value.var='chu', na.rm=TRUE)
chu <- as.vector(chu - mean(chu)) # Center the covariate
dat$yr <- as.numeric(dat$year)
yridx <- as.vector(acast(dat, year+env~., mean, value.var='yr', na.rm=TRUE))
dat$loc <- as.numeric(dat$env)
locidx <- acast(dat, year+env~., mean, value.var='loc', na.rm=TRUE)
locidx <- as.vector(locidx)
jdat <- list(nVar = 10, nYear = 5, nLoc = 7, nYL = 29, yield = ymat,
chu = chu, year = yridx, loc = locidx)
libs(rjags)
m1 <- jags.model(file=system.file(package="agridat", "files/theobald.covariate.jag"),
data=jdat, n.chains=2)
# Table 3, Variety deviations from means (Expert prior)
c1 <- coda.samples(m1, variable.names=(c('alpha')),
n.iter=10000, thin=10)
s1 <- summary(c1)
effs <- s1$statistics[,'Mean']
# Perfect match (different order?)
rev(sort(round(effs - mean(effs), 2)))
}
}
Run the code above in your browser using DataLab