`lme`

and `gls`

models`Standardized-mean-differences-in-multi-level-models.Rmd`

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

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`

```
## 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
```

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