# NOT RUN {
data(stroup.nin)
dat <- stroup.nin
# Experiment layout. All "Buckskin" plots are near left side
if(require(desplot)){
desplot(yield~col*row, dat,
aspect=47.3/26.4, out1="rep", num=gen, cex=0.6, # true aspect
main="stroup.nin - yield heatmap (true shape)")
}
# ----- nlme -----
require(nlme)
# Random RCB model
lme1 <- lme(yield ~ 0 + gen, random=~1|rep, data=dat, na.action=na.omit)
# Linear (Manhattan distance) correlation model
lme2 <- gls(yield ~ 0 + gen, data=dat,
correlation = corLin(form = ~ col + row, nugget=TRUE),
na.action=na.omit)
# Random block and spatial correlation.
# Note: corExp and corSpher give nearly identical results
lme3 <- lme(yield ~ 0 + gen, data=dat,
random = ~ 1 | rep,
correlation = corExp(form = ~ col + row),
na.action=na.omit)
AIC(lme1,lme2,lme3) # lme2 is lowest
## df AIC
## lme1 58 1333.702
## lme2 59 1189.135
## lme3 59 1216.704
# Compare the estimates from the two methods
eff = data.frame(ranblock=fixef(lme1), corLin = coef(lme2),
corExp = fixef(lme3))
rownames(eff) <- gsub("gen", "", rownames(eff))
plot(eff$ranblock, eff$corLin, xlim=c(13,37), ylim=c(13,37),
main="stroup.nin - Gen estimates (intercepts differ)",
xlab="RCB (random block)", ylab="corLin",type='n')
text(eff$ranblock, eff$corLin, rownames(eff), cex=0.5)
abline(0,1)
# ----------------------------------------------------------------------------
# Note, the 'sommer' package has dropped support for ar1
# ----------------------------------------------------------------------------
# }
# NOT RUN {
if(require(SpATS)){
dat <- transform(dat, yf = as.factor(row), xf = as.factor(col))
sp1 <- SpATS(response = "yield",
spatial = ~ SAP(col, row, nseg = c(10,20), degree = 3, pord = 2),
genotype = "gen",
#fixed = ~ colcode + rowcode,
random = ~ yf + xf,
data = dat,
control = list(tolerance = 1e-03))
plot(sp1)
eff <- cbind(eff, spats=predict(sp1, which="gen")$predicted.values)
plot(eff$ranblock, eff$spats, xlim=c(13,37), ylim=c(13,37),
main="stroup.nin - Gen estimates",
xlab="RCB (random block)", ylab="SpATS",type='n')
text(eff$ranblock, eff$spats, rownames(eff), cex=0.5)
abline(0,1)
}
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# }
# NOT RUN {
# asreml3
require(asreml)
# RCB analysis
dat.rcb <- asreml(yield ~ gen, random = ~ rep, data=dat,
na.method.X="omit")
pred.rcb <- predict(dat.rcb, data=dat, classify="gen")$predictions
# Two-dimensional AR1xAR1 spatial model
dat <- transform(dat, xf=factor(col), yf=factor(row))
dat <- dat[order(dat$xf, dat$yf),]
dat.sp <- asreml(yield~gen, data=dat,
rcov=~ar1(xf):ar1(yf),
na.method.X='ignore')
pred.sp <- predict(dat.sp, data=dat, classify="gen")$predictions
require(lucid)
vc(dat.sp)
## effect component std.error z.ratio constr
## R!variance 48.7 7.155 6.8 pos
## R!xf.cor 0.6555 0.05638 12 unc
## R!yf.cor 0.4375 0.0806 5.4 unc
# Compare the estimates from the two methods
plot(pred.rcb$pvals[,2], pred.sp$pvals[,2], xlim=c(13,37), ylim=c(13,37),
xlab="RCB",ylab="AR1xAR1",type='n')
title("stroup.nin: Comparison of predicted values")
text(pred.rcb$pvals[,2],pred.sp$pvals[,2],
as.character(pred.rcb$pvals[,1]),cex=0.5)
abline(0,1)
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# }
# NOT RUN {
## require(asreml4)
## data(stroup.nin)
## dat <- stroup.nin
## # Note, ASREML4 stand-alone might not need a completely-filled rectangle
## # of plots to perform AR1xAR1 analysis. It might be able to fill in missing
## # combinations of rows/columns. The R version asreml4 still seems to require
## # rectangles.
## # RCB analysis
## dat.rcb <- asreml(yield ~ gen, random = ~ rep, data=dat,
## na.action=na.method(x="omit"))
## pred.rcb <- predict(dat.rcb, data=dat, classify="gen")
## # Two-dimensional AR1xAR1 spatial model
## dat <- transform(dat, xf=factor(col), yf=factor(row))
## #dat <- dat[order(dat$xf, dat$yf),]
## dat.sp <- asreml(yield~gen, data=dat,
## residual = ~ar1(xf):ar1(yf),
## na.action=na.method(x="omit"))
## pred.sp <- predict(dat.sp, data=dat, classify="gen")
## require(lucid)
## vc(dat.sp)
## ## effect component std.error z.ratio constr
## ## R!variance 48.7 7.155 6.8 pos
## ## R!xf.cor 0.6555 0.05638 12 unc
## ## R!yf.cor 0.4375 0.0806 5.4 unc
## # Compare the estimates from the two methods
## plot(pred.rcb$pvals[,2], pred.sp$pvals[,2], xlim=c(13,37), ylim=c(13,37),
## xlab="RCB", ylab="AR1xAR1", type='n')
## title("stroup.nin: Comparison of predicted values")
## text(pred.rcb$pvals[,2],pred.sp$pvals[,2],
## as.character(pred.rcb$pvals[,1]),cex=0.5)
## abline(0,1)
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# }
Run the code above in your browser using DataLab