[R] Two-factorial Huynh-Feldt-Test

[Peter Dalgaard]

Bela Bauer <belabauer at cbs.mpg.de> writes:

> Peter Dalgaard wrote:
> 
> >>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 ...
> >
> 
> Yes, of course, the entire SAS script is available to me:
> PROC GLM;
> MODEL col1 col3 col5 col7 col9 col11 col13 col15 col17 col19 col21
> col23 =/nouni;
> repeated roi 6, ord 2/nom mean;
> TITLE 'ABDERUS lat ACC 300-500';
> 
> This script wasn't written by me and unfortunately, I don't know
> anything about SAS scripting, which makes it hard for me to follow the
> script.

OK, so it is basically isomorphic to the structure I had above. You
just have no grouping and no covariate but a 2x6 substructure in the
12 responses.
 
> >http://www.psych.upenn.edu/~baron/rpsych.pdf
> >
> Thanks for that link, it seems like a quite useful paper!
> 
> Now, I've tried to follow it and apply the steps to my own problem,
> but I came up with a H-F and G-G value that is still slightly
> different from the one that SAS calculates (I'll attach my code at the
> end of this mail). With my old code (see my last email), I get
> 0.4682799. With the new code, I get 0.4494588 as G-G epsilon and
> 0.6063626 for the H-F
> correction. I think that this corresponds (or should do so, rather) to
> the SAS value of  0.6333.
> Now, as you can see, the differences are small, but still there.
> The main difference between my two versions of the code are that in
> the old code (based on that book), a covariance matrix is calculcated
> for every condition and then these are averaged for the covariance
> matrix that is used in the H-F formula. The new code, on the other
> hand, averages the data over all conditions and then calculcates the
> covariance matrix for that, which could probably explain the
> differences. There also seems to be a difference between the original
> H-F correction and the algorithm that SAS uses; this is mentioned in
> the PDF, but they don't explain it any further. I suspect that this
> difference could be exactly the difference between my two code
> versions, but I don't know for sure.
> 
> Now, with the difference between my value and the one from SAS being
> so small, I suspect that there's only a very slight difference between
> the algorithms. Do you have any hints what these could be, or how I
> could go about investigating it?
> 
> Thanks a lot for your help!

My suspicion would be that it has something to do with calculating the
correction terms before or after contrast transformations (there must
be a coordinate-free version of the corrections?), but I can't grok
the details that easily. However, if you peek over on the R-devel
list, you will see that I was just about to get serious with
programming some of this stuff as methods for the "mlm" class. I think
you just volunteered to test the code...
 
        -p

> Bela Bauer
> 
> mtx <- NULL
> for (iROI in 1:length(unique( roi )) ) {
>   for (iSubj in 1:length(unique (subj )) ) {
>     mtx <- c(mtx,
>              mean(asa[subj==unique(subj)[iSubj] & roi==unique(roi)[iROI]],
>                   aoa[subj==unique(subj)[iSubj] & roi==unique(roi)[iROI]]))
>   }
> }
> mtx <- matrix(mtx,ncol=length(unique( roi )),byrow=F)
> 
> S <- var(mtx)
> 
> k <- 6
> D <- (k^2 * (mean(S) - mean(diag(S)))^2)
> N1 <- sum(S^2)
> N2 <- 2 * k * sum(apply(S, 1, mean)^2)
> N3 <- k^2 * mean(S)^2
> epsiGG <- D / ((k - 1) * (N1 - N2 + N3))
> epsiHF <- (10 * (k-1) * epsiGG - 2) / ((k-1) * ((10-1) - (k-1)*epsiGG))
> print(epsiGG)
> print(epsiHF)
> 
> ______________________________________________
> 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