if (FALSE) {
library(agridat)
data(stroup.nin)
dat <- stroup.nin
# Experiment layout. All "Buckskin" plots are near left side and suffer
# from poor fertility in two of the reps.
libs(desplot)
desplot(dat, yield~col*row,
aspect=47.3/26.4, out1="rep", num=gen, cex=0.6, # true aspect
main="stroup.nin - yield heatmap (true shape)")
# Dataframe to hold model predictions
preds <- data.frame(gen=levels(dat$gen))
# -----
# nlme
libs(nlme)
# Random RCB model
lme1 <- lme(yield ~ 0 + gen, random=~1|rep, data=dat, na.action=na.omit)
preds$lme1 <- fixef(lme1)
# Linear (Manhattan distance) correlation model
lme2 <- gls(yield ~ 0 + gen, data=dat,
correlation = corLin(form = ~ col + row, nugget=TRUE),
na.action=na.omit)
preds$lme2 <- coef(lme2)
# 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)
preds$lme3 <- fixef(lme3)
# AIC(lme1,lme2,lme3) # lme2 is lowest
## df AIC
## lme1 58 1333.702
## lme2 59 1189.135
## lme3 59 1216.704
# -----
# SpATS
libs(SpATS)
dat <- transform(dat, yf = as.factor(row), xf = as.factor(col))
# what are colcode and rowcode???
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)
preds$spats <- predict(sp1, which="gen")$predicted.value
# -----
# Template Model Builder
# See the ar1xar1 example:
# https://github.com/kaskr/adcomp/tree/master/TMB/inst/examples
# This example uses dpois() in the cpp file to model a Poisson response
# with separable AR1xAR1. I think this example could be used for the
# stroup.nin data, changing dpois() to something Normal.
# -----
if(require("asreml", quietly=TRUE)){
libs(asreml,lucid)
# RCB analysis
as1 <- asreml(yield ~ gen, random = ~ rep, data=dat,
na.action=na.method(x="omit"))
preds$asreml1 <- predict(as1, data=dat, classify="gen")$pvals$predicted.value
# Two-dimensional AR1xAR1 spatial model
dat <- transform(dat, xf=factor(col), yf=factor(row))
dat <- dat[order(dat$xf, dat$yf),]
as2 <- asreml(yield~gen, data=dat,
residual = ~ar1(xf):ar1(yf),
na.action=na.method(x="omit"))
preds$asreml2 <- predict(as2, data=dat, classify="gen")$pvals$predicted.value
lucid::vc(as2)
## 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 asreml models.
# We see that Buckskin has correctly been shifted upward by the spatial model
plot(preds$as1, preds$as2, xlim=c(13,37), ylim=c(13,37),
xlab="RCB", ylab="AR1xAR1", type='n')
title("stroup.nin: Comparison of predicted values")
text(preds$asreml1, preds$asreml2, preds$gen, cex=0.5)
abline(0,1)
}
# -----
# sommer
# Fixed gen, random row, col, 2D spline
libs(sommer)
dat <- stroup.nin
dat <- transform(dat, yf = as.factor(row), xf = as.factor(col))
so1 <- mmer(yield ~ 0+gen,
random = ~ vs(xf) + vs(yf) + spl2Db(row,col),
data=dat)
preds$so1 <- coef(so1)[,"Estimate"]
# spatPlot
# -----
# compare variety effects from different packages
lattice::splom(preds[,-1], main="stroup.nin")
}
Run the code above in your browser using DataLab