# Which are the correct sum of squares for repeated measures ? (balanced example)

Below are three ways to fit a MANOVA model in R and to extract ANOVA tables from the fitted model. The design is balanced. All methods give different results.

The data are 5 groups of 4 individuals with three repeated measures shown in figure below

and simulated with the following R code:

R code:

library(ez)
library(nlme)
library(lattice)

##################################
### BALANCED DESGIN ##############
##################################
groups <- c("A","B","C","D","E")
n.groups <- length(groups)
n.individuals <- rep(4,n.groups) # four individuals per group
id <- unlist(lapply(1:n.groups, function(i) paste0(groups[i],1:n.individuals[i]))) # individual id
n <- length(id)
set.seed(666) # diabolic seed
dat <- data.frame(id=id,
group=rep(groups, each=4),
x=round(rnorm(n,2),1),
y=round(rnorm(n,4),1),
z=round(rnorm(n,6),1))
ldat <- melt(dat, id=c("id","group")) # data in long format
xtabs( ~ group+variable, data=ldat) # check design
# plot data:
xyplot(value ~ variable | group, groups=id, data=ldat,
panel = function(...){
panel.grid(h=-1,v=-3)
panel.superpose(..., type = "o", lty=2, pch=16)
}
)

#### 1st method : car package ####
mfit <- lm( cbind(x,y,z)~group, data=dat )
idata <- data.frame(variable=c("x","y","z"))
aov.carII <- Anova(mfit, idata=idata, idesign=~variable, type="II")
aov.carIII <- Anova(mfit, idata=idata, idesign=~variable, type="III")

#### 2nd method : ez package ####
aov.ezI <- ezANOVA(data = ldat,
dv = value,
wid = id,
within = .(variable),
between = group,
type=1
)
aov.ezII <- ezANOVA(data = ldat,
dv = value,
wid = id,
within = .(variable),
between = group,
type=2
)
aov.ezIII <- ezANOVA(data = ldat,
dv = value,
wid = id,
within = .(variable),
between = group,
type=3
)

### 3rd method : generalized least-squares fitting  ###
gfit <- gls(value ~ group*variable, data=ldat,
correlation=corSymm(form= ~ 1 | id),
weights=varIdent(form = ~1 | variable))


Comparison of ANOVA tables:

With gls() there are two possible ANOVA tables:

>  two possible ANOVA tables with gls :
> anova(gfit, type="sequential")
Denom. DF: 45
numDF   F-value p-value
(Intercept)        1 1401.9971  <.0001
group              4    2.3793  0.0658
variable           2   79.5687  <.0001
group:variable     8    1.4759  0.1929
> anova(gfit, type="marginal") # differs from sequential except for interaction
Denom. DF: 45
numDF   F-value p-value
(Intercept)        1 23.527654  <.0001
group              4  1.383566  0.2548
variable           2 14.166517  <.0001
group:variable     8  1.475904  0.1929


The type I ANOVA with ez (type I is not available with car):

> #########################
> ### TYPE I STATISTICS ###
> #########################
> aov.ezI$ANOVA # close to sequential anova(gfit) Effect DFn DFd F p p<.05 ges 1 group 4 15 1.996615 1.467733e-01 0.1651121 2 variable 2 30 80.595661 8.605897e-13 * 0.7715479 3 group:variable 8 30 1.093567 3.945923e-01 0.1549055 > # car: no Type I available  The type II ANOVA: > ########################## > ### TYPE II STATISTICS ### > ########################## > aov.ezII$ANOVA # identical to Type I ez
Effect DFn DFd         F            p p<.05       ges
2          group   4  15  1.996615 1.467733e-01       0.1651121
3       variable   2  30 80.595661 8.605897e-13     * 0.7715479
4 group:variable   8  30  1.093567 3.945923e-01       0.1549055

> aov.carII # car is close to ez but differs except for the first factor

Type II Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df    Pr(>F)
(Intercept)     1   0.98072   763.01      1     15 2.810e-14 ***
group           4   0.34744     2.00      4     15    0.1468
variable        1   0.91386    74.26      2     14 3.519e-08 ***
group:variable  4   0.51812     1.31      8     30    0.2758
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The type III ANOVA :

> ##########################
> ### TYPE III STATISTICS ###
> ##########################

## Conclusion

Is it a "bug" with SAS, something wrong ? That depends whether the SAS developers have an interpretation of their Type II table but it seems to be a naive generalization of the classical linear model to more general linear models, whereas J. Fox and G. Monette provide in their slides a well intepretable general approach to Type II hypotheses (is it an old or a recent approach ? I'd like to know).

### About the difference between Type II and Type III with the car package:

In spite of the balanced design there is a difference because Type III tests depend on the specified contrasts. With orthogonal constrats there's no difference:

> mfit <- update(mfit, contrasts=list(group=contr.sum))
> Anova(mfit, idata=idata, idesign=~variable, type="III")

Type III Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df    Pr(>F)
(Intercept)     1   0.98072   763.01      1     15 2.810e-14 ***
group           4   0.34744     2.00      4     15    0.1468
variable        1   0.91386    74.26      2     14 3.519e-08 ***
group:variable  4   0.51812     1.31      8     30    0.2758
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Finally with this example we could also suspect that SAS provides the Pillai statistic for factor group and the Wald statistic for factor variable. Yet I don't really know what's exactly going on... I will edit my posts in case of further understanding.