Effect size for multilevel models
$$\rho = \frac{\mbox{variance(between subjects)}}{\mbox{(variance(between subjects) + variance(within subjects)}}$$
In the wikipedia it states that this "design effect" is used with cluster observations. If we fit a mixed model with a random subject-specific intercept, the clusters are the observations within a participant e. g. the 7 times the participant chose a fruit or a snack. The "design effect" is
D_{eff} = 1 + (m-1) $$\rho$$
where m is the number of observations in each cluster (e.g. number of repeated measures per subject).
The Snijders and Bosker effect size formed from pooling the intercept and residual variances may be evaluated using the below R syntax incorporating a function written by Adam Wagner.
T statistics associated with fixed regression terms may also be used (see here.) which can be converted into Cohen's d and r (correlation) effect sizes for pairwise group and covariate predictors respectively.
#install.packages("stringr")
library(stringr)
library (foreign) # make available data import/export commands
library (graphics) # makes graphics available
library (stats) # makes some important stats available, eg nlm
library (nlme) # makes nlme and lme available
#library (lme4) #not needed as only used with lme which has a VarCorr function of its own; the mod function does not work with VarCorr produced in the lme4 library
BDIData <- read.spss("U:\\My Documents\\BDItest.sav")
BDIData <- data.frame(BDIData)
BDIData$BDI[BDIData$BDI==999] <- NA
BDIData <- na.omit(BDIData)
attach(BDIData)
bdiModelA <- lme(BDI~grp*time,
random=~1|id,
BDIData,
method="ML",
na.action=na.omit)
summary(bdiModelA)
############################################################################################
# Function calculates the proportion of variance explained by each term in an nlme LME model
# Assumes:
# - stringr package loaded
# - nlme package loaded (VarCorr is key for the following code to work)
# Not sure what is coded works in the prescence of interactions of greater order than 2
calcPropExplain <- function(mod){
# Inner function for calculating pooled variance - ignores any slope
pooledVar <- function(model)sum(as.numeric(VarCorr(model)[c("(Intercept)","Residual"), "Variance"]))
# Extract the names of the terms in the model
termLabels <- attr(terms(mod), "term.labels")
# Refit the model length(termLabels) times, where each fit excludes; lapply lists the output from fitting separate models and puts these in modDroppedTerms
# one term (and any interactions associated with that term)
modDroppedTerms <- lapply(termLabels,
function(term){
termsToDrop <- termLabels[str_detect(termLabels, term)]
update(mod, paste(".~.", paste("-", termsToDrop, sep="", collapse="")))})
# Find comparison models for looking at change in proportion variance explained
# Essentially, Where there are interactions involving the term being tested, drop them from the model
# checks for any terms in the model with names a: or :a e.g. grp:time will be flagged as an interaction to drop for grp since 'grp:' is part of grp:time and also for time since
# :time is a part of the string grp:time
compModels <- lapply(termLabels,
function(term){
intsToDrop <- c(termLabels[str_detect(termLabels, paste(term, ":", sep=""))],
termLabels[str_detect(termLabels, paste(":", term, sep=""))])
if(length(intsToDrop)==0) return(mod)
else return(update(mod,
paste(".~.", paste("-", intsToDrop, sep="", collapse=""))))
})
# Find the residual variance for each term; sapply produces a vector of values and puts thses in residVars and compModelVars
residVars <- sapply(modDroppedTerms, pooledVar)
compModelVars <- sapply(compModels, pooledVar)
# Calculate the variance explained by each term
propExplained <- 1-compModelVars/residVars
names(propExplained) <- termLabels
propExplained
}
calcPropExplain(bdiModelA)