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)which gives R-squareds as output for each model effect (usually > 0)
grp time grp:time 0.0031146434 0.1078555113 0.0005971626
You can also compute Cohen's d effect sizes directly using the eff_size function which works in conjunction with the emmeans which one can output from a regression, in this case from the mixed model routine, lmer. An example of its use is below. Thanks to Andrea Kusec for this.
modelbadsgroup2 <-lmer(badst ~ time + group + time*group + (1 | id), data = longdatat1t2) m.bads2 <-emmeans(modelbadsgroup2, "group","time") m.bads2 > m.bads2 time = 1: group emmean SE df lower.CL upper.CL 1 93.9 4.41 104 85.1 102.6 2 84.2 4.48 105 75.3 93.1 3 90.5 4.94 127 80.8 100.3 time = 2: group emmean SE df lower.CL upper.CL 1 104.2 4.65 113 95.0 113.4 2 90.9 4.74 114 81.5 100.3 3 93.6 5.15 126 83.5 103.8 > eff_size(m.bads2, sigma=sigma(modelbadsgroup2), edf = df.residual(modelbadsgroup2)) time = 1: contrast effect.size SE df lower.CL upper.CL 1 - 2 0.739 0.441 104 -0.1350 1.612 1 - 3 0.255 0.400 104 -0.5387 1.049 2 - 3 -0.484 0.415 105 -1.3068 0.340 time = 2: contrast effect.size SE df lower.CL upper.CL 1 - 2 1.018 0.474 113 0.0783 1.957 1 - 3 0.806 0.429 113 -0.0434 1.655 2 - 3 -0.212 0.442 114 -1.0885 0.664
Reference
Snijders TAB and Bosker RJ (1999) Multilevel Analysis: An Introduction to Basic and Advanced Multilevel Modeling. Sage:London.