# NOT RUN {
data(gilmour.serpentine)
dat <- gilmour.serpentine
if(require(desplot)){
desplot(yield~ col*row, data=dat,
num=gen, show.key=FALSE, out1=rep,
aspect = 16.5/90, # true aspect
main="gilmour.serpentine")
}
# Extreme field trend. Blocking insufficient--needs a spline/smoother
# xyplot(yield~col, data=dat, main="gilmour.serpentine")
# ----------------------------------------------------------------------------
# }
# NOT RUN {
# asreml3
require(asreml)
dat <- transform(dat, rowf=factor(row), colf=factor(10*(col-8)))
dat <- dat[order(dat$rowf, dat$colf), ] # Sort order needed by asreml
# RCB
m0 <- asreml(yield ~ gen, data=dat, random=~rep)
# Add AR1 x AR1
m1 <- asreml(yield ~ gen, data=dat, rcov = ~ar1(rowf):ar1(colf))
# Add spline
m2 <- asreml(yield ~ gen + col, data=dat,
random= ~ spl(col) + colf,
rcov = ~ar1(rowf):ar1(colf))
# Figure 4 shows serpentine spraying
p2 <- predict(m2, data=dat, classify="colf")$predictions$pvals
plot(p2$predicted, type='b', xlab="column number", ylab="BLUP")
# Define column code (due to serpentine spraying)
# Rhelp doesn't like double-percent modulus symbol, so compute by hand
dat <- transform(dat, colcode = factor(dat$col-floor((dat$col-1)/4)*4 -1))
m3 <- asreml(yield ~ gen + lin(colf) + colcode, data=dat,
random= ~ colf + rowf + spl(colf),
rcov = ~ar1(rowf):ar1(colf))
# Figure 6 shows serpentine row effects
p3 <- predict(m3, data=dat, classify="rowf")$predictions$pvals
plot(p3$predicted, type='l', xlab="row number", ylab="BLUP")
text(1:22, p3$predicted, c('L','L','M','R','R','M','L','L',
'M','R','R','M','L','L','M','R','R','M','L','L','M','R'))
# Define row code (due to serpentine planting). 1=middle, 2=left/right
dat <- transform(dat, rowcode = factor(row))
levels(dat$rowcode) <- c('2','2','1','2','2','1','2','2','1',
'2','2','1','2','2','1','2','2','1','2','2','1','2')
m6 <- asreml(yield ~ gen + lin(colf) + colcode +rowcode, data=dat,
random= ~ colf + rowf + spl(col),
rcov = ~ar1(rowf):ar1(colf))
plot(variogram(m6), xlim=c(0:17), ylim=c(0,11), zlim=c(0,4000),
main="gilmour.serpentine")
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# }
# NOT RUN {
## require(asreml4)
## dat <- transform(dat, rowf=factor(row), colf=factor(10*(col-8)))
## dat <- dat[order(dat$rowf, dat$colf), ] # Sort order needed by asreml
## # RCB
## m0 <- asreml(yield ~ gen, data=dat, random=~rep)
## # Add AR1 x AR1
## m1 <- asreml(yield ~ gen, data=dat, resid = ~ar1(rowf):ar1(colf))
## # Add spline
## m2 <- asreml(yield ~ gen + col, data=dat,
## random= ~ spl(col) + colf,
## resid = ~ar1(rowf):ar1(colf))
## # Figure 4 shows serpentine spraying
## p2 <- predict(m2, data=dat, classify="colf")$pvals
## plot(p2$predicted, type='b', xlab="column number", ylab="BLUP")
## # Define column code (due to serpentine spraying)
## # Rhelp doesn't like double-percent modulus symbol, so compute by hand
## dat <- transform(dat, colcode = factor(dat$col-floor((dat$col-1)/4)*4 -1))
## m3 <- asreml(yield ~ gen + lin(colf) + colcode, data=dat,
## random= ~ colf + rowf + spl(colf),
## resid = ~ar1(rowf):ar1(colf))
## # Figure 6 shows serpentine row effects
## p3 <- predict(m3, data=dat, classify="rowf")$pvals
## plot(p3$predicted, type='l', xlab="row number", ylab="BLUP")
## text(1:22, p3$predicted, c('L','L','M','R','R','M','L','L',
## 'M','R','R','M','L','L','M','R','R','M','L','L','M','R'))
## # Define row code (due to serpentine planting). 1=middle, 2=left/right
## dat <- transform(dat, rowcode = factor(row))
## levels(dat$rowcode) <- c('2','2','1','2','2','1','2','2','1',
## '2','2','1','2','2','1','2','2','1','2','2','1','2')
## m6 <- asreml(yield ~ gen + lin(colf) + colcode +rowcode, data=dat,
## random= ~ colf + rowf + spl(col),
## resid = ~ar1(rowf):ar1(colf))
## plot(varioGram(m6), xlim=c(0:17), ylim=c(0,11), zlim=c(0,4000),
## main="gilmour.serpentine")
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab