if (FALSE) {
library(agridat)
libs(dplyr,lattice, reshape2, stringr)
data(edwards.oats)
dat <- edwards.oats
dat$env <- paste0(dat$year,".",dat$loc)
dat$eid <- factor(dat$eid)
mat <- reshape2::acast(dat, env ~ gen,
fun.aggregate=mean, value.var="yield", na.rm=TRUE)
lattice::levelplot(mat, aspect="m",
main="edwards.oats",
xlab="environment", ylab="genotype",
scales=list(x=list(rot=90)))
# Calculate BLUEs of gen/env effects
m1 <- lm(yield ~ gen+eid, dat)
gg <- coef(m1)[2:80]
names(gg) <- stringr::str_replace(names(gg), "gen", "")
gg <- c(0,gg)
names(gg)[1] <- "ACStewart"
ee <- coef(m1)[81:113]
names(ee) <- stringr::str_replace(names(ee), "eid", "")
ee <- c(0,ee)
names(ee)[1] <- "1"
# Subtract gen/env coefs from yield values
dat2 <- dat
dat2$gencoef <- gg[match(dat2$gen, names(gg))]
dat2$envcoef <- ee[match(dat2$eid, names(ee))]
dat2 <- dplyr::mutate(dat2, y = yield - gencoef - envcoef)
# Calculate variance for each gen*env. Shape of the graph is vaguely
# similar to Fig 2 of Edwards et al (2006), who used a Bayesian model
dat2 <- group_by(dat2, gen, eid)
dat2sum <- summarize(dat2, stddev = sd(y))
bwplot(stddev ~ eid, dat2sum)
}
Run the code above in your browser using DataLab