3
$\begingroup$

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

ldat dataset

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 ###
> ##########################
> aov.ezIII$ANOVA # identical to Type II except for the second factor
Effect DFn DFd         F            p p<.05       ges
2          group   4  15  1.996615 0.1467732928       0.1651121
3       variable   2  30 11.560705 0.0001895724     * 0.3263455
4 group:variable   8  30  1.093567 0.3945922933       0.1549055

> aov.carIII # identical to Type II except for the second factor and intercept

Type III Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df    Pr(>F)
(Intercept)     1   0.92734  191.436      1     15 6.038e-10 ***
group           4   0.34744    1.997      4     15 0.1467733
variable        1   0.65384   13.222      2     14 0.0005956 ***
group:variable  4   0.51812    1.311      8     30 0.2757596
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here is in addition the SAS output which does not depend on the chosen type (tpes I, II and III yield the same table):

proc mixed data=ldat ;
class id group variable ;
model value = group variable group*variable / s htype=2 ; 
repeated variable / subject=id type=un ; 
run;quit; 

sas table

I am particularly interested in Type II sum of squares. The different ANOVA tables are close from each other but there are some differences (the F value and/or degrees of freedom). I'd like to understand these differences.

For the ez package I don't even know whether the results are given under the sphericity assumption.

$\endgroup$
3
$\begingroup$

I'm afraid SAS gives "wrong" Type II test statistics.

I have learnt a lot of things.

Type II hypothesis

According to these slides by G. Monette and J. Fox the Type II test of an "effect" corresponds to a test for the nullity of a multivariate linear hypothesis about the parameters, say $H_0\colon\{L_{1\mid 2}\beta=0\}$ for a certain matrix $L_{1 \mid 2}$.

Actually $H_0$ is a ``conditional'' hypothesis (see the slides) and $L_{1\mid 2}$ depends on the asymptotic covariance matrix $V$ of $\hat\beta$.

For mixed models and/or models with a repeated structure the $L_{1\mid 2}$ matrix is obtained by estimating $V$, but for classical linear models the $V$ matrix is proportional to ${(X'X)}^{-1}$ and there's no need to estimate it.

Classical linear model example

Consider for instance a classical two-way ANOVA model:

> library(car)
> dat <- Soils[-c(1,2),] # delete two rows of the Soils dataset
> xtabs(~Contour+Depth,data=dat) # the design is unbalanced
            Depth
Contour      0-10 10-30 30-60 60-90
  Depression    4     4     4     4
  Slope         4     4     4     4
  Top           2     4     4     4
> fit <- lm(pH~Contour*Depth, data=dat)

I have written a function which computes the $L_{1 \mid 2}$ matrix which is valid for any model with two crossed fixed factors facAand facB as long as the vcov()function returns the estimated asymptotic covariance matrix :

L12matrix <- function(fit, facA, facB, V=vcov(fit)){
  coefs <- names(coef(fit))
  facA.match <- str_detect(coefs,facA)
  L1 <- NULL
  p <- length(facA.match)
  for(i in 1:p){
    if(facA.match[i]){
      l <- rep(0,p)
      l[i] <- 1
      L1 <- rbind(L1,l)
    }  
  }
  rownames(L1) <- coefs[facA.match]
  facAB.match <- facA.match & str_detect(coefs,facB)
  L2 <- NULL
  for(i in 1:p){
    if(facAB.match[i]){
      l <- rep(0,p)
      l[i] <- 1
      L2 <- rbind(L2,l)
    }  
  }
  rownames(L2) <- coefs[facAB.match]
  L12 <- t(car:::ConjComp(t(L2),t(L1), V))
return(L12)
}

And indeed, with the linearHypothesis() function we get the same results as the Anova() function:

> L12 <- L12matrix(fit,"Contour","Depth")
> linearHypothesis(fit,L12) 
Linear hypothesis test

Hypothesis: ...

Model 1: restricted model
Model 2: pH ~ Contour * Depth

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     36 5.6592                           
2     34 5.3003  2   0.35891 1.1511 0.3283


> L12 <- L12matrix(fit,"Depth","Contour")
> linearHypothesis(fit,L12)
Linear hypothesis test

Hypothesis: .... 

Model 1: restricted model
Model 2: pH ~ Contour * Depth

  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1     37 18.3119                                  
2     34  5.3003  3    13.012 27.822 2.847e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> Anova(fit)
Anova Table (Type II tests)

Response: pH
               Sum Sq Df F value    Pr(>F)    
Contour        0.3589  2  1.1511    0.3283    
Depth         13.0116  3 27.8220 2.847e-09 ***
Contour:Depth  0.5268  6  0.5633    0.7564    
Residuals      5.3003 34                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

A GLS example

Consider the simulated data and the fitted gls model of my previous post:

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 

##
glsfit <- gls(value ~ group*variable, data=ldat,  
            correlation=corSymm(form= ~ 1 | id),
            weights=varIdent(form = ~1 | variable))

The Type II test of the group factor is performed as follows:

> L12 <- L12matrix(glsfit, "group", "variable")
> linearHypothesis(glsfit, L12)
Linear hypothesis test

Hypothesis: ...

Model 1: restricted model
Model 2: value ~ group * variable

  Df  Chisq Pr(>Chisq)  
1                       
2  4 9.5172    0.04939 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Here the linearHypothesis()functin has returned the Wald $\chi^2$-test statistics for $H_0$, and the Wald $F$-test statistics is simply obtained by dividing by the degrees of freedom:

> 9.5172/4
[1] 2.3793

Unfortunately the Anova() function from the car package does not work for gls fitted models. But here the Type II Wald $F$-test statistic is equal to the Type I Wald $F$-test statistic because the design is balanced, and this statistic is returned by the anova() function:

> anova(glsfit)
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

If the design were not balanced then the Type II Wald $F$-test statistic can be obtained by inverting the order of the factors:

gfit.reverse <- update(gfit, model = value ~ variable*group)
anova(gfit.reverse)

Note that the L12matrix(glsfit, ...) command only uses the fitted model glsfit to get the parameters names and to get the asymptotic covariance matrix from the default value of the V argument. Thus we would get the same results by typing :

L12 <- L12matrix(fit, "group", "variable", V=vcov(glsfit))

The SAS bug ?

Now, what's the problem with SAS ? For a classical linear model such as fit, no problem: it gives the same results as the Anova() function.

But for a model with a repeated structure, such as glsfit, the SAS output is obtained as follows with the R tools introduced above:

> # SAS result :
> lmfit <- lm(value ~ group*variable, data=ldat)
> L12 <- L12matrix(glsfit, "group", "variable", V=vcov(lmfit))
> linearHypothesis(glsfit, L12)
Linear hypothesis test

Hypothesis: ...

Model 1: restricted model
Model 2: value ~ group * variable

  Df  Chisq Pr(>Chisq)  
1                       
2  4 7.9864    0.09208 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Divide 7.9864 by 4 and you find the SAS output : enter image description here

In other words, SAS gives the Wald $F$-test statistic for $H_0:\colon\{L_{1\mid 2}\beta=0\}$ but by taking the $L_{1\mid 2}$ matrix obtained with the covariance matrix of the classical linear model obtained by ignoring the repeated structure.

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).

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.