# [S] Split-plot analysis

Laffont, Jean-Louis (Laffont@phibred.com)
Mon, 20 Apr 1998 07:53:46 -0500

Dear S-users,

I am using S-Plus 4.0, release 3 and I have some questions about
split-plot analysis.
I will use data frame oats provided in Venables & Ripley (p.178 - table
6.2) to illustrate my questions.

1. Standard error of differences when specifying the model in different
ways

Specifying the model as in Venables & Ripley works well for getting
standard error of differences:
> oats.aov1<-aov(Y~N*V+Error(B/V),data=oats,qr=T)
> se.contrast(oats.aov1,list(V=="Victory",V=="Marvellous"))
[,1]
[1,] 7.078904

However, specifying the model as below, I get an error message:
> oats.aov2<-aov(Y~B+N*V+Error(paste(B,V)+paste(B,V,N)),data=oats,qr=T)
> summary(oats.aov2)
Error: paste(B, V)
Df Sum of Sq Mean Sq F Value Pr(F)
B 5 15875.28 3175.056 5.28005 0.0124404
V 2 1786.36 893.181 1.48534 0.2723869
Residuals 10 6013.31 601.331

Error: paste(B, V, N)
Df Sum of Sq Mean Sq F Value Pr(F)
N 3 20020.50 6673.500 37.68565 0.0000000
N:V 6 321.75 53.625 0.30282 0.9321988
Residuals 45 7968.75 177.083
> se.contrast(oats.aov2,list(V=="Victory",V=="Marvellous"))
Error in c.qr[e.assign[[strata.nm]], , drop = F]: Array subscript (89)
out of bounds, should be at most 72
Dumped

Am I doing something wrong when calling se.contrast?

2. Anova table when missing data

Suppose I have 1 missing data in oats:

> oats.1NA<-oats[-1,]
> oats.1NA.aov1<-aov(Y~N*V+Error(B/V),data=oats.1NA,qr=T)
> summary(oats.1NA.aov1)
Error: B
Df Sum of Sq Mean Sq F Value Pr(F)
N 1 14758.53 14758.53 34.47047 0.004203461
Residuals 4 1712.60 428.15

Error: V %in% B
Df Sum of Sq Mean Sq F Value Pr(F)
N 1 1072.061 1072.061 2.097827 0.1814346
V 2 2847.696 1423.848 2.786211 0.1143401
Residuals 9 4599.304 511.034

Error: Within
Df Sum of Sq Mean Sq F Value Pr(F)
N 3 18732.82 6244.275 34.71880 0.0000000
N:V 6 299.32 49.886 0.27737 0.9446197
Residuals 44 7913.52 179.853

As I am not familiar with split-plot analysis with missing data, could
anyone confirm me this sum of squares partitioning for the whole-plot
analysis?

3. Analysis using lme

I tried to use lme for the analysis following "Complements to Modern
Applied Statistics with S+ - Venables & Ripley" that is:

> oats\$sp<-model.matrix(~V-1,oats)
> options(contrasts=c("contr.treatment","contr.poly"))
>
oats.lme<-lme(Y~N*V,random=~sp,cluster=~B,data=oats,re.block=list(1,2:4)
,
re.structure=c("unrestricted","identity"),
control=lme.control(tol=1e-10,ms.tol=1e-10))
> summary(oats.lme)
.....
Fixed Effects Estimate(s):
Value Approx. Std.Error z ratio(C)
(Intercept) 80.0000000 9.106976 8.7844743
N0.2cwt 18.5000000 7.682954 2.4079280
N0.4cwt 34.6666667 7.682954 4.5121533
N0.6cwt 44.8333333 7.682954 5.8354291
VMarvellous 6.6666667 9.715025 0.6862223
VVictory -8.5000000 9.715025 -0.8749334
......
The standard errors for variety differences don't agree with what I get
in §1 (that is 7.078904).
What am I doing wrong in writting the model and is it possible to have
some clarifications about the argument re.block?