<?xml version="1.0" encoding="utf-8"?><!DOCTYPE article  PUBLIC '-//OASIS//DTD DocBook XML V4.4//EN'  'http://www.docbook.org/xml/4.4/docbookx.dtd'><article><articleinfo><title>FAQ/reffs</title><revhistory><revision><revnumber>17</revnumber><date>2021-09-23 08:20:16</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>16</revnumber><date>2021-09-23 08:18:48</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>15</revnumber><date>2021-09-23 08:10:19</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>14</revnumber><date>2016-02-26 11:07:52</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>13</revnumber><date>2014-09-04 14:37:37</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>12</revnumber><date>2014-09-04 08:51:07</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>11</revnumber><date>2014-09-04 08:50:51</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>10</revnumber><date>2014-09-04 08:33:01</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>9</revnumber><date>2014-09-04 08:30:52</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>8</revnumber><date>2014-09-04 08:30:34</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>7</revnumber><date>2014-09-04 08:27:59</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>6</revnumber><date>2013-03-08 10:17:13</date><authorinitials>localhost</authorinitials><revremark>converted to 1.6 markup</revremark></revision><revision><revnumber>5</revnumber><date>2012-12-06 12:57:47</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>4</revnumber><date>2012-12-06 12:57:27</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>3</revnumber><date>2012-12-06 12:57:00</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>2</revnumber><date>2012-12-06 12:56:27</date><authorinitials>PeterWatson</authorinitials></revision><revision><revnumber>1</revnumber><date>2012-12-06 12:55:42</date><authorinitials>PeterWatson</authorinitials></revision></revhistory></articleinfo><section><title>Effect size for multilevel models</title><para>$$\rho = \frac{\mbox{variance(between subjects)}}{\mbox{(variance(between subjects) + variance(within subjects))}}$$ </para><para>In the wikipedia it states that this &quot;design effect&quot; is used with cluster observations. If we fit a mixed model with a <emphasis role="strong">random subject-specific intercept</emphasis>, the clusters are the observations within a participant e. g. the 7 times the participant chose a fruit or a snack. The &quot;design effect&quot; is </para><para>D_{eff} = 1 + (m-1) $$\rho$$ </para><para>where m is the number of observations in each cluster (e.g. number of repeated measures per subject). </para><para>The <ulink url="https://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/reffs/statswiki/FAQ/reffs?action=AttachFile&amp;do=get&amp;target=sandb.pdf">Snijders and Bosker</ulink> 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. </para><para>T statistics associated with fixed regression terms may also be used (see <ulink url="https://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/reffs/statswiki/FAQ/td#">here.</ulink>) which can be converted into Cohen's d and r (correlation) effect sizes for pairwise group and covariate predictors respectively. </para><screen><![CDATA[#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
]]><![CDATA[
BDIData <- read.spss("U:\\My Documents\\BDItest.sav")
BDIData <- data.frame(BDIData)
BDIData$BDI[BDIData$BDI==999] <- NA
BDIData <- na.omit(BDIData) 
attach(BDIData)
]]><![CDATA[
bdiModelA   <- lme(BDI~grp*time,
                   random=~1|id,
                   BDIData,
                   method="ML",
                   na.action=na.omit)
]]><![CDATA[
summary(bdiModelA)
]]><![CDATA[
############################################################################################
# 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)]]></screen><para>which gives R-squareds as output for each model effect (usually &gt; 0) </para><screen><![CDATA[         grp         time      grp:time 
0.0031146434 0.1078555113 0.0005971626]]></screen><para>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. </para><screen><![CDATA[modelbadsgroup2 <-lmer(badst ~ time + group + time*group + (1 | id), data = longdatat1t2)
m.bads2 <-emmeans(modelbadsgroup2, "group","time")
m.bads2
]]><![CDATA[
> 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
]]><![CDATA[
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
]]><![CDATA[
> 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
]]><![CDATA[
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]]></screen><para><emphasis role="underline">Reference</emphasis> </para><para>Snijders TAB and Bosker RJ (1999) Multilevel Analysis: An Introduction to Basic and Advanced Multilevel Modeling. Sage:London. </para></section></article>