[R] library(car): Anova and repeated measures without between subjects factors

[John Fox]

Dear Ralf,

I've now had a chance to take a closer look at the problems that you
reported with Anova.mlm(), and uploaded a new version of the car package to
CRAN that deals with them. Please see below:

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of Ralf Goertz
> Sent: Wednesday, October 17, 2007 6:00 AM
> To: John Fox
> Cc: r-help at r-project.org
> Subject: Re: [R]library(car): Anova and repeated measures 
> without between subjects factors
> 
> John Fox, Dienstag, 16. Oktober 2007:
> > Dear Ralf,
> > 
> > Unfortunately, Anova.mlm(), and indeed Anova() more 
> generally, won't 
> > handle a model with only a constant. As you point out, this isn't 
> > reasonable for repeated-measures ANOVA, where it should be 
> possible to 
> > have only within-subjects factors. When I have a chance, 
> I'll see what 
> > I can do to fix the problem -- my guess is that it shouldn't be too 
> > hard.
> >
> > Thanks for pointing out this limitation in Anova.mlm()

Anova.mlm() was capable of dealing with models with only an intercept in the
between-subject part of the model if the argument type="III" was specified.
Because there is no distinction between "type II" and "type III" tests for
models with only an intercept [Anova.mlm() insures that the coding of
different terms in the within-subjects part of the model is orthogonal],
I've simply let the former invoke the latter. 

Actually, there's no really good reason that Anova() should include a line
for the intercept when type="III" but not when type="II"; it was just easier
to do things this way. Eventually, I'll include a test for the intercept
when type="II".

Here's how your example works out currently [comparing Anova() with
anova()]:

> anova(lm(cbind(mat.1, mat.2, mat.3) ~ 1, data=data), 
+     idata=data.frame(zeit=factor(1:3)), X=~1)
Analysis of Variance Table


Contrasts orthogonal to
~1

            Df Pillai approx F num Df den Df   Pr(>F)   
(Intercept)  1 0.3488   7.4976      2     28 0.002468 **
Residuals   29                                          
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> Anova(lm(cbind(mat.1, mat.2, mat.3) ~ 1, data=data), 
+     idata=data.frame(zeit=factor(1:3)),
+     idesign=~zeit)

Type III Repeated Measures MANOVA Tests: Pillai test statistic
            Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)  1      0.99  2077.34      1     29 < 2.2e-16 ***
zeit         1      0.35     7.50      2     28  0.002468 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
Warning message:
In Anova.mlm(lm(cbind(mat.1, mat.2, mat.3) ~ 1, data = data), idata =
data.frame(zeit = factor(1:3)),  :
  the model contains only an intercept: equivalent Type III test substituted

(You can, of course, get the test assuming sphericity as well.)

> 
> Dear John,
> 
> I am looking forward to your having a chance. There is one 
> thing that I would like to request, though. 
> Greenhouse-Geisser and Huyn-Feldt eps corrections have 
> already been implemented but how about Mauchly's sphericity 
> test? I know this can be done with mauchly.test() but it 
> would be nice to have it in the summary of Anova().

I've added some code to include sphericity tests in the output from
summary.Anova.mlm(), but haven't yet had a chance to test it adequately. If
it checks out, I'll include it in the next version of the car package.

> 
> However, there is one more thing. Look at the following data
> 
> > c1<-c(-6.0,-10.3,-2.9,-8.3,-10.0,5.3,-7.7,-0.8,9.1,-6.2)
> > mat<-matrix(c(c1,c1),10,2)
> > mat
>        [,1]  [,2]
>  [1,]  -6.0  -6.0
>  [2,] -10.3 -10.3
>  [3,]  -2.9  -2.9
>  [4,]  -8.3  -8.3
>  [5,] -10.0 -10.0
>  [6,]   5.3   5.3
>  [7,]  -7.7  -7.7
>  [8,]  -0.8  -0.8
>  [9,]   9.1   9.1
> [10,]  -6.2  -6.2
> 
> > bf<-ordered(rep(1:2,5))
> > bf
>  [1] 1 2 1 2 1 2 1 2 1 2
> Levels: 1 < 2
> 
> Since the two columns of mat are equal:
> 
> > t.test(mat[,1],mat[,2],paired=T)
> 
>         Paired t-test
> 
> data:  mat[, 1] and mat[, 2]
> t = NaN, df = 9, p-value = NA
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
>  NaN NaN
> sample estimates:
> mean of the differences 
>                       0 
> 
> I would assume to either get a warning or a F-value of 0 for 
> the repeated factor zeit but actually:
> 
> > Anova(lm(mat~bf),idata=data.frame(zeit=ordered(1:2)),idesign=~zeit)
> 
> Type II Repeated Measures MANOVA Tests: Pillai test statistic
>         Df test stat approx F num Df den Df Pr(>F)
> bf       1    0.0020   0.0163      1      8 0.9016
> zeit     1    0.2924   3.3059      1      8 0.1065
> bf:zeit  1    0.0028   0.0221      1      8 0.8854
> 
> whereas
> 
> > anova.mlm(lm(mat~bf),X=~1,idata=data.frame(zeit=ordered(1:2)))
> 
> Error in anova.mlm(...) :
>  residuals have rank 1 < 2
> 
> This is quite dangerous. In a real data situation I 
> accidentally used the same column twice and I got a 
> significant effect for the factor zeit! I hope it wouldn't be 
> too hard to fix this. too.
> 

linear.hypothesis.mlm() now tries to find the rank of the error SSP matrix
and throws an error if it is deficient. For your example:

> Anova(lm(mat ~ bf), idata=data.frame(zeit=ordered(1:2)), idesign=~zeit)
Error in linear.hypothesis.mlm(mod, hyp.matrix.1, SSPE = SSPE, V = V,  : 
  The error SSP matrix is apparently of deficient rank = 0 < 1

(I found it interesting that, in the previous version, defining bf as an
ordered factor resulted in the erroneous output that you reported, since
using orthogonal-polynomial contrasts produced small rounding errors, but
defining bf as a factor did not, since the computations using "sum"
contrasts proved to be exact. Of course, relying on an exact rank deficiency
to reveal an error is a very bad idea!)

Thanks again for bringing these issues to my attention.

John