## Standardized mean difference effect sizes in linear mixed effects models

The 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).

## Demonstration

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               12

The estimated variance components from the fitted model can be obtained using extract_varcomp():

extract_varcomp(Bryant2016_RML)
## $Tau ##$Tau$case ## case.var((Intercept)) ## 255.4297 ## ##$Tau$school ## school.var((Intercept)) ## 158.24 ## ## ##$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 <- g_mlm(Bryant2016_RML, p_const = c(0,1), r_const = c(1,1,1), infotype = "expected")
Bryant2016_g
##                           est    se
## unadjusted effect size  1.799 0.340
## adjusted effect size    1.721 0.325
## degree of freedom      17.504

A summary() method is also included, which includes more detail about the model parameter estimates and effect size estimate:

summary(Bryant2016_g)
##                                            est      se
## Tau.case.case.var((Intercept))         255.430 238.653
## Tau.school.school.var((Intercept))     158.240 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.786

## References

Bryant, 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