# NOT RUN {
data(stroup.nin)
dat <- stroup.nin
dat <- transform(dat, xf=factor(x), yf=factor(y))
# Show experiment layout
# Note: all "Buckskin" plots are near left side
desplot(yield~x*y, dat, out1="rep", num=gen, cex=1, main="stroup.nin")
require(nlme)
# Random block 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, correlation = corLin(form = ~ x+y,nugget=TRUE), data=dat,
na.action=na.omit)
# Compare the estimates from the two methods
eff = data.frame(ranblock=fixef(lme1), spat = coef(lme2))
rownames(eff) <- gsub("gen", "", rownames(eff))
plot(eff$ranblock, eff$spat, xlim=c(13,37), ylim=c(13,37),
main="stroup.nin", xlab="RCB (random block)", ylab="corLin",type='n')
text(eff$ranblock, eff$spat, rownames(eff), cex=0.5)
abline(0,1)
# }
# NOT RUN {
require(asreml)
# RCB analysis
dat.rcb <- asreml(yield ~ gen, random = ~ rep, data=dat,
na.method.X="omit")
pred.rcb <- predict(dat.rcb,classify="gen")$predictions
# Two-dimensional AR1xAR1 spatial model
dat <- dat[order(dat$xf, dat$yf),]
dat.sp <- asreml(yield~gen, rcov=~ar1(xf):ar1(yf),data=dat)
pred.sp <- predict(dat.sp,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(16,37), ylim=c(16,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