[R] Two-factorial Huynh-Feldt-Test

[Bela Bauer]

Hi,

I'm currently working on porting some SAS scripts to R, and hence need 
to do the same calculation (and get the same results) as SAS in order to 
make the transition easier for users of the script.
In the script, I'm dealing with a two-factorial repeated-measures anova. 
I'll try to give you a short overview of the setup:

- two between-cell factors: facBetweenROI (numbering regions of 
interest, 6 levels) and facBetweenCond (numbering conditions, 2 levels).
- one within-cell factor: facWithin (which numbers subjects, while 
subjects are considered repetitions of the same measurement). This is 
basically the repeatedness of the data.

For this data, I do anovas for several linear models. SAS also 
calculates the Huynh-Feldt-Test, which is in this case very important to 
the users and cannot be replaced with nlme or something of the kind (as 
recommended in http://maths.newcastle.edu.au/~rking/R/help/03b/0813.html.

The models I use for the anovas are the following:

aov(vecData ~ (facWithin + facBetweenROI + facBetweenCond)^2)
aov(vecData ~ facBetweenROI + facBetweenCond %in% facWithin + 
Error(facBetweenROI %in% facWithin))
aov(vecData ~ facBetweenCond + facBetweenROI %in% facWithin + 
Error(facBetweenCond %in% facWithin))

SAS seems to calculate the Huynh-Feldt test for the first and the second 
model. The SAS output is (for the second aov)

Source                    DF   Type III SS   Mean Square  F Value  Pr > F

 roi                        5   258.7726806    51.7545361    21.28  <.0001
 Error(roi)                95   231.0511739     2.4321176

                                           Adj Pr > F
                  Source                 G - G     H - F

                  roi                   <.0001    <.0001
                  Error(roi)


                   Greenhouse-Geisser Epsilon    0.5367
                   Huynh-Feldt Epsilon           0.6333

and for the first one:

Source                    DF   Type III SS   Mean Square  F Value  Pr > F

 ord                        1     2.2104107     2.2104107     0.24  0.6276
 Error(ord)                19   172.7047994     9.0897263


 Source                    DF   Type III SS   Mean Square  F Value  Pr > F

 roi*ord                    5   10.37034918    2.07406984     2.89  0.0180
 Error(roi*ord)            95   68.25732255    0.71849813

                                           Adj Pr > F
                  Source                 G - G     H - F

                  roi*ord               0.0663    0.0591
                  Error(roi*ord)

Now, to do this in R, I have a short (and not very thorough) description 
of it in "Lehrbuch der Statistik", Bortz and I have the script from 
http://maths.newcastle.edu.au/~rking/R/help/03b/7891.html
Still, I can't figure out how to do this correctly in the two-factorial 
way (and with the different models that SAS seems to use.) I'll attach 
my code at the end, in case there's something that makes sense about it. 
I've tried to do it in several ways, but this is the only one that gives 
a somewhat reasonable result.

Can anyone give me a suggestion of how I could do this, where I could 
find information about it etc? Google doesn't help much except more SAS 
examples...:-(

hf <- function(mtxCov,ncol,nrow) {
  X <- mtxCov*(nrow-1)
  r <- length(X[,1])
  D <- 0
  for (i in 1: r) D<- D+ X[i,i]
  D <- D/r
  SPm  <- mean(X)
  SPm2 <- sum(X^2)
  SSrm <- 0
  for (i in 1: r) SSrm<- SSrm + mean(X[i,])^2
  epsilon <- (ncol^2*(D-SPm)^2) / ((ncol-1) * (SPm2 - 2*ncol*SSrm + 
ncol^2*SPm^2))
  print(epsilon)
}

# 2. do variance-covariance matrices for conditions first
avCov2 <- matrix(rep(0,36),ncol=length(unique( roi )))
for (iCond in 1:length(preparedData[1,3:4])) {
  mtx <- NULL
  for (iROI in 1:length(unique( roi ))) {
    mtx <- c(mtx,
             preparedData[roi==unique(roi)[iROI],iCond + 4])
  }
  mtx <- matrix(mtx, ncol=length(unique( roi )), byrow=F)
  avCov2 <- avCov2 + cov(mtx)
}
avCov2 <- avCov2 / length(preparedData[1,3:4])

hf(avCov2,
   length(unique( roi )),
   length(unique( subj )))