Quantile Regression-Based Multiple Imputation of Missing Values — An Evaluation and Application to Corporal Punishment Data

Multiple Imputation

Multiple imputation based on Rubin’s (1987) theory is a three-step process: First, each missing value is imputed multiple times by plausible values, so that researchers end up with $M$ complete datasets. Plausible in this context refers to any value that leads to unbiased statistical inferences.¹ Second, these $M$ datasets are analysed separately by any complete data technique such as OLS regression, or QR and yield $M$ sets of modelling results. In the final step, the $M$ sets of modelling results are combined into an overall result based on combination rules such as the large-sample formula given in Rubin (1987) or adaptations for small samples (Barnard & Rubin, 1999). For a detailed discussions, see Kleinke, Reinecke, Salfrán et al. (2020, Chapter 4). The combination formulae yield point estimates and standard errors for each model parameter ${\mathit{\theta }}_p$, $p=1,\dots ,P$. The combined point estimate is simply the average of the $M$ estimates. Corresponding standard errors are obtained by combining a within-imputation and a between-imputation variance component. The latter component reflects the amount of missing information with respect to ${\mathit{\theta }}_p$: Fewer information to predict the unobserved values leads to higher uncertainty in the imputations and a larger between variance. Note that if the missing values could be predicted without any uncertainty, the estimator ${\widehat{\mathit{\theta }}}_{p,m}$ obtained from the $m=1,\dots ,M$ imputed datasets would not vary over the datasets – resulting in a between variance of zero.

More comprehensive introductions to MI and helpful resources for applied researchers are given in: Kleinke, Reinecke, Salfrán et al. (2020), who provide a general introduction and review recent developments and various applications of MI. Chapters 1, 4, and 5 in Little and Rubin (2020) summarize the problem of missing values, missing value patterns and mechanisms; additionally, a taxonomy of methods for imputing missing values and a discussion of the rationale of MI in comparison to other methods to account for uncertainty in parameter estimates due to missing values are given. Applications of MI in various contexts including linear and generalized linear models, categorical data analysis, survival analysis, structural equation modelling, longitudinal data analysis, and complex survey settings are discussed in Raghunathan et al. (2018).

QR and QR-Based MI of Missing Values

The QR model assumes a linear additive functional form for the $\tau$-quantiles of the response $y_i$ conditional on the $P$ predictors ${\mathit{x}}_i$, where $i=1,\dots ,n$ are the sample observations. Define the conditional $\tau$ quantile as

where ${\mathit{\beta }}_{\tau }$ denotes the parameter corresponding to the $\tau$th quantile. Parameter estimates ${\widehat{\mathit{\beta }}}_{\tau }$ are obtained by minimizing the asymmetrically weighted loss function where ${\rho }_{\tau }\left(\cdot \right)$ governs the weighting. Choosing $\tau =0.5$ yields the conditional median (for further details, see Koenker, 2005).

The following data generating process illustrates the exact fit property of QR:

where $F_{\tau ,i}$ is some distribution function and $F_{\tau ,i}^{-1}$ its inverse. Note that if $q$ is the conditional $\tau$ quantile, such that $\tau =F_{\tau ,i}\left(q\right)$, it holds that $q=F_{\tau ,i}^{-1}\left(\tau \right)$ and the conditional $\tau$ quantile can be recovered by inverting the distribution function at $\tau$. The weak assumptions about the distribution of the errors and the fact that quantiles of the conditional distribution of $y_i$ are modelled implies advantages over restrictive parametric methods when distributional assumptions of the latter are violated (e.g., in case of outliers or heteroscedasticity). A disadvantage of QR is that a sufficiently large sample size is required for the model parameters to be identified because all parameters are estimated for each quantile.²

To illustrate QR-based MI of missing values, denote the response vector by $\mathit{y}$, the matrix of $P$ predictors by $\mathit{X}$, and assume that all variables contain missing values (or incomplete cases). Define the dataset with missing values by $\mathit{U}=\left(\mathit{y},\mathit{X}\right)$, where each of the $i=1,\dots ,n$ lines of the dataset is at least partially observed. The columns of the dataset $\mathit{U}$ are referred to by $q=1,\dots ,Q$, with $Q=P+1$. Note that the total number of observations in the dataset is $N=n\cdot Q$. Denote the total number of missing values as $N_{mis}$ and the number of observed values (or complete cases) as $N_{obs}=N-N_{mis}$. The number of missing values in each column of the dataset is referred to as $n_{q,mis}$, such that $N_{mis}=n_{1,mis}+\dots +n_{Q,mis}$.

We aim to impute the missing values in $\mathit{U}$ and use them in the analysis model. The matrix of imputed values is represented by $\mathit{V}$. To simplify the exposition, assume that $\mathit{U}$ is sorted such that the lines with the observed values are given first and refer to the individual lines via ${\mathit{u}}_j$, $j=1,\dots ,J$ ( $J=min\left(n_{1,obs},\dots ,n_{Q,obs}\right)$). The remaining lines with the missing values are denoted by ${\mathit{u}}_k$, $k=1,\dots ,K$ ( $K=N-J$). The definition of ${\mathit{v}}_j$ and ${\mathit{v}}_k$ is analogous.

In the following, we describe function mice.impute.rq from package Qtools for QR-based MI in greater detail. Note that we do not consider approaches restricting regression quantiles to avoid quantile crossing. Based on the defined notation, function mice.impute.rq proceeds according to the Algorithm 1 (Geraci, 2016a).

Algorithm 1 details the steps required for QR-based MI of missing values. First, choose the number of iterations $C$ and the number of bootstrap samples $M$. For bootstrap iteration $m=1,\dots ,M$, proceed as follows: Define a temporary dataset ${\mathit{D}}_{m,c}$; in the first iteration $c=1$, the lines of the dataset without missing values ${\mathit{u}}_j$ are used. At subsequent iterations $2,\dots ,C$, the imputed values from the previous iteration ${\mathit{V}}_{m,c-1}$ are employed and the following five steps are used on all columns $q=1,\dots ,Q$: In (1), draw a bootstrap sample ${\mathit{B}}_{m,c}$ of size $n_{q,mis}$ from the lines of ${\mathit{D}}_{m,c}$ with replacement. In (2), partition ${\mathit{B}}_{m,c}$ into its $q$ th column ${\mathit{b}}_{m,c,q}$ and all other columns ${\mathit{B}}_{m,c,-q}$. In (3), draw $K$ random values ${\tau }_1,\dots ,{\tau }_K$ from the uniform distribution $\mathcal{U}\left(\epsilon ,1-\epsilon \right)$, with $0<\epsilon <0.5$ (note that the function default is $\epsilon =0.001$). In (4), fit a QR model based on the partitioned bootstrap sample $Q_{b_{i,q}\mid {\mathit{b}}_{i,-q}}\left(\tau \right)={\mathit{b}}_{i,-q}^T{\mathit{\beta }}_{\tau }$ based on the bootstrap sample ${\mathit{B}}_{m,c}$. In (5), use the parameter estimates ${\mathit{\beta }}_{\tau }$ from (4) to impute the missing values for each column of the dataset with missing values ${\mathit{U}}_K$. Collect the imputed values for all $q=1,\dots ,Q$ columns and assign them to ${\mathit{V}}_m$. Then, iterate the five steps until the pre-specified number of iterations $C$ is reached.³ The algorithm yields $M$ datasets with imputed values ${\mathit{V}}_1,\dots ,{\mathit{V}}_M$.

GAMLSS and GAMLSS-Based MI of Missing Values

De Jong (2012) and de Jong et al. (2016) have proposed a MI procedure based on a class of semi-parametric GAMLSS (Rigby & Stasinopoulos, 2005; Stasinopoulos & Rigby, 2007), which is implemented in R package ImputeRobust (Salfrán & Spiess, 2018) and extends the functionality of the mice package. Compared to standard normal model-based MI, which only considers modelling the location of the distribution, GAMLSS-based MI also allows to account for the scale, and the shape of the distribution and can be better tailored to a given empirical dataset.

In the following, we slightly extend the notation from the previous subsection to illustrate GAMLSS-based MI briefly.⁴ Assume that the conditional distribution of the $q$th column ${\mathit{u}}_q$ of dataset $\mathit{U}$ given the observed values ${\mathit{u}}_{J,q}$ and the other columns ${\mathit{U}}_{-q}$ is

Parameters ${\mathit{\theta }}^{\left(1\right)}$ and ${\mathit{\theta }}^{\left(2\right)}$ are the location and scale parameters, and ${\mathit{\theta }}^{\left(3\right)},\dots ,{\mathit{\theta }}^{\left(S\right)}$ are shape parameters – depending on the specified distribution $\mathcal{D}$. Imputations for the missing values in column $q$ of the dataset ${\mathit{U}}_K$ can be obtained by where ${\mathit{\theta }}_q$ is a parameter vector and $g_q\left(\cdot \right)$ a link function relating the location, scale, and shape parameter of the missing values in column $q$ to the dataset $\mathit{U}$. Equation 5 is additively decomposed into a linear parametric part with design matrix ${\mathit{U}}_q$ and corresponding coefficient vector ${\mathit{\gamma }}_q$ and a sum of smooth functions $h_{j,q}\left(\cdot \right)$ of the observed values ${\mathit{u}}_{J,q}$.

Since package gamlss provides different approaches to model $h\left(\cdot \right)$, such as cubic smoothing splines, penalized splines, and local regression, the following algorithm provides a general description of GAMLSS-based MI (Salfrán & Spiess, 2018).

Algorithm 2 details GAMLSS-based MI of missing values. First, choose the number of bootstrap samples $M$, set $m=1$, use the lines of the dataset without missing values ${\mathit{U}}_J$ as initial bootstrap sample ${\mathit{B}}_1$, and choose starting values ${\mathit{\theta }}_{0,k,q}^{\left(1\right)},{\mathit{\theta }}_{0,k,q}^{\left(2\right)},{\mathit{\theta }}_{0,k,q}^{\left(3\right)},\dots ,{\mathit{\theta }}_{0,k,q}^{\left(S\right)}$ for the parameters of the distribution $\mathcal{D}$. In (1), the penalized likelihood function ${\mathit{\ell }}_{\text{p}\text{e}\text{n}}$ is specified based on the distribution $\mathcal{D}$. In (2), maximise ${\mathit{\ell }}_{\text{p}\text{e}\text{n}}$ based on ${\mathit{B}}_1$ to obtain parameter estimates ${\widehat{\mathit{\theta }}}_{1,k,q}^{\left(1\right)},{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(2\right)},{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(3\right)},\dots ,{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(S\right)}$ and fitted values for column $q$ of the dataset ${\widehat{\mathit{u}}}_{J,q}$. Next, set $m=2$. In (3), the bootstrap sample ${\mathit{B}}_m$ is updated by drawing $n_{q,mis}$ observations with replacement from the fitted values from the previous bootstrap iteration ${\widehat{\mathit{u}}}_{m-1,q}$ and the corresponding observed columns of the dataset besides the $q$th column ${\mathit{U}}_{J,-q}$. In (4), maximise ${\mathit{\ell }}_{\text{p}\text{e}\text{n}}$ based on ${\mathit{B}}_m$ and update the parameter estimates ${\widehat{\mathit{\theta }}}_{1,k,q}^{\left(1\right)},{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(2\right)},{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(3\right)},\dots ,{\widehat{\mathit{\theta }}}_{1,k,q}^{\left(S\right)}$. Then, repeat (3) and (4) until $m=M$. The algorithm yields $M$ datasets with imputed values ${\mathit{V}}_1,\dots ,{\mathit{V}}_M$.

Imputation Based on PMM

Another approach that was recommended for a wide range of scenarios is PMM (cf. van Buuren, 2018, Chapter 3): The method samples from a pool of observed cases, which are ‘similar’ to the incomplete case. This so-called donor is then used to fill in the corresponding missing value of the donee. The similarity of observations is defined by similarity in predicted means for the observed and unobserved part of the variable based on a linear regression model. Details about the implementation in package mice (function mice.impute.pmm) are given in van Buuren (2018, Chapter 3).

PMM is considered a semi-parametric imputation procedure: The parametric linear model is used only as a means for matching the incomplete case to a potential donor. Then, actually observed values are imputed. This preserves many properties of the observed part of the data, which makes PMM robust against mild to moderate violations of distributional assumptions of the normal linear regression model. This was systematically evaluated by Kleinke (2017) for different degrees of skewness, different sample sizes, and varying percentages of missing values. Their results indicated that: 1) especially when sample size was small, small donor pools yielded better results in terms of unbiased parameter estimates and standard errors compared to large donor pools; 2) PMM generally worked well, unless the data were highly skewed and/or more than roughly 20% to 30% of the data were imputed; 3) PMM generally performed better when sample size was large.

Disadvantages of PMM are that the approach can neither extrapolate beyond the range of the observed data, nor interpolate when regions of the data are sparsely populated. Additionally, PMM fails when there are substantial violations of the assumptions of the normal linear regression model.

Overview

In this section, we evaluate QR-based MI by Monte Carlo Simulation based on a typical setting in the social sciences: Skewed and heteroskedastic data. The underlying data were taken from Haupt et al. (2014), who investigated risk factors for corporal punishment by fathers with a QR model. We proceed as follows: First, draw 1000 bootstrap samples with replacement from the complete cases and fit a QR model to each of the bootstrap samples. Second, use the average parameter estimates as a proxy for the unknown ‘true’ population parameters. Third, introduce that values are missing at random (MAR) to the bootstrap samples by deleting a certain fraction of values. Finally, evaluate QR-based MI regarding bias in point estimates and coverage rates (CR) of the 95% confidence intervals of the parameters and contrast the results with the complete data estimates and results obtained from other established MI methods (standard normal model-based MI, PMM, and GAMLSS-based MI).

Hypotheses

The advantages and disadvantages of QR in comparison to GAMLSS are discussed in Kneib (2013) and Rigby et al. (2013). QR, for example, is distribution-free, while GAMLSS introduce model structure and are more parsimonious regarding the data requirements – while also bearing a higher risk of mis-specification. To our best knowledge, QR- and GAMLSS-based MI have not been compared systematically in Monte Carlo simulations, yet. We expect the preferred method to depend on the specific setting. In our simulation, we expect QR- and GAMLSS-based MI to yield comparable results. We assume that both approaches will outperform PMM and standard normal model-based MI – especially when quantiles farther away from the center of the distribution of the response are considered: PMM and normal model-based MI create imputations based on predictions for the conditional mean. Thus, both normal model-based MI and PMM should perform well in the vicinity of the center of the distribution.

Design

Our simulation has a $2\times 2\times 4$ (percentage of missing values $\times$ missing value pattern $\times$ imputation method) design: From the bootstrap samples, we deleted 25% or 40% of values either from the response only, or from the response and one of the predictors following a MAR pattern (for details, see Section “Procedure”). We imputed missing values by QR- and GAMLSS-based MI, PMM, and standard normal model-based MI.

Evaluation Criteria

We evaluated the approaches based on the goal to yield unbiased statistical inference: Bias in the parameter estimate was defined as the difference between the true parameter $\mathit{\theta }$ and the average estimate $\hat{\mathit{\theta }}$ across the bootstrap samples. Bias was reported relative to the true population parameter as

Coverage Rate (CR) was defined as the percentage of 95% confidence intervals that included the true population parameter. Schafer and Graham (2002) deem CR below 90% as seriously low. Undercoverage may result from biased point estimates and/or an underestimated variance. The confidence interval could either be too narrow to include the true parameter, and/or it can be shifted too much to the left or to the right to include the true parameter. Together, bias and CR allow to assess the accuracy of statistical inferences.

Sample and Measures

Our simulation is based on a sample taken from the longitudinal Erlangen-Nuremberg Development and Prevention Study. The data include 609 families (675 preschool children) from the North-Bavarian cities of Erlangen and Nuremberg (for more details, see Lösel et al., 2009).

The response in our model corporal punishment by fathers (CP), as well as two risk factors inconsistent discipline (ID) and other discipline practices (OD) were assessed using a German version of the Alabama Parenting Questionnaire (APQ; Lösel et al., 2003). Items were rated on five-point scales ranging from 1 (“never”) to 5 (“always”). Scale scores were obtained by calculating the means of the respective items. The corporal punishment scale was based on three items such as “You spank your child with your hand when he/she has done something wrong”. The empirical distribution of the corporal punishment scale is summarized in Figure 1.⁵ Inconsistent discipline (six items) summarizes parenting behaviors like “The punishment you give your child depends on your mood”. Other discipline practices (seven items) refer to actions such as “You take away privileges or money from your child as a punishment”. A further risk factor of interest was socioeconomic status (SES). The computation of the SES index was based on Geißler (1994). Finally, our model included the child’s gender by the gender indicator (FEMALE).

Procedure

We simulated 1000 bootstrap samples from the original $N=485$ complete cases. Each bootstrap sample had a size corresponding to the original sample size $N=675$. We then introduced MAR patterns (Rubin, 1976). Depending on the respective condition, we deleted either 25% or 40% of the values of the response CP, and/or the predictor OD. Values were deleted using the following rule: For each observation $i$ we computed the probability of a missing value $P{r}_i^{mis}$ with the logit model

Coefficients in Equation 7 were chosen to ensure probabilities of missing values ranging from about 0.2 to 0.8 – depending on the individual SES scores. In our simulation, participants from lower social classes were assigned larger probabilities of missing values.

Depending on the respective condition, we sampled the necessary number of cases with selection probabilities $P{r}_i^{mis}$ and deleted the respective observations from the datasets. An overview of the resulting missing value patterns is displayed in Table 1.

Missing values were imputed with QR- and GAMLSS-based MI, PMM, and normal model-based MI by using functions implemented in packages mice, Qtools, and ImputeRobust. Unless stated otherwise, we used default settings of the functions. All other columns of the dataset were used to predict missing values in a particular column. As common practice (e.g., Kleinke, 2018), the number of imputations was set to the percentage of missing values.

We then compared parameter estimates based on the imputed datasets with the means of the 1000 complete data estimates. We obtained coefficient estimates from a QR model, in which we regressed the response CP on SES, FEMALE, ID, and OD regarding quantiles $\tau \in \left\{.2,.3,.4,.5,.6,.7,.8\right\}$. Heteroscedasticity-robust standard errors were computed as in Haupt et al. (2014).

Results

Figure 2 summarizes the results on bias in the point estimates and displays the average amount of relative absolute bias across all model parameters. Following Forero and Maydeu Olivares (2009), bias was considered substantial if it exceeded 10% of the true population parameter. Overall, bias generally increased: (i) when the response and additionally predictor values were missing; (ii) with the number of missing values.

When considering the different approaches for MI, both PMM and normal model-based MI had only a small bias in the center part of the distribution of the response (i.e., $\tau \in \left\{.4,.5,.6\right\}$), on average. As expected, bias increased when moving towards the tails of the distribution for both approaches. QR- and GAMLSS-based MI produced widely unbiased results across all quantiles and Monte Carlo designs.

Figure 3 displays the corresponding CR. All approaches produced large enough confidence intervals (some even slighly above the nominal $\left(1-\alpha \right)$ percent), such that a bias in point estimates was alleviated by a wider interval. CR was sufficiently large across all parameter estimates and quantiles, even when the percentage of missing values was 40%.

Discussion

Normal model-based MI produced unbiased inferences in the center region of the response. While point estimates in the off-center regions were biased, CR were nevertheless above 90%. The observed behaviour of normal model-based MI is due to the so-called self-correcting property of MI: A mild to moderate bias in the parameter estimate resulting from the mis-specification of the imputation model by assuming normality and homoscedasticity of the error distribution is alleviated by an increase in the corresponding variance estimate. This also lead to accurate CR produced by PMM. This result should be of considerable interest to applied researchers facing datasets with incomplete cases: In practice, valid inference in off-center regions of the response are often of paramount interest – for example, when investigating risk factors for corporal punishment. Typically, applied researchers are not only interested in a suitable CR, but also in obtaining accurate point estimates regarding the tails of the distribution. To achieve this, both, the imputation and analysis model need to be able to capture the underlying properties of the data generating process. Our simulations show that this is widely the case for QR- and GAMLSS-based MI: Valid inference can be obtained from both approaches under different missing value patterns and percentages. Following previous recommendations (cf. Kleinke, Reinecke, Salfrán et al., 2020) that the imputation method and the subsequent analysis model should coincide, we recommend to use QR-based MI when the analysis model is a QR model.

Substantive Models

Haupt et al. (2014) considered further predictors in their QR model that we did not include in our simulations above. These were 1) fathers’ aggressiveness (AG), measured by the aggressiveness scale of the FPI-R personality questionnaire 2) parental distress (PD), measuring the perceived stress of being a parent, 3) dysfunctional interaction (DI), assessing the parents’ perception that the child does not meet their expectations or that they are manipulated by the child, and 4) the age of the father (AGE). For further details on the predictors, see Haupt et al. (2014). CP was regressed on the aforementioned predictors and heteroscedasticity-robust standard errors were computed (Powell, 1991). Regression coefficients were estimated for quantiles $\tau \in \left\{.2,.3,.4,.5,.6,.7,.8\right\}$.

In addition to Haupt et al. (2014), Fritsch et al. (2019) explored the relevance of further predictors (for an overview and a detailed description, see Table 1 and Section 3 in Fritsch et al., 2019). These variables were age difference of the parents (AGE.DIFF), as well as categorical indicators if the mother or father were older than the other parent respectively, parental involvement, positive parenting, monitoring and supervision by the parents, a measure of the ‘difficulty’ of the child, fathers’ life satisfaction, social orientation, performance orientation, excitability, inhibition, levels of stress, health concerns and physical complaints, candor, extraversion, and emotionality. Fritsch et al. (2019) focused on a prediction-centered approach and modelled the conditional mean of CP. We considered changes in the predictive performance when using the dataset obtained from QR-based MI.

Missing Values

Complete case analysis can be expected to yield unbiased statistical inferences in certain setting only, for example when only values in the response are missing or when the missing values are missing completely at random (MCAR) in the sense of Rubin (1976). MCAR means the probability of a missing value in a variable neither depends on observed information in the dataset (e.g., on SES), nor on unobserved information (e.g., on the unobserved CP values). The MCAR assumption is violated for example, when variables in the dataset have non-zero correlations with the missing value indicators of incompletely observed variables. We computed correlations of each variable with the missing value indicators of other variables in the dataset and found only small correlations between –.16 and .15. Effects of MAR should therefore be rather minor and we expect only marginal differences from the results reported in Haupt et al. (2014).

For the present analysis we included data from 575 children from families where both parents were present at the beginning of the study and for whom we had complete information regarding AGE, AGE.DIFF, and SES. Missing value patterns (e.g., van Buuren & Groothuis-Oudshoorn, 2011) are displayed in Table 2. About 82% of all cases are complete cases. No answers about CP were available from 79 fathers (about 14%). As many of these fathers also did not answer the other questionnaires (e.g., about parental distress or their aggressiveness), we could use mainly sociodemographic information (e.g., SES and FEMALE) to impute the missing values about CP. For some cases, we could also make use of observed information from the Freiburg Personality Inventory.

We imputed missing values using package mice. Missing values in the response were imputed by function mice.impute.rq from R package Qtools (Geraci, 2016a, 2016b). All other variables were imputed using mice’s default settings. Imputation models were automatically generated using mice’s default settings (for details, see van Buuren, 2018; van Buuren & Groothuis Oudshoorn, 2011). We then ran 20 iterations of mice’s Gibbs sampler and created a total of $m=100$ imputations. We fitted the same QR model as in Haupt et al. (2014) to the imputed data and used standard combination rules (Barnard & Rubin, 1999; Rubin, 1987) to obtain the multiple imputation inferences.

MI Results

The top part of Table 3 displays point estimates and corresponding standard errors, the bottom part deviations from the complete case results reported in Haupt et al. (2014). Coefficient estimates and standard errors are similar and the qualitative results across the two analysis are identical: Predictors SES, FEMALE, and PD have negative coefficient estimates, which are close to zero for PD for small values of $\tau$. The strength of the absolute effects relative to the standard errors is highest for SES and smaller for FEMALE. For the latter predictor, the relative strength of the absolute effect decreases with increasing $\tau$. For PD, the relative strength of the effect increases with increasing $\tau$ until $\tau =0.5$ and slightly decreases for larger $\tau$. Positive coefficient estimates are obtained for all other predictors. The largest coefficient estimates relative to the standard errors are obtained for DI, ID, and OD. For AG and AGE, coefficient estimates are generally smaller or only slightly larger than the standard errors. This confirms the risk factors for CP identified by Haupt et al. (2014) for the dataset imputed by QR-based MI: Low SES, DI, ID, and OD – while the limitations of the cross-sectional analysis and the bidirectional design of the analysis noted by Haupt et al. (2014) remain unchanged.

Additionally, we reproduced the results of the prediction-centered modelling approach in Fritsch et al. (2019). The paper models the conditional expectation of the response CP with multiple linear regression, regression trees and random forests.⁶ An out-of-sample evaluation of the three approaches is carried out based on root mean square error, mean square error, and mean absolute error. The paper reports superior predictive performance of random forest-based modelling. We reproduced the analysis with the 100 datasets obtained from QR-based MI as follows: First we split the datasets into training and test data; for the test data, we set aside the exact same observations as in the original analysis across all imputed datasets. Second, we fitted a multiple linear regression model, a regression tree, and a random forest to each of the 100 imputed datasets. Note that for the regression trees and the random forests, we tuned the meta-parameters of the algorithm for each of the 100 datasets. Third, we predicted the response and evaluated the predictive performance of all methods by the three prediction error criteria for each of the 100 datasets. Overall, 1) there were only minor differences across the 100 datasets in predictive performance; 2) the predictive performance of modelling based on random forests and multiple linear regression improved slightly; 3) the predictive performance of modelling based on regression trees decreased slightly. We do not report more detailed results here, as there were only minor differences compared to the results reported in Fritsch et al. (2019).

Discussion

Note that in this application, the true distribution of observed and missing values, as well as the mechanism generating the missing values is unknown. This is typical for empirical datasets, where applied researchers are generally unaware of the ‘true’ impact of missing values on the obtained inferences. We compared estimates based on the complete cases – under the assumption of MCAR with estimates obtained by QR-based MI under the assumption of MAR. Although we cannot rule out MNAR mechanisms, our analyses provide empirical evidence that MAR processes were negligible and did not substantially affect the modelling results. Given that the MAR assumption holds – at least approximately – and that both, the substantive model and the imputation model are correctly specified (an untestable assumption in the analysis of empirical datasets), the reported MI inferences can be considered as unbiased.