Analysis of Multiply Imputed Data Sets

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.

Example data (`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.

Generating imputations

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

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

The completed data are then extracted with mitmlComplete.

implist <- mitmlComplete(imp, "all")

Transforming the imputed data sets

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.

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

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.
implist <- with(implist,{
  df <- data.frame(as.list(environment()))
  df <- ... # dplyr commands
  df
})
implist <- as.mitml.list(implist)
Advanced users may also consider using lapply for a similar workaround.`

Fitting the analysis model

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.

fit <- with(implist, {
  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.

Pooling

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.

Parameter estimates

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.

Multiple parameters and model comparisons

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₁), χ² test statistics (D₂), and LRTs (D₃ and D₄; 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₁).

fit.null <- with(implist, {
  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₁) 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₃), 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.

Constraints on parameters

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.

c1 <- c("I.SES", "G.SES")
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.

c2 <- c("G.SES - I.SES")
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₁ and D₂) 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).

References

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 ([email protected])
# Date:   2024-11-02