if (FALSE) {
library(agridat)
data(federer.diagcheck)
dat <- federer.diagcheck
dat$check <- ifelse(dat$gen == "G121" | dat$gen=="G122", "C","N")
# Show the layout as in Federer 1998.
libs(desplot)
desplot(dat, yield ~ col*row,
text=gen, show.key=FALSE, # aspect unknown
shorten='no', col=check, cex=.8, col.text=c("yellow","gray"),
main="federer.diagcheck")
# Now reproduce the analysis of Federer 2003.
# Only to match SAS results
dat$row <- 16 - dat$row
dat <- dat[order(dat$col, dat$row), ]
# Add row / column polynomials to the data.
# The scaling factors sqrt() are arbitrary, but used to match SAS
nr <- length(unique(dat$row))
nc <- length(unique(dat$col))
rpoly <- poly(dat$row, degree=10) * sqrt(nc)
cpoly <- poly(dat$col, degree=10) * sqrt(nr)
dat <- transform(dat,
c1 = cpoly[,1], c2 = cpoly[,2], c3 = cpoly[,3],
c4 = cpoly[,4], c6 = cpoly[,6], c8 = cpoly[,8],
r1 = rpoly[,1], r2 = rpoly[,2], r3 = rpoly[,3],
r4 = rpoly[,4], r8 = rpoly[,8], r10 = rpoly[,10])
dat$trtn <- ifelse(dat$gen == "G121" | dat$gen=="G122", dat$gen, "G999")
dat$new <- ifelse(dat$gen == "G121" | dat$gen=="G122", "N", "Y")
dat <- transform(dat, trtn=factor(trtn), new=factor(new))
m1 <- lm(yield ~ c1 + c2 + c3 + c4 + c6 + c8
+ r1 + r2 + r4 + r8 + r10
+ c1:r1 + c2:r1 + c3:r1 + gen, data = dat)
# To get Type III SS use the following
# libs(car)
# car::Anova(m1, type=3) # Matches PROC GLM output
## Sum Sq Df F value Pr(>F)
## (Intercept) 538948 1 159.5804 3.103e-16 ***
## c1 13781 1 4.0806 0.0494940 *
## c2 51102 1 15.1312 0.0003354 ***
## c3 45735 1 13.5419 0.0006332 ***
## c4 24670 1 7.3048 0.0097349 **
## ...
# lmer
libs(lme4,lucid)
# "group" for all data
dat$one <- factor(rep(1, nrow(dat)))
# lmer with bobyqa (default)
m2b <- lmer(yield ~ trtn + (0 + r1 + r2 + r4 + r8 + r10 +
c1 + c2 + c3 + c4 + c6 +
c8 + r1:c1 + r1:c2 + r1:c3 || one) +
(1|new:gen),
data = dat,
control=lmerControl(check.nlev.gtr.1="ignore"))
vc(m2b)
## grp var1 var2 vcov sdcor
## new.gen (Intercept) 2869 53.57
## one r1:c3 5532 74.37
## one.1 r1:c2 58230 241.3
## one.2 r1:c1 128000 357.8
## one.3 c8 6456 80.35
## one.4 c6 1400 37.41
## one.5 c4 1792 42.33
## one.6 c3 2549 50.49
## one.7 c2 5942 77.08
## one.8 c1 0 0
## one.9 r10 1133 33.66
## one.10 r8 1355 36.81
## one.11 r4 2269 47.63
## one.12 r2 241.8 15.55
## one.13 r1 9200 95.92
## Residual 4412 66.42
# lmer with Nelder_Mead gives 'wrong' results
## m2n <- lmer(yield ~ trtn + (0 + r1 + r2 + r4 + r8 + r10 +
## c1 + c2 + c3 + c4 + c6 + c8 + r1:c1 + r1:c2 + r1:c3 || one) +
## (1|new:gen)
## , data = dat,
## control=lmerControl(optimizer="Nelder_Mead",
## check.nlev.gtr.1="ignore"))
## vc(m2n)
## groups name variance stddev
## new.gen (Intercept) 3228 56.82
## one r1:c3 7688 87.68
## one.1 r1:c2 69750 264.1
## one.2 r1:c1 107400 327.8
## one.3 c8 6787 82.38
## one.4 c6 1636 40.45
## one.5 c4 12270 110.8
## one.6 c3 2686 51.83
## one.7 c2 7645 87.43
## one.8 c1 0 0.0351
## one.9 r10 1976 44.45
## one.10 r8 1241 35.23
## one.11 r4 2811 53.02
## one.12 r2 928.2 30.47
## one.13 r1 10360 101.8
## Residual 4127 64.24
if(require("asreml", quietly=TRUE)) {
libs(asreml,lucid)
m3 <- asreml(yield ~ -1 + trtn, data=dat,
random = ~ r1 + r2 + r4 + r8 + r10 +
c1 + c2 + c3 + c4 + c6 + c8 +
r1:c1 + r1:c2 + r1:c3 + new:gen)
## coef(m3)
## # REML cultivar means. Very similar to Federer table 2.
## rev(sort(round(coef(m3)$fixed[3] + coef(m3)$random[137:256,],0)))
## ## gen_G060 gen_G021 gen_G011 gen_G099 gen_G002
## ## 974 949 945 944 942
## ## gen_G118 gen_G058 gen_G035 gen_G111 gen_G120
## ## 938 937 937 933 932
## ## gen_G046 gen_G061 gen_G082 gen_G038 gen_G090
## ## 932 931 927 927 926
## vc(m3)
## ## effect component std.error z.ratio constr
## ## r1!r1.var 9201 13720 0.67 pos
## ## r2!r2.var 241.7 1059 0.23 pos
## ## r4!r4.var 2269 3915 0.58 pos
## ## r8!r8.var 1355 2627 0.52 pos
## ## r10!r10.var 1133 2312 0.49 pos
## ## c1!c1.var 0.01 0 4.8 bound
## ## c2!c2.var 5942 8969 0.66 pos
## ## c3!c3.var 2549 4177 0.61 pos
## ## c4!c4.var 1792 3106 0.58 pos
## ## c6!c6.var 1400 2551 0.55 pos
## ## c8!c8.var 6456 9702 0.67 pos
## ## r1:c1!r1.var 128000 189700 0.67 pos
## ## r1:c2!r1.var 58230 90820 0.64 pos
## ## r1:c3!r1.var 5531 16550 0.33 pos
## ## new:gen!new.var 2869 1367 2.1 pos
## ## R!variance 4412 915 4.8 pos
}
}
Run the code above in your browser using DataLab