Rutils/maybe-not-useful/F.test.r

In manfredo89/ED2io: Manages Input/Output for Ecosystem Demography 2 Model

#==========================================================================================#
#==========================================================================================#
#     F.test.                                                                              #
#     This function performs a F-test with the samples, given as multiple elements in a    #
# list.                                                                                    #
#------------------------------------------------------------------------------------------#
F.test <<- function(x,...){
   #----- Check that x is a list. ---------------------------------------------------------#
   if (! is.list(x)) stop("x is not a list!")
   #---------------------------------------------------------------------------------------#


   #----- Eliminate NA from lists. --------------------------------------------------------#
   x  = lapply(X=x,FUN=na.exclude)
   ux = unlist(x)
   fx = as.factor(unlist(mapply(FUN=rep,x=seq_along(x),each=sapply(X=x,FUN=length))))
   #---------------------------------------------------------------------------------------#




   #---------------------------------------------------------------------------------------#
   #     Perform a Levene's test to check for homoscedasticity.                            #
   #---------------------------------------------------------------------------------------#
   l.test = leveneTest(y=ux,group=fx)
   l.test = list( statistic = c(l.test[["F value"]])[[1]][1]
                , p.value   = c(l.test[["Pr(>F)" ]])[[1]][1]
                )#end list
   #---------------------------------------------------------------------------------------#



   #---------------------------------------------------------------------------------------#
   #      Initialise the ANOVA table.                                                      #
   #---------------------------------------------------------------------------------------#
   aovtab           = data.frame(sumsq=rep(NA,3),df=rep(NA,3),msq=rep(NA,3))
   rownames(aovtab) = c("Between groups","Within groups","Total")
   #---------------------------------------------------------------------------------------#



   #---------------------------------------------------------------------------------------#
   #     Find the SoSq and degrees of freedom of the full model.                           #
   #---------------------------------------------------------------------------------------#
   aovtab$sumsq[2] = sum((unlist(lapply(X=x,FUN=function(x) x - mean(x))))^2)
   aovtab$df   [2] = sum(sapply(X=lapply(X=x,FUN=length),FUN="-",1))
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #     Find the SoSq and degrees of freedom of the reduced model.                        #
   #---------------------------------------------------------------------------------------#
   aovtab$sumsq[3] = sum((ux-mean(ux))^2)
   aovtab$df   [3] = length(ux)-1
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #      Find the SoSq and degrees of freedom between groups.                             #
   #---------------------------------------------------------------------------------------#
   aovtab$sumsq[1] = aovtab$sumsq[3] - aovtab$sumsq[2]
   aovtab$df   [1] = aovtab$df   [3] - aovtab$df   [2]
   #---------------------------------------------------------------------------------------#



   #---------------------------------------------------------------------------------------#
   #      Find the mean square.                                                            #
   #---------------------------------------------------------------------------------------#
   aovtab$msq = aovtab$sumsq / aovtab$df
   #---------------------------------------------------------------------------------------#

   #---------------------------------------------------------------------------------------#
   #    The F statistic is the ratio between msq between groups and within groups.         #
   #---------------------------------------------------------------------------------------#
   sizes   = sapply(X=x,FUN=length)
   means   = sapply(X=x,FUN=mean)
   vars    = sapply(X=x,FUN=var)
   F.value = aovtab$msq[1] / aovtab$msq[2]
   p.value = 1. - pf(q=F.value,df1=aovtab$df[1],df2=aovtab$df[2])
   #---------------------------------------------------------------------------------------#



   return( list( statistic = F.value
               , mean      = means
               , var       = vars
               , sizes     = sizes
               , p.value   = p.value
               , aovtab    = aovtab
               , levene    = l.test
               )#end list
         )#end return
}#end F.test
#==========================================================================================#
#==========================================================================================#