lmeInfo provides analytic derivatives and information matrices for fitted linear mixed effects (lme) models and generalized least squares (gls) models estimated using
nlme::gls(), respectively. The package includes functions for estimating the sampling variance-covariance of variance component parameters using the inverse Fisher information. The variance components include the parameters of the random effects structure (for lme models), the variance structure, and the correlation structure. The expected and average forms of the Fisher information matrix are used in the calculations, and models estimated by full maximum likelihood or restricted maximum likelihood are supported. The package also includes a function for estimating standardized mean difference effect sizes (Pustejovsky et al., 2014) based on fitted
You can install the released version of lmeInfo from CRAN with:
You can install the development version from GitHub with:
# install.packages("devtools") devtools::install_github("jepusto/lmeInfo")
We use a dataset from a multiple baseline study conducted by Bryant and colleagues (2016) to demonstrate how
lmeInfo estimates the sampling variance-covariance of variance component parameters of fitted LME models. We also show 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 three schools. After an initial baseline period in each school, 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 school and each student. Here we fit the model using
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(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 #>  338.4864 #> #> attr(,"class") #>  "varcomp"
The sampling variance-covariance of variance component parameters of model
Bryant2016_RML can be estimated using
type = "expected" will calculate the sampling variance-covariance of variance component parameters using the inverse expected Fisher information. Setting
type = "average" will calculate the inverse of the average Fisher information matrix (Gilmour, Thompson, & Cullis, 1995).
varcomp_vcov(Bryant2016_RML, type = "expected") #> Tau.school.var((Intercept)) Tau.case.var((Intercept)) sigma_sq #> Tau.school.var((Intercept)) 5.695541e+04 -4547.41772 0.06132317 #> Tau.case.var((Intercept)) -4.547418e+03 16045.57184 -32.27967543 #> sigma_sq 6.132317e-02 -32.27968 801.20973947
The package also includes a function,
g_mlm(), for estimating a standardized mean difference effect size from a fitted multi-level model. The estimation methods follow Pustejovsky, Hedges, and Shadish (2014). A standardized mean difference effect size parameter can be defined as the ratio of a linear combination of the model’s fixed effect parameters over the square root of a linear combination of the model’s variance components. The
g_mlm() function takes as inputs a fitted multi-level model and the vectors
r_const, which define the linear combinations of fixed effects and variance components, respectively. The function calculates an effect size estimate by first substituting maximum likelihood or restricted maximum likelihood estimates in place of the corresponding parameters, then applying a small-sample correction. The small-sample correction and the standard error are based on approximating the distribution of the estimator by a t distribution, with degrees of freedom given by a Satterthwaite approximation (Pustejovsky, Hedges, & Shadish, 2014). The
g_mlm() function includes an option allowing use of the expected or average form of the Fisher information matrix in the calculations. The
g_mlm() function also includes an option allowing returning the fitted model parameters in addition to effect size estimate.
In our model for the Bryant data, we use the treatment effect in the numerator of the effect size and the sum of the school-level, student-level, and within-student variance components in the denominator of the effect size. The constants are therefore given by
p_const = c(0, 1) and
r_const = c(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
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.786
varcomp_vcov() functions operate on fitted
gls() models. Models fitted using
lme() can include three types of variance component parameters: random effects variances and covariances, correlation structure parameters, and variance structure parameters. Models fitted using
gls() can include correlation structure parameters and variance structure parameters. The
nlme package provides many different forms for each of these components, not all of which are supported in
lmeInfo. The package can handle the following classes of variance components:
pdSymmmatrices, including in the
corARMAfor MA(1) models only
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. https://doi.org/10.1177/0022219414538516
Gilmour, A. R., Thompson, R., & Cullis, B. R. (1995). Average information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics, 51(4), 1440–1450. https://doi.org/10.2307/2533274
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. https://doi.org/10.3102/1076998614547577
Wang T, Merkle EC (2018). merDeriv: Derivative Computations for Linear Mixed Effects Models with Application to Robust Standard Errors. Journal of Statistical Software, Code Snippets, 87(1), 1-16. https://doi.org/10.18637/jss.v087.c01.