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