# NOT RUN {
data(kempton.barley.uniformity)
dat <- kempton.barley.uniformity
if(require(desplot)){
desplot(yield~col*row, dat,
aspect=140/98, tick=TRUE, # true aspect
main="kempton.barley.uniformity")
}
# ----------------------------------------------------------------------------
# }
# NOT RUN {
# Kempton estimated auto-regression coefficients b1=0.10, b2=0.91
dat <- transform(dat, xf = factor(col), yf=factor(row))
# asreml3
require(asreml)
m1 <- asreml(yield ~ 1, data=dat, rcov=~ar1(xf):ar1(yf))
require(lucid)
vc(m1)
## effect component std.error z.ratio constr
## R!variance 0.1044 0.022 4.7 pos
## R!xf.cor 0.2458 0.07484 3.3 uncon
## R!yf.cor 0.8186 0.03822 21 uncon
# asreml estimates auto-regression correlations of 0.25, 0.82
# Kempton estimated auto-regression coefficients b1=0.10, b2=0.91
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# }
# NOT RUN {
# Kempton estimated auto-regression coefficients b1=0.10, b2=0.91
## dat <- transform(dat, xf = factor(col), yf=factor(row))
## require(asreml4)
## m1 <- asreml(yield ~ 1, data=dat, residual = ~ar1(xf):ar1(yf))
## require(lucid)
## vc(m1)
## ## effect component std.error z.ratio bound <!-- %ch -->
## ## xf:yf(R) 0.1044 0.022 4.7 P 0
## ## xf:yf!xf!cor 0.2458 0.07484 3.3 U 0
## ## xf:yf!yf!cor 0.8186 0.03822 21 U 0
## # asreml estimates auto-regression correlations of 0.25, 0.82
## # Kempton estimated auto-regression coefficients b1=0.10, b2=0.91
# }
# NOT RUN {
# ----------------------------------------------------------------------------
# Kempton defines 4 blocks, randomly assigns variety codes 1-49 in each block, fits
# RCB model, computes mean squares for variety and residual. Repeat 40 times.
# Kempton's estimate: variety = 1032, residual = 1013
# Our estimate: variety = 825, residual = 1080
fitfun <- function(dat){
dat <- transform(dat, block=factor(ceiling(row/7)),
gen=factor(c(sample(1:49),sample(1:49),sample(1:49),sample(1:49))))
m2 <- lm(yield*100 ~ block + gen, dat)
anova(m2)[2:3,'Mean Sq']
}
set.seed(251)
out <- replicate(50, fitfun(dat))
rowMeans(out) # 826 1079
# ----------------------------------------------------------------------------
# }
Run the code above in your browser using DataCamp Workspace