= 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 [[attachment:sandb.pdf | 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 [[FAQ/td|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.