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