This vignette is intended to provide an overview of the analysis of
multiply imputed data sets with `mitml`

. Specifically, this
vignette addresses the following topics:

- Working with multiply imputed data sets
- Rubin’s rules for pooling individual parameters
- Model comparisons
- Parameter constraints

Further information can be found in the other vignettes and the package documentation.

`studentratings`

)For the purposes of this vignette, we make use of the
`studentratings`

data set, which contains simulated data from
750 students in 50 schools including scores on reading and math
achievement, socioeconomic status (SES), and ratings on school and
classroom environment.

The package and the data set can be loaded as follows.

```
library(mitml)
library(lme4)
data(studentratings)
```

As evident from its `summary`

, most variables in the data
set contain missing values.

`summary(studentratings)`

```
# ID FedState Sex MathAchiev MathDis
# Min. :1001 B :375 Length:750 Min. :225.0 Min. :0.2987
# 1st Qu.:1013 SH:375 Class :character 1st Qu.:440.7 1st Qu.:1.9594
# Median :1513 Mode :character Median :492.7 Median :2.4350
# Mean :1513 Mean :495.4 Mean :2.4717
# 3rd Qu.:2013 3rd Qu.:553.2 3rd Qu.:3.0113
# Max. :2025 Max. :808.1 Max. :4.7888
# NA's :132 NA's :466
# SES ReadAchiev ReadDis CognAbility SchClimate
# Min. :-9.00 Min. :191.1 Min. :0.7637 Min. :28.89 Min. :0.02449
# 1st Qu.:35.00 1st Qu.:427.4 1st Qu.:2.1249 1st Qu.:43.80 1st Qu.:1.15338
# Median :46.00 Median :490.2 Median :2.5300 Median :48.69 Median :1.65636
# Mean :46.55 Mean :489.9 Mean :2.5899 Mean :48.82 Mean :1.73196
# 3rd Qu.:59.00 3rd Qu.:558.4 3rd Qu.:3.0663 3rd Qu.:53.94 3rd Qu.:2.24018
# Max. :93.00 Max. :818.5 Max. :4.8554 Max. :71.29 Max. :4.19316
# NA's :281 NA's :153 NA's :140
```

In the present example, we investigate the differences in mathematics achievement that can be attributed to differences in SES when controlling for students’ sex. Specifically, we are interested in the following model.

\[ \mathit{MA}_{ij} = \gamma_{00} + \gamma_{10} \mathit{Sex}_{ij} + \gamma_{20} (\mathit{SES}_{ij}-\overline{\mathit{SES}}_{\bullet j}) + \gamma_{01} \overline{\mathit{SES}}_{\bullet j} + u_{0j} + e_{ij} \]

Note that this model also employs group-mean centering to separate the individual and group-level effects of SES.

In the present example, we generate 20 imputations from the following imputation model.

```
<- ReadDis + SES ~ 1 + Sex + (1|ID)
fml <- panImpute(studentratings, formula = fml, n.burn = 5000, n.iter = 200, m = 20) imp
```

The completed data are then extracted with
`mitmlComplete`

.

`<- mitmlComplete(imp, "all") implist `

In empirical research, the raw data rarely enter the analyses but
often require to be transformed beforehand. For this purpose, the
`mitml`

package provides the `within`

function,
which applies a given transformation directly to each data set.

In the following, we use this to (a) calculate the group means of SES and (b) center the individual scores around their group means.

```
<- within(implist, {
implist <- clusterMeans(SES, ID) # calculate group means
G.SES <- SES - G.SES # center around group means
I.SES })
```

This method can be used to apply arbitrary transformations to all of the completed data sets simultaneously.

Note regarding`dplyr`

:Due to how it is implemented,`within`

cannot be used directly with`dplyr`

. Instead, users may use`with`

instead of`within`

with the following workaround.`<- with(implist,{ implist <- data.frame(as.list(environment())) df <- ... # dplyr commands df df })<- as.mitml.list(implist) implist`

Advanced users may also consider using

`lapply`

for a similar workaround.`

In order to analyze the imputed data, each data set is analyzed using
regular complete-data techniques. For this purpose, `mitml`

offers the `with`

function. In the present example, we use it
to fit the model of interest with the R package `lme4`

.

```
<- with(implist, {
fit lmer(MathAchiev ~ 1 + Sex + I.SES + G.SES + (1|ID))
})
```

This results in a list of fitted models, one for each of the imputed data sets.

The results obtained from the imputed data sets must be pooled in order to obtain a set of final parameter estimates and inferences. In the following, we employ a number of different pooling methods that can be used to address common statistical tasks, for example, for (a) estimating and testing individual parameters, (b) model comparisons, and (c) tests of constraints about one or several parameters.

Individual parameters are commonly pooled with the rules developed by
Rubin (1987). In `mitml`

, Rubin’s rules are implemented in
the `testEstimates`

function.

`testEstimates(fit)`

```
#
# Call:
#
# testEstimates(model = fit)
#
# Final parameter estimates and inferences obtained from 20 imputed data sets.
#
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 433.015 28.481 15.203 1.081e+03 0.000 0.153 0.134
# SexGirl 3.380 7.335 0.461 2.794e+05 0.645 0.008 0.008
# I.SES 0.692 0.257 2.690 2.334e+02 0.008 0.399 0.291
# G.SES 1.296 0.597 2.173 1.097e+03 0.030 0.152 0.133
#
# Unadjusted hypothesis test as appropriate in larger samples.
```

In addition, the argument `extra.pars = TRUE`

can be used
to obtain pooled estimates of variance components, and
`df.com`

can be used to specify the complete-data degrees of
freedom, which provides more appropriate (i.e., conservative) inferences
in smaller samples.

For example, using a conservative value for the complete-data degrees of freedom for the fixed effects in the model of interest (Snijders & Bosker, 2012), the output changes as follows.

`testEstimates(fit, extra.pars = TRUE, df.com = 46)`

```
#
# Call:
#
# testEstimates(model = fit, extra.pars = TRUE, df.com = 46)
#
# Final parameter estimates and inferences obtained from 20 imputed data sets.
#
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 433.015 28.481 15.203 36.965 0.000 0.153 0.176
# SexGirl 3.380 7.335 0.461 43.752 0.647 0.008 0.051
# I.SES 0.692 0.257 2.690 27.781 0.012 0.399 0.332
# G.SES 1.296 0.597 2.173 37.022 0.036 0.152 0.175
#
# Estimate
# Intercept~~Intercept|ID 168.506
# Residual~~Residual 8092.631
# ICC|ID 0.020
#
# Hypothesis test adjusted for small samples with df=[46]
# complete-data degrees of freedom.
```

Oftentimes, statistical inference concerns more than one parameter at a time. For example, the combined influence of SES (within and between groups) on mathematics achievement is represented by two parameters in the model of interest.

Multiple pooling methods for Wald and likelihood ratio tests (LRTs)
are implemented in the `testModels`

function. This function
requires the specification of a full model and a restricted model, which
are then compared using (pooled) Wald tests or LRTs. Specifically,
`testModels`

allows users to pool Wald tests (\(D_1\)), \(\chi^2\) test statistics (\(D_2\)), and LRTs (\(D_3\) and \(D_4\); for a comparison of these methods,
see also Grund, Lüdtke, & Robitzsch, 2016b).

To examine the combined influence of SES on mathematics achievement, the following restricted model can be specified and compared with the model of interest (using \(D_1\)).

```
<- with(implist, {
fit.null lmer(MathAchiev ~ 1 + Sex + (1|ID))
})
testModels(fit, fit.null)
```

```
#
# Call:
#
# testModels(model = fit, null.model = fit.null)
#
# Model comparison calculated from 20 imputed data sets.
# Combination method: D1
#
# F.value df1 df2 P(>F) RIV
# 6.095 2 674.475 0.002 0.275
#
# Unadjusted hypothesis test as appropriate in larger samples.
# Models fitted with REML were used as is.
```

Note regarding the order of arguments:Please note that`testModels`

expects that the first argument contains the full model, and the second argument contains the restricted model. If the order of the arguments is reversed, the results will not be interpretable.

Similar to the test for individual parameters, smaller samples can be
accommodated with `testModels`

(with method \(D_1\)) by specifying the complete-data
degrees of freedom for the denominator of the \(F\) statistic.

`testModels(fit, fit.null, df.com = 46)`

```
#
# Call:
#
# testModels(model = fit, null.model = fit.null, df.com = 46)
#
# Model comparison calculated from 20 imputed data sets.
# Combination method: D1
#
# F.value df1 df2 P(>F) RIV
# 6.095 2 39.812 0.005 0.275
#
# Hypothesis test adjusted for small samples with df=[46]
# complete-data degrees of freedom.
# Models fitted with REML were used as is.
```

The pooling method used by `testModels`

is determined by
the `method`

argument. For example, to calculate the pooled
LRT corresponding to the Wald test above (i.e., \(D_3\)), the following command can be
issued.

`testModels(fit, fit.null, method="D3")`

```
#
# Call:
#
# testModels(model = fit, null.model = fit.null, method = "D3")
#
# Model comparison calculated from 20 imputed data sets.
# Combination method: D3
#
# F.value df1 df2 P(>F) RIV
# 5.787 2 519.143 0.003 0.328
#
# Models originally fitted with REML were refitted using ML.
```

Finally, it is often useful to investigate functions (or constraints)
of the parameters in the model of interest. In complete data sets, this
can be achieved with a test of linear hypotheses or the delta method.
The `mitml`

package implements a pooled version of the delta
method in the `testConstraints`

function.

For example, the combined influence of SES on mathematics achievement
can also be tested without model comparisons by testing the constraint
that the parameters pertaining to `I.SES`

and
`G.SES`

are both zero. This constraint is defined and tested
as follows.

```
<- c("I.SES", "G.SES")
c1 testConstraints(fit, constraints = c1)
```

```
#
# Call:
#
# testConstraints(model = fit, constraints = c1)
#
# Hypothesis test calculated from 20 imputed data sets. The following
# constraints were specified:
#
# Estimate Std. Error
# I.SES: 0.692 0.245
# G.SES: 1.296 0.628
#
# Combination method: D1
#
# F.value df1 df2 P(>F) RIV
# 6.095 2 674.475 0.002 0.275
#
# Unadjusted hypothesis test as appropriate in larger samples.
```

This test is identical to the Wald test given in the previous section. Arbitrary constraints on the parameters can be specified and tested in this manner, where each character string denotes an expression to be tested against zero.

In the present example, we are also interested in the
*contextual* effect of SES on mathematics achievement (e.g.,
Snijders & Bosker, 2012). The contextual effect is simply the
difference between the coefficients pertaining to `G.SES`

and
`I.SES`

and can be tested as follows.

```
<- c("G.SES - I.SES")
c2 testConstraints(fit, constraints = c2)
```

```
#
# Call:
#
# testConstraints(model = fit, constraints = c2)
#
# Hypothesis test calculated from 20 imputed data sets. The following
# constraints were specified:
#
# Estimate Std. Error
# G.SES - I.SES: 0.605 0.644
#
# Combination method: D1
#
# F.value df1 df2 P(>F) RIV
# 0.881 1 616.380 0.348 0.166
#
# Unadjusted hypothesis test as appropriate in larger samples.
```

Similar to model comparisons, constraints can be tested with
different methods (\(D_1\) and \(D_2\)) and can accommodate smaller samples
by a value for `df.com`

. Further examples for the analysis of
multiply imputed data sets with `mitml`

are given by Enders
(2016) and Grund, Lüdtke, and Robitzsch (2016a).

Enders, C. K. (2016). Multiple imputation as a flexible tool for
missing data handling in clinical research. *Behaviour Research and
Therapy*. doi: 10.1016/j.brat.2016.11.008 (Link)

Grund, S., Lüdtke, O., & Robitzsch, A. (2016a). Multiple
imputation of multilevel missing data: An introduction to the R package
pan. *SAGE Open*, *6*(4), 1–17. doi:
10.1177/2158244016668220 (Link)

Grund, S., Lüdtke, O., & Robitzsch, A. (2016b). Pooling ANOVA
results from multiply imputed datasets: A simulation study.
*Methodology*, *12*, 75–88. doi: 10.1027/1614-2241/a000111
(Link)

Rubin, D. B. (1987). *Multiple imputation for nonresponse in
surveys*. Hoboken, NJ: Wiley.

Snijders, T. A. B., & Bosker, R. J. (2012). *Multilevel
analysis: An introduction to basic and advanced multilevel
modeling*. Thousand Oaks, CA: Sage.

```
# Author: Simon Grund (simon.grund@uni-hamburg.de)
# Date: 2023-01-05
```