lme and gls modelsStandardized-mean-differences-in-multi-level-models.RmdThe lmeInfo package includes a function,
g_mlm(), for estimating a standardized mean difference
effect size from a multi-level model fitted using lme() or
gls() from the nlme package. The estimation
methods follow Pustejovsky, Hedges, and Shadish (2014). Let \(\boldsymbol{\beta}\) denote the fixed
effect parameters and \(\boldsymbol{\theta}\) denote the variance
components (in the natural, variance-covariance parameterization) from a
multi-level model. A standardized mean difference effect size can then
be defined as \[
\delta =
\frac{\mathbf{p}'\boldsymbol\beta}{\sqrt{\mathbf{r}'\boldsymbol\theta}},
\] for some fixed vectors \(\mathbf{p}\) and \(\mathbf{r}\).
The g_mlm() function takes as inputs a fitted
lme() model and the vectors \(\mathbf{p}\) and \(\mathbf{r}\). It computes an estimate of
\(\delta\) as \[
d = \left(1 - \frac{3}{4 \nu - 1}\right)
\frac{\mathbf{p}'\boldsymbol{\hat\beta}}{\sqrt{\mathbf{r}'\boldsymbol{\hat\theta}}},
\] where \(\boldsymbol{\hat\beta}\) and \(\boldsymbol{\hat\theta}\) are maximum
likelihood or restricted maximum likelihood estimates of the fixed
effects and variance components, respectively. Here, \(\nu\) is the Satterthwaite degrees of
freedom given by \[
\nu = \frac{2
\left(\mathbf{r}'\boldsymbol{\hat\theta}\right)^2}{\text{Var}\left(\mathbf{r}'\boldsymbol{\hat\theta}\right)},
\] and the variance of \(\mathbf{r}'\boldsymbol{\hat\theta}\) is
calculated using the inverse of the Fisher information matrix. An
approximate standard error for the standardized mean difference estimate
is \[
\text{SE}(d) \approx \left(1 - \frac{3}{4 \nu - 1}\right)
\sqrt{\frac{\nu}{\nu - 2} \times
\frac{\text{Var}\left(\mathbf{p}'\boldsymbol{\hat\beta}\right)}{\mathbf{r}'\boldsymbol{\hat\theta}}
+ d^2 \times \frac{8 \nu^2 - \nu + 2}{16 (\nu - 2) (\nu - 1)^2}}
\] (Pustejovsky, Hedges, & Shadish, 2014).
We use a dataset from a multiple baseline study conducted by Bryant and colleagues (2016) to demonstrate how to calculate a design-comparable standardized mean difference effect size based on the fitted model.
The study by Bryant and colleagues (2016) involved collecting
repeated measures of math performance on multiple students, in each of
several classrooms. After an initial baseline period in each classroom,
an intervention was introduced and its effects on student math
performance were observed over time. For sake of illustration, we use a
very simple model for these data, consisting of a simple change in
levels coinciding with the introduction of treatment. We include random
effects for each classroom and each student. Here we fit the model using
nlme::lme():
library(lmeInfo)
library(nlme)
data(Bryant2016)
Bryant2016_RML <- lme(fixed = outcome ~ treatment,
random = ~ 1 | school/case,
data = Bryant2016)
summary(Bryant2016_RML)## Linear mixed-effects model fit by REML
## Data: Bryant2016
## AIC BIC logLik
## 2627.572 2646.041 -1308.786
##
## Random effects:
## Formula: ~1 | school
## (Intercept)
## StdDev: 12.57935
##
## Formula: ~1 | case %in% school
## (Intercept) Residual
## StdDev: 15.98217 18.398
##
## Fixed effects: outcome ~ treatment
## Value Std.Error DF t-value p-value
## (Intercept) 56.14333 9.019268 286 6.224821 0
## treatmenttreatment 49.33454 2.399065 286 20.564070 0
## Correlation:
## (Intr)
## treatmenttreatment -0.194
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.15430670 -0.65412204 0.02112957 0.61701090 2.90624593
##
## Number of Observations: 299
## Number of Groups:
## school case %in% school
## 3 12The estimated variance components from the fitted model can be
obtained using extract_varcomp():
extract_varcomp(Bryant2016_RML)## $Tau
## $Tau$school
## school.var((Intercept))
## 158.24
##
## $Tau$case
## case.var((Intercept))
## 255.4297
##
##
## $cor_params
## numeric(0)
##
## $var_params
## numeric(0)
##
## $sigma_sq
## [1] 338.4864
##
## attr(,"class")
## [1] "varcomp"In our model for the Bryant data, we use the treatment effect in the numerator of the effect size and the sum of the classroom-level, student-level, and within-student variance components in the denominator of the effect size. The constants are therefore given by \(\mathbf{p} = (0, 1)'\) and \(\mathbf{r} = (1, 1, 1)'\). The effect size estimate can be calculated as:
Bryant2016_g## est se
## unadjusted effect size 1.799 0.340
## adjusted effect size 1.721 0.325
## degree of freedom 17.504A summary() method is also included, which includes more
detail about the model parameter estimates and effect size estimate:
summary(Bryant2016_g)## est se
## Tau.school.school.var((Intercept)) 158.240 238.653
## Tau.case.case.var((Intercept)) 255.430 126.671
## sigma_sq 338.486 28.306
## total variance 752.156 254.250
## (Intercept) 56.143 9.019
## treatmenttreatment 49.335 2.399
## treatment effect at a specified time 49.335 2.399
## unadjusted effect size 1.799 0.340
## adjusted effect size 1.721 0.325
## degree of freedom 17.504
## constant kappa 0.087
## logLik -1308.786Bryant, B. R., Bryant, D. P., Porterfield, J., Dennis, M. S., Falcomata, T., Valentine, C., … & Bell, K. (2016). The effects of a Tier 3 intervention on the mathematics performance of second grade students with severe mathematics difficulties. Journal of Learning Disabilities, 49(2), 176-188. doi:10.1177/0022219414538516
Pustejovsky, J. E., Hedges, L. V., & Shadish, W. R. (2014). Design-comparable effect sizes in multiple baseline designs: A general modeling framework. Journal of Educational and Behavioral Statistics, 39(5), 368–393. doi:10.3102/1076998614547577