stefan.premke at email.de
2006-May-15 13:40 UTC
[R] extracting the F-Values from a aov-call to calculate a green-gei-corrected p-values
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