[R] Two-factorial Huynh-Feldt-Test

Bela Bauer
2005-02-18T08:05:15Z

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