[R] extracting the F-Values from a aov-call to calculate a green-gei-corrected p-values

stefan.premke at email.de stefan.premke at email.de
Mon May 15 15:40:25 CEST 2006


Hello all!
considering the following data

group=factor(rep(rep(c(1,2),c(7,7)),8))
subj=factor(rep(seq(1,14,1),8))
cond=factor(rep(rep(c(1,2),c(14,14)),4))
roi=factor(rep(c(1:4),rep(28,4)))
signal=round(rnorm(112,1,2),2)
data=data.frame(group,subj,cond,roi,signal)

I want to calculate ANOVA with group between subject factor and cond, roi within subject factors so I type

t=aov(signal~group*cond*roi+Error(subj/(cond*roi)),data=data)

now I want to result of the significance tests so I type

summary(t) 
which results in a output of F and p-values (I did not copy it here, because I gave you code to produce the output on your own computer)

now I come up with the Idea to calculate a Greenhouse Geisser Correction (or Box adjustment or both)
as the ?anova.mlm results in a relativly cryptic and not too helpful help for me, I decide to simply calculate the corrected p-values by myself. 
I can take the code given in the R-Archivs written by John Christie  shortly requoted here and adapt it to my dataset
# This returns the Huynh-Feldt or "Box Correction" for degrees of
freedom

m=data[,-1]
hf <- function(m){
        # m is a matrix with subjects as rows and conditions as columns
        # note that checking for worst case scenarios F correction first might
        # be a good idea using J/(J-1) as the df correction factor
        n<- length(m[,1])
        J<-length(m[1,])
        X<-cov(m)*(n-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<- (J^2*(D-SPm)^2) / ((J-1) * (SPm2 - 2*J*SSrm + J^2*SPm^2))
        epsilon
} 
box=hf(m)
box
[1] 0.5434483

So this gives me the epsilon. All I have to do is to multiply this with the num_df and den_df and type

1-pf(F-value,new_num_df,new_den_df)

which leaves me with what I thought to be very simple in the first place. 
I just need to extract the F-values out of the aov-call, ie from the list it returns 
First I expected (intuitively) that there will be a list-element named "F-Values" or so, which just stores the F-values. but I did not see that. 
Ok I can calculate the F-Value by meanSq(factor_of_interest) / meanSqres(from the proper Error stratum)
(my Factors of interest are group, cond, group:cond, roi, group:roi, cond:roi, group:cond:roi)
But I also do not find the meanSq. 
So I can calculate them by sumSq / df
I can find the the sumSqres (sum(t[[i]]$res^2)) and the res_df (t[[i]]$df) (i in 2:5) but this only calculates the denominator of the "F-Value".
I could not find out how to get the numerator and his df (other then with eye, paper and pencil)

by try and error I found (but do not understand) that
sum(t[[i]]$eff^2) gives the total SS (which would solve the problem for t[[2]] if it is correct and not incidence) but not for the other
sum(x[[3]]$fit^2) which gives something like total SS - Res SS (not sure if it makes sense). but again this only helps (if it helps anyway) for t[[2]]

so to finish with a summary of my question
how to extract the F-values and the dfs out of aov-call other then by eye, paper and pencil, ie automatic

sincerely
stefan




More information about the R-help mailing list