agridat (version 1.12)

besag.met: Multi-environment trial of corn laid out in incomplete-blocks

Description

Multi-environment trial of corn laid out in incomplete-blocks

Format

A data frame with 1152 observations on the following 7 variables.

county

County factory, 1-6

row

Row ordinate

col

Column ordinate

rep

Rep factor, 1-3

block

Incomplete block factor, 1-8

yield

Yield

gen

Genotype factor, 1-64

Details

Multi-environment trial of 64 corn hybrids in six counties in North Carolina. Each location had 3 replicates in in incomplete-block design.

Note: In the original data, each county had 6 missing plots. This data has rows for each missing plot that uses the same county/block/rep to fill-out the row, sets the genotype to G01, and sets the yield to missing. These missing values were added to the data so that asreml could more easily do AR1xAR1 analysis using rectangular fields.

Examples

# NOT RUN {
data(besag.met)
dat <- besag.met

desplot(yield ~ col*row|county, out1=rep, out2=block, dat,
        main="besag.met")

# }
# NOT RUN {
# Heteroskedastic variance model (separate variance for each variety)
# asreml takes 1 second, lme 73 seconds, SAS PROC MIXED 30 minutes

# Average reps
datm <- aggregate(yield ~ county + gen, data=dat, FUN=mean)

# asreml Using 'rcov' ALWAYS requires sorting the data
require(asreml)
datm <- datm[order(datm$gen),]
m1a <- asreml(yield ~ gen, data=datm,
              random = ~ county,
              rcov = ~ at(gen):units,
              predict=asreml:::predict.asreml(classify="gen"))

require(lucid)
vc(m1a)[1:7,]
##             effect component std.error z.ratio constr
##  county!county.var   1324       838.2      1.6    pos
##   gen_G01!variance     91.93     58.82     1.6    pos
##   gen_G02!variance    210.7     133.9      1.6    pos
##   gen_G03!variance     63.03     40.53     1.6    pos
##   gen_G04!variance    112.1      71.53     1.6    pos
##   gen_G05!variance     28.39     18.63     1.5    pos
##   gen_G06!variance    237.4     150.8      1.6    pos

# lme
require(nlme)
m1l <- lme(yield ~ -1 + gen, data=datm, random=~1|county,
               weights = varIdent(form=~ 1|gen))
m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = F))^2
##           G02    G03    G04    G05    G06    G07    G08
##  91.90 210.75  63.03 112.05  28.39 237.36  72.72  42.97
## ... etc ...


# We get the same results from asreml & lme
plot(m1a$gammas[-1],
     m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = F))^2)

# }
# NOT RUN {
# }
# NOT RUN {
# The following example shows how to construct a GxE biplot
# from the FA2 model.

dat <- besag.met
dat <- transform(dat, xf=factor(col), yf=factor(row))
dat <- dat[order(dat$county, dat$xf, dat$yf), ]

# First, AR1xAR1
m1 <- asreml(yield ~ county, data=dat,
              random = ~ gen:county,
              rcov = ~ at(county):ar1(xf):ar1(yf))
# Add FA1
m2 <- update(m1, random=~gen:fa(county,1))
# FA2
m3 <- update(m2, random=~gen:fa(county,2))

# Use the loadings to make a biplot
vars <- vc(m3)
psi <- vars[grepl(".var$", vars$effect), "component"]
la1 <- vars[grepl(".fa1$", vars$effect), "component"]
la2 <- vars[grepl(".fa2$", vars$effect), "component"]
mat <- as.matrix(data.frame(psi, la1, la2))
rot <- svd(mat[,-1])$v # rotation matrix
lam <- mat[,-1] <!-- %*% rot # Rotate the loadings -->
colnames(lam) <- c("load1", "load2")

co3 <- coef(m3)$random # Scores are the GxE coefficients
ix1 <- grepl("_Comp1$", rownames(co3))
ix2 <- grepl("_Comp2$", rownames(co3))
sco <- matrix(c(co3[ix1], co3[ix2]), ncol=2, byrow=FALSE)
sco <- sco <!-- %*% rot # Rotate the scores -->
dimnames(sco) <- list(levels(dat$gen) , c('load1','load2'))
rownames(lam) <- levels(dat$county)
sco[,1] <- -1 * sco[,1]
lam[,1] <- -1 * lam[,1]
biplot(sco, lam, cex=.5, main="FA2 coefficient biplot")
# G variance matrix
gvar <- lam <!-- %*% t(lam) + diag(mat[,1]) -->

# Now get predictions and make an ordinary biplot
p3 <- predict(m3, classify="county:gen")
p3 <- p3$pred$pval
bi3 <- gge(predicted.value ~ gen*county, data=p3, scale=FALSE)
windows()
# Very similar to the coefficient biplot
biplot(bi3, stand=FALSE, title="SVD biplot of FA2 predictions")
# }