[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