[R] Two-factorial Huynh-Feldt-Test

[Peter Dalgaard]

Bela Bauer <bela_b at gmx.net> writes:

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

[There's no "Huynh-Feldt-Test", there's a H-F epsilon, which is a
correction to the F test]
 
> 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)

Would you happen to know what logic SAS uses to recognize cases  where
the corrections apply?  I thought it only did it in the case of
repeated measurements, as in

proc glm;
        model bmc2-bmc7=bmc1 grp / nouni;
        repeated time ...


 
> 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...:-(

I went looking recently and all *I* got was SPSS examples... As it
turned out, Jon Baron/Yuelin Li has the H-F and G-G calculations
neatly outlined with R code and everything in

http://www.psych.upenn.edu/~baron/rpsych.pdf

> 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 )))
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907