# NOT RUN {
library(agridat)
data(butron.maize)
dat <- butron.maize
libs(reshape2)
mat <- acast(dat, gen~env, value.var='yield')
mat <- sweep(mat, 2, colMeans(mat))
mat.svd <- svd(mat)
# Calculate PC1 and PC2 scores as in Table 4 of Butron
# Comment out to keep Rcmd check from choking on '<!-- %*%' -->
# round(mat.svd$u[,1:2] <!-- %*% diag(sqrt(mat.svd$d[1:2])) %*% diag(c(-1,1)),3) -->
biplot(princomp(mat), main="butron.maize", cex=.7) # Figure 1 of Butron
if(0){
# Fixme: This section is broken because synbreed has been removed from CRAN
# Here we see if including pedigree information is helpful for a
# multi-environment model
# Including the pedigree provided little benefit
# Create the pedigree
ped <- dat[, c('gen','male','female')]
ped <- ped[!duplicated(ped),] # remove duplicates
unip <- unique(c(ped$male, ped$female)) # Unique parents
unip <- unip[!is.na(unip)]
# We have to define parents at the TOP of the pedigree
ped <- rbind(data.frame(gen=c("Dent","Flint"), # genetic groups
male=c(0,0),
female=c(0,0)),
data.frame(gen=c("A509","A637","A661","CM105","EP28",
"EP31","EP42","F7","PB60","Z77016"),
male=rep(c('Dent','Flint'),each=5),
female=rep(c('Dent','Flint'),each=5)),
ped)
ped[is.na(ped$male),'male'] <- 0
ped[is.na(ped$female),'female'] <- 0
}
# asreml 3
if(0){
libs(asreml)
if( utils::packageVersion("asreml") < "4") {
# asreml3
ped.ainv <- asreml.Ainverse(ped)$ginv
m0 <- asreml(yield ~ 1+env, random = ~ gen, data=dat)
m1 <- asreml(yield ~ 1+env, random = ~ ped(gen),
ginverse=list(gen=ped.ainv), data=dat)
m2 <- update(m1, random = ~ id(env):ped(gen))
m3 <- update(m2, random = ~ diag(env):ped(gen))
m4 <- update(m3, random = ~ fa(env,1):ped(gen))
## AIC(m0,m1,m2,m3,m4)
## df AIC
## m0 2 229.4037
## m1 2 213.2487
## m2 2 290.6156
## m3 6 296.8061
## m4 11 218.1568
p0 <- predict(m0, data=dat, classify="gen")$pred$pvals
p4 <- predict(m4, data=dat, classify="gen")$pred$pvals
p4par <- p4[1:12,] # parents
p4 <- p4[-c(1:12),] # hybrids
# Careful! Need to manually sort the predictions
p0 <- p0[order(as.character(p0$gen)),]
p4 <- p4[order(as.character(p4$gen)),]
# lims <- range(c(p0$pred, p4$pred)) * c(.95,1.05)
lims <- c(6,8.25) # zoom in on the higher-yielding hybrids
plot(p0$predicted.value, p4$predicted.value,
pch="", xlim=lims, ylim=lims, main="butron.maize",
xlab="BLUP w/o pedigree", ylab="BLUP with pedigree")
abline(0,1,col="lightgray")
text(x=p0$predicted.value, y=p4$predicted.value, p0$gen, cex=.5, srt=-45)
text(x=min(lims), y=p4par$predicted.value, p4par$gen, cex=.5)
}
# asreml 4
libs(asreml)
if( utils::packageVersion("asreml") > "4") {
# asreml4
ped.ainv <- ainverse(ped)
m0 <- asreml(yield ~ 1+env, data=dat, random = ~ gen)
m1 <- asreml(yield ~ 1+env, random = ~ vm(gen, ped.ainv), data=dat)
m2 <- update(m1, random = ~ idv(env):vm(gen, ped.ainv))
m3 <- update(m2, random = ~ diag(env):vm(gen, ped.ainv))
m4 <- update(m3, random = ~ fa(env,1):vm(gen, ped.ainv))
#summary(m0)$aic
#summary(m4)$aic
## df AIC
## m0 2 229.4037
## m1 2 213.2487
## m2 2 290.6156
## m3 6 296.8061
## m4 11 218.1568
p0 <- predict(m0, data=dat, classify="gen")$pvals
p1 <- predict(m1, data=dat, classify="gen")$pvals
p1par <- p1[1:12,] # parents
p1 <- p1[-c(1:12),] # remove parents
# Careful! Need to manually sort the predictions
p0 <- p0[order(as.character(p0$gen)),]
p1 <- p1[order(as.character(p1$gen)),]
# lims <- range(c(p0$pred, p1$pred)) * c(.95,1.05)
lims <- c(6,8.25) # zoom in on the higher-yielding hybrids
plot(p0$predicted.value, p1$predicted.value,
pch="", xlim=lims, ylim=lims, main="butron.maize",
xlab="BLUP w/o pedigree", ylab="BLUP with pedigree")
abline(0,1,col="lightgray")
text(x=p0$predicted.value, y=p1$predicted.value,
p0$gen, cex=.5, srt=-45)
text(x=min(lims), y=p1par$predicted.value, p1par$gen, cex=.5, col="red")
round( cor(p0$predicted.value, p1$predicted.value), 3)
# Including the pedigree provided very little change
}
}
# }
Run the code above in your browser using DataLab