Nested data

In nested data, every observation belongs to one and only one higher-level unit. Two-level data are a simple example for this type data, as shown below for six clusters with \(n\) observations.

More deeply nested structures are possible. For example, in three-level data, the clusters themselves are nested in even-higher-level units (e.g., students nested in classrooms nested in schools).

Here, observations \(p\) are nested within clusters \(C\), and clusters are nested within higher-level clusters \(D\).

Cross-classified data

In cross-classified data, every observation belongs directly to two or more clusters at once (e.g., experimental data with observations clustered within subjects and stimuli). However, the clusters are not themselves nested within one another but “crossed” as shown below.

In contrast to nested data, there is no clear hierarchy of the two cluster variables. Differently put, although both \(A\) and \(B\) have observations clustered within them, neither of the two is itself nested within the other.

Why bother?

For the treatment of missing data, the hierarchical structure must be taken into account when using model-based methods such as MI (Enders, Mistler, & Keller, 2016; Lüdtke, Robitzsch, & Grund, 2017). This means that we need to acknowledge that, in hierarchical data, variables can vary both within and between clusters, and multiple variables can be related at each level of the structure.

Several articles have considered the case with two-level data (e.g., the two above). In the following, I show two examples for how to conduct MI for three-level and cross-classified data in R.

Example data

For this example, I simulated data with a three-level structure consisting of 50 schools (level 3), five classrooms per school (level 2), and 10 students per classroom (level 1, total \(n\) = 2500). The data can be downloaded here (CSV, Rdata).

Every row corresponds to one student, and the classrooms and schools are numbered consecutively. All variables in the data set contain missing data.

Multiple imputation

To perform MI, I use the R packages mice and miceadds. The mice package treats missing data by iterating through a sequence of imputation models, thus treating variable after variable in a step-by-step manner (for a general introduction to mice, see van Buuren & Groothuis-Oudshoorn, 2011).

The imputation models are set up by defining (a) a method for each variable, naming the type of model to be used, and (b) a predictor matrix, naming which predictors (columns) should be used for each variable (rows). Extracting the defaults provides a good starting point.

library(mice)
library(miceadds)

# predictor matrix and imputation method (defaults)
predMatrix <- make.predictorMatrix(data = dat)
impMethod <- make.method(data = dat)

By default, mice uses methods intended for non-hierarchical data. For multilevel data, we need to ensure that the imputation model takes the multilevel structure into account such that the models will need to include variance components at higher levels and allow for different relations between variables at different levels.

# method for lower-level variables (x, y, and z)
impMethod[c("x", "y", "z")] <- "ml.lmer"

# method for variables at top level (w)
impMethod["w"] <- "2lonly.norm"

# remove indicator variables from predictor matrix
predMatrix[, c("class", "school")] <- 0

# ... specify cluster indicator (2lonly.norm)
predMatrix["w", "school"] <- -2

# ... specify levels of higher-level variables
level <- character(ncol(dat))
names(level) <- colnames(dat)

level["w"] <- "school"
level["z"] <- "class"

# ... specify cluster indicators (as list)
cluster <- list()

cluster[["x"]] <- c("class", "school")
cluster[["y"]] <- c("class", "school")
cluster[["z"]] <- c("school")

# run mice
imp <- mice(dat, method = impMethod, predictorMatrix = predMatrix, maxit = 20,
            m = 20, levels_id = cluster, variables_levels = level)

Example analysis

I use the R packages mitml and lme4 to analyze the imputed data. First, I extract a list imputed data sets and calculate the cluster means that we need in order to fit the analysis model.

library(mitml)

# create list of imputed data sets
implist <- mids2mitml.list(imp)

# calculate group means
implist <- within(implist, {
  x.cls <- clusterMeans(x, class)
  x.sch <- clusterMeans(x, school)
})

The analysis model is then fitted with the lme4 package, and the results are pooled with mitml with the following lines of code.

library(lme4)

# fit model
fit <- with(implist,{
  lmer(y ~ 1 + x + x.cls + x.sch + z + w + (1|class) + (1|school))
})

# pool results
testEstimates(fit, var.comp = TRUE)

# 
# Call:
# 
# testEstimates(model = fit, var.comp = TRUE)
# 
# Final parameter estimates and inferences obtained from 20 imputed data sets.
# 
#              Estimate Std.Error   t.value        df   P(>|t|)       RIV       FMI 
# (Intercept)    -0.025     0.082    -0.306  3462.095     0.760     0.080     0.075 
# x               0.190     0.019    10.186   113.399     0.000     0.693     0.419 
# x.cls           0.460     0.058     7.968   220.257     0.000     0.416     0.300 
# x.sch           0.352     0.149     2.359  3807.532     0.018     0.076     0.071 
# z              -0.135     0.028    -4.774   363.762     0.000     0.296     0.233 
# w              -0.071     0.083    -0.860   446.933     0.390     0.260     0.210 
# 
#                             Estimate 
# Intercept~~Intercept|class     0.084 
# Intercept~~Intercept|school    0.247 
# Residual~~Residual             0.312 
# 
# Unadjusted hypothesis test as appropriate in larger samples.

These results are very close to the parameters I used to generate the data. In the next example, we move on to clustered data with a cross-classified structure.

Example data

For this example, I simulated data with a cross-classified structure and a total of \(n\) = 5000 responses (level 1) from 100 subjects (level 2b) to 50 items (level 2a). The experimental condition (\(a\) = 1) comprised all even-numbered items; the control (\(a\) = 0) all the odd-numbered items. The data can be downloaded here (CSV, Rdata).

Some of the responses (\(y\)) are sporadically missing. In addition, a number of subjects failed to provide data on the subject-level covariate (\(z\)).

Multiple imputation

The main strategy for MI remains the same as in the previous example. In order to accommodate the multilevel structure, we again need to ensure that the imputation model allows for the variables to have variance and relations with each other at different levels, with the exception that aggregated variables are not used here (see Footnote 1).

We again start with the default setup and adjust it the way we need to.

# create default predictor matrix and imputation methods
predMatrix <- make.predictorMatrix(data = dat)
impMethod <- make.method(data = dat)

# method for lower-level variables (y)
impMethod["y"] <- "ml.lmer"

# ... for variables at top level (z)
impMethod["z"] <- "2lonly.norm"

# remove indicator variables from set of predictors
predMatrix[, c("subject", "item")] <- 0

# specify cluster indicator (2lonly.norm)
predMatrix["z", "subject"] <- -2

# specify levels of higher-level variables
level <- character(ncol(dat))
names(level) <- colnames(dat)

level["a"] <- "item"
level["z"] <- "subject"

# specify cluster indicators (as list)
cluster <- list()

cluster[["y"]] <- c("subject", "item")

# remove effect of (average) a on z
predMatrix["z", "a"] <- 0

# run mice
imp <- mice(dat, method = impMethod, predictorMatrix = predMatrix, maxit = 20, 
            m = 20, levels_id = cluster, variables_levels = level,
            aggregate_automatically = FALSE)

Example analysis

The analysis of the data is done with lme4 and mitml as before. First, we extract the imputed data sets as a list.

# create list of imputed data sets
implist <- mids2mitml.list(imp)

Then we fit the analysis model with lme4 and pool the results with mitml.

# fit model
fit <- with(implist,{
  lmer(y ~ 1 + a + z + (1|item) + (1|subject))
})

# pool results
testEstimates(fit, var.comp = TRUE)

# 
# Call:
# 
# testEstimates(model = fit, var.comp = TRUE)
# 
# Final parameter estimates and inferences obtained from 20 imputed data sets.
# 
#                Estimate   Std.Error     t.value          df     P(>|t|)         RIV         FMI 
# (Intercept)      -0.161       0.106      -1.514  154490.706       0.130       0.011       0.011 
# a                 0.659       0.124       5.324 2581059.039       0.000       0.003       0.003 
# z                -0.248       0.072      -3.469     361.838       0.001       0.297       0.233 
# 
#                              Estimate 
# Intercept~~Intercept|subject    0.341 
# Intercept~~Intercept|item       0.186 
# Residual~~Residual              0.462 
# 
# Unadjusted hypothesis test as appropriate in larger samples.

The results are again close to the true values I used to generate the data.