## 7.4 Multilevel imputation by joint modeling

The joint model specifies a single model for all incomplete variables in data. Suppose that we have missing values in both outcome $$y_c$$ and the level-1 predictors $$X_c$$ in the linear mixed model (7.6) with random intercepts. A limitation of the standard model is that both $$X_c$$ and $$Z_c$$ need to be completely observed. Schafer and Yucel (2002) showed that it is possible to reformulate the model by placing all variables in the model as an outcome on the left-hand side of model (7.6), which gives rise to multilevel models with multivariate outcomes. Suppose that $$Y_{ic}$$ is the $$1 \times p$$ matrix of all incomplete variables from unit $$i$$ in class $$c$$. Under a random intercept model, the multivariate multilevel model decomposes the data into a between-group part and a within-group part as

$Y_{ic} = \mu + Y_c^\mathrm{between} + Y_{ic}^\mathrm{within}. \tag{7.7}$

The relations between the $$p$$ variables at the group level are captured by the between component $$Y_c^\mathrm{between}$$, while relations at the individual level are captured by $$Y_{ic}^\mathrm{within}$$. Lüdtke, Robitzsch, and Grund (2017) explain how this decomposition can be applied to draw imputations for the missing elements in $$Y_c$$.

The approach has been implemented in the pan package (Zhao and Schafer 2016). and Yucel (2011) discussed extensions for categorical data using the generalized linear mixed model. Goldstein et al. (2009) described a joint model for mixed continuous-categorical data with a multilevel structure. Carpenter and Kenward (2013) and Goldstein, Carpenter, and Browne (2014) proposed extensions for ordinal and unordered categorical data, which are implemented in the REALCOM-IMPUTE software (Carpenter, Goldstein, and Kenward 2011). Quartagno and Carpenter (2016) extended this work by allowing for heteroscedasticity of the imputation model, combined with imputation of binary (and more generally categorical) variables, which is available through the jomo package (Quartagno and Carpenter 2017). Related work has been done by Asparouhov and Muthén (2010) for Mplus.

Because of the decomposition (7.7), imputation by joint modeling is very natural when the complete-data analysis focuses on different within- and between-cluster associations (Enders, Mistler, and Keller 2016; Grund, Lüdtke, and Robitzsch 2016a). Multilevel imputation is not without problems. Except for jomo, most models assume a homoscedastic error structure in the level-1 residuals, which implies no random slope variation between $$Y_{ic}$$ (Carpenter and Kenward 2013; Enders, Mistler, and Keller 2016). Imputations created by jomo reflect pairwise linear relationships in the data and ignore higher-orders interaction and non-linearities. Joint modeling may also experience difficulties with smaller samples and the default inverse Wishart prior (Grund, Lüdtke, and Robitzsch 2018b; Audigier et al. 2018).

Imputation models can also be formulated as latent class models (Vermunt et al. 2008; Vidotto, Vermunt, and Kaptein 2015). proposed a Bayesian multilevel latent class model that is designed to capture heterogeneity in the data at both levels through local independence and conditional independence assumptions. This class of models is quite flexible. As the method is very recent, there is not yet much practical experience.