Authors: Shahab Jolani (https://ror.org/02jz4aj89Care and Public Health Research Institute (CAPHRI), Maastricht University, Maastricht, The Netherlands)
Categories: Research Article, individual participants data (IPD) meta-analysis, missing data, multilevel model, multiple imputation (MI), multivariate imputation using chained equations (MICE)
Source: Research Synthesis Methods
Authors: Shahab Jolani
Modern quantitative evidence synthesis methods often combine patient-level data from different sources, known as individual participants data (IPD) sets. A specific challenge in meta-analysis of IPD sets is the presence of systematically missing data, when certain variables are not measured in some studies, and sporadically missing data, when measurements of certain variables are incomplete across different studies. Multiple imputation (MI) is among the better approaches to deal with missing data. However, MI of hierarchical data, such as IPD meta-analysis, requires advanced imputation routines that preserve the hierarchical data structure and accommodate the presence of both systematically and sporadically missing data. We have recently developed a new class of hierarchical imputation methods within the MICE framework tailored for continuous variables. This article discusses the extensions of this methodology to categorical variables, accommodating the simultaneous presence of systematically and sporadically missing data in nested designs with arbitrary missing data patterns. To address the challenge of the categorical nature of the data, we propose an accept–reject algorithm during the imputation process. Following theoretical discussions, we evaluate the performance of the new methodology through simulation studies and demonstrate its application using an IPD set from patients with kidney disease.
Systematically and sporadically missing data are very common in individual participant data meta-analysis (IPD-MA).Multiple imputation of missing data in IPD-MA should account for clustering and heterogeneity between studies.Multilevel imputation methods are available for continuous variables with both systematically and sporadically missing data.
We develop a new class of imputation methods for categorical variables in IPD-MA while addressing the challenge of simultaneous presence of systematically and sporadically missing data.
Imputation of missing values in IPD-MA allows using full data potential which can provide additional insight for evidence synthesis.The proposed imputation methodology can directly be applied to other settings of clustered (or multilevel) data such as longitudinal studies.
Quantitative evidence synthesis is an important aspect of contemporary clinical research and medical decision making. It is widely used to summarize the effectiveness of medical interventions, ^1^ the accuracy of diagnostic tests, ^2^ the association of prognostic factors, ^3^ or even the performance of published prediction models. ^4^ Historically, evidence synthesis originates from the desire to summarize the results from multiple related studies, often through meta-analysis techniques that aggregate published data (e.g., estimates of relative treatment effect) into a weighted average. However, there is a growing trend toward combining patient-level data from multiple studies and performing a so-called individual participant data meta-analysis (IPD-MA). ^5^
^,^
^6^
IPD-MA represents a new era in evidence synthesis, offering several advantages over traditional meta-analyses that are solely based on aggregate data. Most notably, the use of IPD allows for more precise tailoring of diagnostic strategies, identification of risk and prognostic factors, and personalizing risk prediction. ^7^ Moreover, pooling of IPD is particularly sui3 for identifying modifiers of relative treatment effect ^8^ and developing models to predict absolute outcome risk or absolute treatment effects in individual patients. ^9^
A key challenge in combining multiple sources of IPD is the presence of between-study heterogeneity. Briefly, this implies that studies differ in aspects such as outcome occurrence (baseline risk) or the magnitude of predictor-outcome associations. Such heterogeneity can arise from differences in participant eligibility criteria, variable and outcome definitions, measurement methods, or treatment protocols across studies. The presence of between-study heterogeneity is an important concern in any meta-analysis, as it may substantially affect the interpretation of pooled summary estimates. ^10^
Another common issue in meta-analysis of IPD is the presence of missing data. This situation typically occurs when variables of interest have not been measured in one or more studies (resulting in completely or systematically missing) or in some participants within studies (resulting in partially or sporadically missing). For example, in the GREAT Network study, ^11^ the biomarker brain natriuretic peptide (BNP) was used to explain left ventricular ejection fraction. However, BNP is a relatively recent technique which was unavailable (thus, systematically missing), in several studies within the IPD-MA of the GREAT study. While it is possible to restrict the analysis to participants with complete data, the potential loss of efficiency and validity is generally undesirable. ^12^ Therefore, it is generally recommended to adopt multiple imputation (MI) ^13^ methods when one or more studies are affected by missing data. Briefly, MI generates multiple versions of the original dataset(s) by replacing missing values with imputations that are based on observed data. A common approach is to use chained equations to impute each variable sequentially, conditional on the other variables. ^14^
^,^
^15^ This method allows for a great deal of flexibility as each variable (e.g., binary, ordinal, continuous, etc.) can be imputed using an appropriate functional form, such as logistic or probit models for binary variables (see van Buuren ^16^ ). Although MI can be relatively straightforward to implement, it requires careful consideration when participants are clustered within different studies or centers. Imputation models should account for variability within and across studies; otherwise, the imputed values and their uncertainty may no longer be valid. ^17^
^,^
^18^
For clustered data such as IPD, the multivariate imputation using chained equations (MICE) algorithm naturally extends to multilevel settings. This requires specifying conditional imputation models that incorporate random effects, known as multilevel or hierarchical imputation models. Although several multilevel imputation methods (MLMIs) have been proposed (for a recent overview, see Audigier et al. ^19^ ) they have at least two (i) they are primarily developed for continuous (and normally distributed) variables and (ii) they are not designed to address the combined presence of sporadically and systematically missing data.
Several multilevel imputation routines are available for continuous variables within the MICE framework. Van Buuren ^20^ and Yucel et al., ^21^ among others, developed multilevel imputation approaches based on linear mixed-effects models to handle sporadically missing data in continuous variables. While these methods effectively account for the multilevel structure of the data, they are not designed to handle systematically missing data. Further, Resche-Rigon and White ^22^ and Jolani ^23^ developed multilevel imputation approaches that address both systematically and sporadically missing data, but these are limited to continuous variables. In addition, Resche-Rigon and White ^24^ proposed a two-stage multilevel imputation approach for systematically and sporadically missing data in continuous variables.
For categorical data, MLMIs are scarce and primarily limited to binary variables. Within the MICE framework, Yucel et al. ^21^ developed a MLMIs for binary variables using the generalized linear mixed-effects models. However, this approach only addresses sporadically missing data and is not suitable for systematically missing data. Conversely, Jolani et al. ^25^ proposed a general class of MLMIs for categorical variables with systematically missing data, but their method does not accommodate sporadically missing data. Motivated by the one-stage approach of Jolani et al. ^25^ and two-stage approach of Resche-Rigon and White, ^24^ Audigier et al. ^19^ suggested treating sporadically missing values in binary variables as systematically missing and imputing them accordingly using the one-stage or two-stage multilevel approaches, respectively.
The combined presence of systematically and sporadically missing data in categorical variables is a common challenge in IPD-MA. To date, the MICE framework still lacks principled MLMIs that address both systematically and sporadically missing data in categorical variables. To bridge this gap, we propose a new class of multilevel imputation methodologies for IPD-MA and beyond, specifically designed to address the challenges posed by both systematically and sporadically missing data in categorical variables.
This manuscript is organized as follows. Section 2 introduces an innovative multilevel imputation methodology for categorical variables by utilizing an accept–reject sampling procedure. We mainly focus on the development of MLMIs for a broad class of families, including models for binary and count variables. Section 3 evaluates the performance of the proposed methodologies through extensive simulation studies across a wide range of scenarios, including varying degrees of between-study heterogeneity. Section 4 demonstrates the implementation of the new methodologies in an IPD-MA of chronic kidney disease (CKD) patients. Section 5 concludes with a discussion and final remarks.
Suppose we aim to perform an IPD-MA from n studies. Let
denote a value of the random variable Y for subject
in study
. We assume the following generalized linear mixed-effects model for the IPD (2.1)where
denotes a link function,
is a
vector of potential predictors of outcome
,
is typically a subset of vector
,
is a vector of fixed-effects parameters, and
is a
vector of random effects following a multivariate normal distribution,
with a
covariance matrix
. We denote the diagonal elements of
with
and its off-diagonal elements with
,
and
. We refer to such variances and covariances as random-effects parameters, and
. Furthermore, we assume that the vectors
for all
are independent and identically distributed and that the vectors
are their realizations.
The analysis model 2.1 encompasses a general family of distributions. For continuous outcome
, for instance, the link function is identity (
) so we use For binary outcome
, as another example, the link function is logit (
) and the analysis model 2.1 is represented by
MICE is blind to the role of variables in the analysis, and each variable is imputed in turn conditional on the other variables. Let
denote a full data matrix containing a set of explanatory variables
and outcome Y, and
. Without loss of generality, we assume that each column of
has some missing entries. The MICE (or fully conditional specification) approach specifies an appropriate conditional regression model for each variable (i.e., each column of
) given the other variables, that is, where
denotes the unknown parameters of each conditional distribution. These conditional models are the basis for approximating the posterior predictive distribution of the missing data. To account for the hierarchical data structure (i.e., individuals within studies), we apply a linear mixed-effects model for continuous variables, a logistic mixed-effects model for binary variables, a Poisson mixed-effects model for count variables, and so on. The set of specified conditional multilevel models forms a cycle, a few of which typically should be repeated successively to achieve an adequate approximation to the marginal posterior predictive distribution of the missing data. Specifically, missing values of
are imputed conditional on the other variables. Subsequently, missing values of
are imputed using the recent imputations of
and other variables. The imputed values from the last cycle are eventually considered as one set of imputations. Replicating the whole process multiple times produces MI sets.
We throughout assume the missing data mechanism is missing at random (MAR) ^26^ implying that the probability that any data value is missing may depend on quantities that are observed but not quantities that are missing. Therefore, the imputation model that we focus on would be appropriate under missing completely at random (MCAR) or MAR assumption.
Because MICE is implemented on a variable-by-variable basis, we focus on the imputation of a single (categorical) variable in this section. For subject j in study i, suppose
represents the Kth (last) incomplete variable in the data matrix
. We define the following multilevel imputation model for
, which is a generalized linear mixed-effects model conditional on the remaining
variables (2.2)where g is the link function and
. For brevity, we assume a random intercept and a random slope model for
only in the above imputation model, so
for
, where
is a
covariance matrix containing the random-effects parameters. The fixed-effects parameters are denoted by
in the imputation model 2.2. Following Rubin ^13^ , the formal procedure to obtain imputations for
consists of the following Estimating the parameters
altogether with the random effects
in model 2.2 using the observed data.Drawing
and subsequently
from their observed-data posterior distributions.Imputing missing values of
from the conditional predictive distribution
Although the implementation of step (1) is rather straightforward, drawing the parameters from the posterior distributions in step (2) is cumbersome as the unconditional distributions for these parameters cannot generally be obtained in closed form. Thus, Markov chain Monte Carlo methods are typically employed by combining steps (1) and (2) to estimate and obtain random draws of
and
(see, among others, Drechsler ^27^ ). Such iterative algorithms are, however, unattractive within the MICE framework because a Gibbs sampler needs to be iterated within each conditional model of a cycle. Jolani et al. ^25^ proposed a simplification over the full Gibbs sampler that relies on the conditional independence between
and
and requires no iteration. In short, inference about
can be separated into two conditionally independent parts assuming that the random effects
are known. Here, standard routines for the (generalized) linear mixed-effects models (e.g., the glmer function in the R package
) are used to estimate the parameters in step (1). Random draws of
and
are then obtained independently conditional of the estimated random effects
from step (1) (see Jolani ^23^ for details).
This imputation methodology was originally developed by Jolani et al. ^25^ for any systematically missing variable where the random effects
are drawn from the unconditional (prior) distribution
. Jolani ^23^ further extended the proposed imputation methodology to systematically and sporadically missing continuous variables. Extensions to categorical variables with sporadically missing data are an uneasy task because the conditional posterior distribution of random effects does not have a standard form as opposed to continuous variables. To elaborate on this point, consider an incomplete variable Z. For study i with sporadically missing values, we need to make a draw from the posterior distribution
, which can be approximated by (2.3)From model 2.2, we know that
. If Z is normally distributed (and thus continuous), it follows that
is normally distributed too, and so a random draw of
can easily be obtained conditional on
. However, it is difficult to simulate directly from this distribution if Z is categorical (e.g., binary). Therefore, we propose an accept–reject sampling method, ^28^ to draw a random value
from the posterior distribution
.
The accept–reject sampling method involves obtaining draws from a proposal density (which is easier to sample from) until a draw satisfies a particular condition. We choose
as a proposal density, from which samples can easily be drawn. The method then requires the ratio of the target density (i.e.,
) to the proposal density be bounded above a constant quantity M. It is easy to show that this ratio is proportional to Following Robert and Casella, ^28^ the bound M can be taken to be the likelihood function in equation (2.3) evaluated at the maximum likelihood estimates. The algorithm is then completed when we sample
from the proposal density and U from the uniform distribution on (
). The drawn sample
is accepted if
. Otherwise, a new pair (
) is drawn. After completing the above steps (1) and (2), imputations are obtained in step (3) using an appropriate generalized mixed-effects model.
The proposed accept–reject method is general and can be applied to any categorical variables as long as these include the family of generalized linear mixed-effects models. Also, it should be emphasized that, for a given incomplete variable, the proposed accept–reject sampling method is required for studies with sporadically missing data. For studies with systematically missing data, the random effects
are drawn from the unconditional distribution
(see Jolani et al. ^25^ ). As a showcase, we provide computational details of the proposed imputation methodology for a count variable with systematically and sporadically missing data.Example 2.1.Multilevel imputation of a count variable
For simplicity suppose the data matrix
contains three variables
, and
. Further, assume that the count variable
is sporadically and systematically missing in m and
studies and
and
are the observed and missing part of
respectively. For
, we define the following Poisson mixed-effects model as the imputation where the link function is the natural logarithm, and the model parameters are defined as in model 2.2. Assuming the parameters are a priori independent (i.e.,
) and specifying the standard prior distributions
and
, the imputation procedure consists of the following Obtain the restricted maximum likelihood estimates
and
using
.Obtain the random effects
, where
and calculate
.Obtain a random draw
.Obtain a random draw
where W is a Wishart distribution with m degrees of freedom and a
scale matrix parameter
. (here
For each study i,
If
is sporadically missing, draw
from the developed accept–reject sampling algorithm.If
is systematically missing, draw
from
.Impute
from the Poisson mixed-effects model It should be noted that the scale matrix parameter
is calculated from studies for which
is sporadically missing and studies with systematically missing do not contribute, and consequently this may underestimate
. Further, the above procedure can easily be modified for imputation of other types of categorical variables.
A set of simulation studies was considered to assess the performance of the proposed multilevel imputation methodology by varying the between-study heterogeneity, the size of studies, the proportion of systematically missing data, the missing data mechanism, and the type of (incomplete) variables.
We began by generating a triple
where
is a binary variable,
is a count variable, and Y is the outcome of interest (continuous or binary). Covariates
and
were generated from a joint mixed-effects model where
and
are the fixed-effects parameters and the random-effects
with The continuous outcome variable Y was subsequently generated from the linear mixed-effects model (3.1)where
represents the regression coefficients of interest (i.e., the fixed-effects parameters), the residual errors
, and the random effects
with representing the random-effects parameters. For the binary outcome Y, we similarly defined the generalized linear mixed-effects model (3.2)with the same assumption for
.
The parameter values were chosen to mimic the structure of an IPD-MA investigating potential factors associated with smoking in patients with CKD. Further details on the data are provided in the Application section. We chose alcohol consumption status and the number of comorbidities (e.g., hypertension, obesity, etc.) to motivate the distribution of the binary and count variables
and
, respectively. For the outcome variable y, smoking status was used to motivate the parameters of the analysis model 3.2 (and 3.1). We considered three levels of between-study heterogeneity—weak, moderate, and strong—by varying elements of the covariance matrix
. Details about the parameter values used for the simulation study are provided in the Appendix.
For each simulated dataset, we introduced systematically and sporadically missing values in variables
and
. These variables were independently systematically missing with probability
under the MCAR assumption, which is more plausible in IPD settings. We considered two proportions of systematically missing
and
, resulting in approximately
and
joint systematically missing data in both
and
. In addition, for studies where a variable was not systematically missing, sporadically missing values were created independently for each variable under the MCAR and MAR assumptions. For the MAR scenarios, the probability that a data value is missing depended on the outcome variable (which was fully observed). Specifically, we used the following models to create sporadically missing data in
and
, where
(
) if
(
) is observed and
(
) if
(
) is sporadically missing. Under MCAR,
, and we fixed the proportion of sporadically missing data to
in each of
and
. Under MAR, we set
and
and adjusted the intercept values accordingly. For example, setting
resulted in about
sporadically missing data in
when the outcome variable y was continuous.
Finally, we considered five study sizes
with a fixed total sample size of 3,000 participants, resulting in 500 subjects
, 300 subjects
, 200 subjects
, 150 subjects
, and 100 subjects
per study. This led to a total of 120 simulation scenarios (2 types of outcomes
3 levels of between-study heterogeneity
5 study sizes
2 proportions of systematically missing data
2 types of missing data mechanism), each of which was replicated 1000 times.
The simulation study evaluates five the newly developed MLMI, the two-stage MLMI (2STG), ^24^ the stratified imputation method (STI), complete case analysis (CCA), and the analysis of original data before the introduction of missing values (REF). The last method (REF) serves as a reference to assess the performance of the methods.
The STI method includes a study indicator as a categorical variable to the imputation model to account for the presence of clustering and uses standard imputation routines for non-clustered data, such as predictive mean matching, to generate imputations. However, for studies with systematically missing data, the intercept cannot be unidentified in the STI imputation method. To address this, the method uses the average of the observed study-intercept terms to obtain imputations for studies with systematically missing data. The CCA method, on the other hand, excludes subjects or studies with missing values from the analysis.
For all imputation methods (MLMI, 2STG, and STI), the conditional imputation models include all available covariates. In the MLMI and 2STG methods, the imputation models additionally include a random intercept and a random slope for the outcome variable y. Each incomplete dataset was imputed five times, with 10 cycles of the MICE algorithm used to obtain each imputed dataset. Afterward, each imputed dataset was fitted to the analysis model 3.1 or 3.2, and the results were pooled using Rubin’s rule. All analyses were conducted using the R packages
and
.
The primary parameters of interest are the fixed-effects parameters
and the random-effects parameters
and
in model 3.1 or 3.2. For each method, we calculated the bias, root mean squared error (RMSE), model-based standard error (Model SE), empirical Monte Carlo standard error (Emp SE), and the coverage rate of the
confidence interval (CR).
This section presents the simulation study results for the binary outcome. The findings for the continuous outcome are included in the Supplementary Material.
Tables 1–3 show the results for the fixed-effects and random-effects parameters with
studies under MAR across the five different methods under weak, moderate, and strong between-study heterogeneity, respectively. Overall, all methods provided satisfactory estimates for the fixed-effects parameters
except for CCA where the bias for the estimate of
was moderate. Notably, both STI and 2STG exhibited a slight bias in the estimates of
and
under strong between-study heterogeneity. Moreover, the CCA method had the largest RMSE, likely due to the proportion of missing data as it rose with the rate of systematically missing data. The RMSEs of the other methods were comparable.Table 1Estimates of the fixed- and random-effects parameters in the simulation study for the binary outcome with n = 10 studies and weak between-study heterogeneitySystematically missing10%30%REFCCASTIMLMI2STGREFCCASTIMLMI2STG
=
1.85 Estimate
1.841
2.119
1.849
1.850
1.844
1.843
2.127
1.859
1.855
1.848 Bias0.009
0.2690.0010.0000.0060.007
0.277
0.009
0.0050.002 Model.SE0.1350.1770.1670.1670.1660.1370.2400.2490.2270.203 Emp.SE0.1410.1780.1630.1620.1600.1400.2460.2050.1900.181 CR%91.466.893.394.294.391.377.496.695.895.7 RMSE0.1410.3220.1630.1620.1600.1400.3700.2050.1900.181
= 1.05 Estimate1.0381.0441.0521.0501.0381.0491.0531.0781.0591.052 Bias
0.012
0.0060.0020.000
0.012
0.0010.0030.0280.0090.002 Model.SE0.1840.2290.1730.2050.2070.1830.2910.1870.2410.236 Emp.SE0.1900.2300.2130.2130.2140.1910.3050.2430.2490.241 CR%91.992.988.490.991.091.990.284.390.391.1 RMSE0.1900.2300.2130.2130.2140.1910.3040.2450.2490.241
=
0.04 Estimate
0.042
0.043
0.040
0.042
0.039
0.041
0.045
0.045
0.041
0.039 Bias
0.002
0.0030.000
0.0020.001
0.001
0.005
0.005
0.0010.001 Model.SE0.0430.0610.0800.0530.0530.0430.0810.1700.0670.067 Emp.SE0.0420.0610.0570.0500.0510.0430.0780.0930.0620.062 CR%95.795.098.095.394.695.396.299.094.494.9 RMSE0.0420.0610.0570.0500.0510.0430.0780.0930.0620.062
= 0.316 Estimate0.2400.2280.2660.3010.3020.2540.2420.3230.3870.372 Bias
0.076
0.088
0.050
0.015
0.014
0.062
0.0740.0060.0710.056 RMSE0.1760.2020.1550.1500.1450.1670.2290.1460.2340.157
= 0.500 Estimate0.4550.4490.3380.4260.4370.4500.4060.3000.4220.423 Bias
0.045
0.051
0.162
0.074
0.063
0.050
0.094
0.200
0.078
0.077 RMSE0.1620.2030.1980.1880.1800.1570.2540.2280.2210.205 Note: REF indicates the results that were obtained before missing data were introduced and can be viewed as a benchmark for comparing the performance of methods that are applied after missingness is complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). The following values are mean of estimates (Estimate), bias (Bias), mean of standard error (Model SE), empirical standard error (Emp SE), the coverage rate of 95% confidence interval (CR), and the root of mean squared error (RMSE). Table 2Estimates of the fixed- and random-effects parameters in the simulation study for the binary outcome with n = 10 studies and moderate between-study heterogeneitySystematically missing10%30%REFCCASTIMLMI2STGREFCCASTIMLMI2STG
=
1.85 Estimate
1.844
2.124
1.824
1.846
1.829
1.851
2.136
1.850
1.860
1.843 Bias0.006
0.2740.0260.0040.021
0.001
0.2860.000
0.0100.007 Model.SE0.2810.3280.2840.3070.2990.2800.4260.3350.3530.335 Emp.SE0.3020.3560.3200.3220.3160.2970.4530.3490.3480.334 CR%90.581.789.791.591.190.783.991.493.093.2 RMSE0.3020.4490.3210.3220.3170.2970.5360.3480.3480.334
= 1.05 Estimate1.0401.0431.0331.0431.0191.0571.0761.0751.0681.045 Bias
0.010
0.007
0.017
0.007
0.0310.0070.0260.0250.018
0.005 Model.SE0.2990.3490.2410.3290.3160.2990.4390.2350.3870.357 Emp.SE0.3110.3760.3370.3480.3360.3150.4730.3780.3910.377 CR%89.689.881.288.988.691.888.675.589.689.0 RMSE0.3110.3750.3380.3480.3370.3150.4730.3780.3910.377
=
0.04 Estimate
0.038
0.039
0.038
0.039
0.036
0.040
0.040
0.039
0.038
0.036 Bias0.0020.0010.0020.0010.0040.0000.0000.0010.0020.004 Model.SE0.0440.0630.0820.0550.0560.0450.0830.1700.0690.071 Emp.SE0.0440.0640.0600.0540.0540.0450.0820.0920.0640.063 CR%95.894.897.794.595.193.994.499.393.595.4 RMSE0.0440.0640.0600.0540.0540.0450.0820.0920.0640.063
= 0.866 Estimate0.7850.7590.7360.8040.7900.7800.7100.7470.8490.825 Bias
0.081
0.107
0.130
0.062
0.076
0.086
0.156
0.119
0.017
0.041 RMSE0.2410.3120.2500.2320.2300.2480.4070.2440.2350.221
= 0.922 Estimate0.8350.8100.5830.8040.7900.8360.7390.4910.7900.753 Bias
0.087
0.112
0.339
0.118
0.132
0.086
0.183
0.431
0.132
0.169 RMSE0.2610.3250.3880.3070.2990.2510.4110.4680.3590.337 Note: REF indicates the results that were obtained before missing data were introduced and can be viewed as a benchmark for comparing the performance of methods that are applied after missingness is complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). The following values are mean of estimates (Estimate), bias (Bias), mean of standard error (Model SE), empirical standard error (Emp SE), the coverage rate of 95% confidence interval (CR), and the root of mean squared error (RMSE). Table 3Estimates of the fixed- and random-effects parameters in the simulation study for the binary outcome with n = 10 studies and strong between-study heterogeneitySystematically missing10 %30 %REFCCASTIMLMI2STGREFCCASTIMLMI2STG
=
1.85 Estimate
1.858
2.125
1.714
1.852
1.754
1.869
2.140
1.708
1.875
1.746 Bias
0.008
0.2750.136
0.0020.096
0.019
0.2900.142
0.0250.104 Model.SE0.5240.5920.4530.5480.5020.5270.7870.4720.6270.543 Emp.SE0.5400.6160.5260.5630.5110.5510.8430.5870.6340.569 CR%90.987.786.291.991.891.985.585.193.492.0 RMSE0.5400.6750.5430.5620.5200.5510.8910.6030.6340.578
= 1.05 Estimate1.0581.0640.9421.0580.9471.0501.0710.9171.0450.904 Bias0.0080.014
0.1080.008
0.1030.0000.021
0.133
0.005
0.146 Model.SE0.5330.5980.3760.5730.5160.5370.8090.3520.7200.597 Emp.SE0.5500.6420.5620.5940.5490.5600.8520.6420.7020.634 CR%91.489.478.691.889.991.186.868.989.886.7 RMSE0.5500.6420.5720.5940.5590.5590.8520.6560.7020.650
=
0.04 Estimate
0.037
0.039
0.035
0.037
0.034
0.042
0.046
0.042
0.039
0.038 Bias0.0030.0010.0050.0030.006
0.002
0.006
0.0020.0010.002 Model.SE0.0460.0640.0830.0570.0600.0460.0860.1680.0720.077 Emp.SE0.0440.0630.0580.0520.0530.0450.0830.0930.0630.066 CR%95.896.398.397.596.195.895.498.994.995.9 RMSE0.0440.0630.0580.0520.0540.0450.0830.0930.0630.066
= 1.658 Estimate1.5331.4961.2851.5061.4041.5431.4261.2451.5411.413 Bias
0.125
0.162
0.373
0.152
0.254
0.115
0.232
0.414
0.117
0.245 RMSE0.4420.5200.5070.4360.4420.4320.7190.5230.4380.429
= 1.688 Estimate1.5551.5041.0271.4861.3701.5671.4630.8951.5431.350 Bias
0.133
0.184
0.661
0.202
0.319
0.121
0.225
0.794
0.145
0.338 RMSE0.4420.5140.7120.4800.5170.4370.7370.8380.5870.587 Note: REF indicates the results that were obtained before missing data were introduced and can be viewed as a benchmark for comparing the performance of methods that are applied after missingness is complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). The following values are mean of estimates (Estimate), bias (Bias), mean of standard error (Model SE), empirical standard error (Emp SE), the coverage rate of 95% confidence interval (CR), and the root of mean squared error (RMSE).
The performance of methods varied with respect to confidence interval coverage. For the coefficient
of the binary covariate, all methods were close to the nominal level since the variance of
was generally underestimated, even for the reference method (REF). Nevertheless, the STI method had the lowest coverage rate, dropping below
under strong between-study heterogeneity. Interestingly, STI resulted in coverage rates that exceeded the nominal level for the coefficient
. These are likely due to improper modeling of the multilevel data structure during imputation. The other methods achieved the nominal
level for the coefficient
of the count variable. Despite not underestimating the variance of
, the CCA method led to severe undercoverage for
.
For the random-effects parameters
and
, all methods, including the reference, yielded downwardly biased estimates, likely due to the shrinkage effect in random-effects models and the limited number of studies, and the bias increased with the degree of between-study heterogeneity. For
, the largest bias was observed with CCA under weak and moderate between-study heterogeneity, whereas it occurred with STI under strong between-study heterogeneity. For
, the largest bias was consistently observed with STI across all between-study heterogeneity scenarios. In contrast, the MLMI and 2STG methods produced estimates close to the reference method, with 2STG exhibiting a slightly larger bias under strong between-study heterogeneity. When comparing RMSE, CCA, and STI had the highest values for
, and
, respectively, followed by the MLMIs, where RMSE increased as the proportion of systematically missing data grew.
To evaluate the performance of the methods across different study sizes (
), Figures 1–3 present the results for the fixed-effects parameters when the rate of systematically missing data was
, and the missing data mechanism was MAR (for sporadically missing data). In Figures 4 and 5, the focus is on the random-effects parameters
and
.Figure 1Bias of the fixed-effects estimates with 10% systematically missingness. Methods include reference (REF-before introducing missing data), complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). Figure 2Coverage rate of the 95% confidence interval for the fixed-effects parameters with 10% systematically missingness. Methods include reference (REF-before introducing missing data), complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). Figure 3Root mean squared error (RMSE) of the fixed-effects estimates with 10% systematically missingness. Methods include reference (REF-before introducing missing data), complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). Figure 4Bias of the random-effects estimates with 10% systematically missingness. Methods include reference (REF-before introducing missing data), complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG). Figure 5Root mean squared error (RMSE) of the random-effects estimates with 10% systematically missingness. Methods include reference (REF-before introducing missing data), complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG).
Figure 1 illustrates the bias in estimates of
and
for the analysis model 3.2. Under weak and moderate between-study heterogeneity, all methods resulted in negligible bias for the coefficient
of the binary covariate, except for 2STG. Surprisingly, the 2STG method tended to underestimate
, with the bias gradually increasing as study size grew. Under strong between-study heterogeneity, both the 2STG and STI methods tended to underestimate
across different study sizes, while the MLMI method showed only a negligible bias in this scenario. For the coefficient
of the count covariate, on the other hand, bias remained trivial across all methods and study sizes.
Figure 2 presents the coverage rates of the
confidence intervals for the fixed-effects parameters
and
. STI exhibited significant coverage issues for the coefficient
of the binary covariate, especially under moderate and strong between-study heterogeneity, where coverage dropped below
. Additionally, STI led to overcoverage for the coefficient
of the count covariate due to inflated standard errors across different study sizes and between-study heterogeneity scenarios. Conversely, the other methods performed reasonably well, maintaining coverage rates close to the nominal
level (except for 2STG, which showed slightly lower coverage rates under strong between-study heterogeneity).
Figure 3 displays the RMSE of
and
across different methods and study sizes. As expected, CCA showed the highest error, particularly for the coefficient
of the binary covariate with strong between-study heterogeneity. STI, 2STG, and MLMI had comparable RMSE values, though marginally higher than REF. We also noticed that RMSE of
increased with the degree of between-study heterogeneity across all methods.
Figure 4 illustrates the bias in estimating the random-effects parameters
and
. As discussed earlier, all methods exhibited marginal underestimation of random-effects parameters. Most notably, STI suffered from a large bias, particularly when estimating the random-effect parameter
or when the between-study heterogeneity was strong. The multilevel imputation methods (MIML and 2STG) produced comparable estimates of the random-effects parameters and were generally in line with the reference method under weak and moderate between-study heterogeneity. However, under strong between-study heterogeneity, 2STG exhibited a noticeable bias, followed by MLMI, though the bias from MLMI remained trivial. Finally, CCA showed the least bias across different scenarios, although it increased with smaller study sizes.
Figure 5 shows the RMSE for the estimates of random-effects parameters
and
. For all methods, RMSE generally decreased as the number of studies increased but increased with greater between-study heterogeneity. CCA and STI consistently exhibited higher RMSE for
and
, respectively. In contrast, MLMI and 2STG tended to have lower RMSE, closely aligning with the reference method.
The findings for the binary outcome with
systematically missing data are presented in Figures A.1–A.5 in the Appendix. Although these results generally align with those in Figures 1–5, differences in model performance are more pronounced in Figures A.1–A.5 due to the higher rate of systematically missing data. Most notably, STI exhibited a substantial decline in coverage rates for the coefficient
of the binary covariate, while it increased sharply for the coefficients
of the count covariate. Additionally, 2STG showed a gradual decline in coverage rates for
under strong between-study heterogeneity (Figure A2). Lastly, CCA showed a modest increase in RMSE of
(Figure A3).
It is also worth noting that the performance of methods changed slightly when the number of studies was small (
) and the rate of systematically missing data was
. Specifically, CCA exhibited a rapid decline in coverage rate for
, dropping below
, whereas the decline was less pronounced for the other methods, including REF (Figure A2). Moreover, CCA produced downward-biased estimates for the random-effects parameters
and
, while MLMI exhibited an upward bias only for
(Figure A4). This may be due to the limited number of studies available, which could affect imputation precision for MLMI and reduce analysis accuracy for CCA.
Finally, we did not report the results for MCAR cases (for sporadically missing data), as the methods showed no substantial performance differences across scenarios.
The results for the continuous outcome are presented in the Supplementary Material. Specifically, Tables S1–S3 in the Supplementary Material provide the results for the fixed-effects and random-effects parameters with
studies. Figures S1–S10 in the Supplementary Material illustrate the findings for
and
systematically missing data, respectively.
Overall, the performance of the methods remained consistent with that observed in the binary case. However, CCA performed poorly. Specifically, it exhibited the largest bias for the coefficient
of the binary covariate across most scenarios, while 2STG had the second largest bias, but only when study size was large and between-study heterogeneity was strong (Figures S1 and S6 in the Supplementary Material). The bias with the other methods was fairly trivial. Furthermore, both CCA and STI showed coverage issues for
, whereas the other methods achieved improved coverage rates (Figures S2 and S7 in the Supplementary Material). For the coefficient
of the count covariate, all methods performed well in terms of bias and coverage rates, similar to the binary case.
For the random-effects parameters
and
, CCA and STI showed biased estimates with CCA showing the largest bias for
and STI for
(Figures S4 and S9 in the Supplementary Material). MLMI and 2STG performed comparably, exhibiting less bias across different scenarios. Nevertheless, 2STG tended to show greater bias for
as study size and between-study heterogeneity increased.
As in the binary case, the methods showed no substantial performance differences across scenarios under the MCAR mechanism, so the results are not presented.
In summary, the findings indicate that naïve methods like CCA and STI have major drawbacks. Specifically, CCA struggles with the reduced number of studies, leading to increased errors and potentially biased estimates, and lower coverage rates, particularly when the proportion of systematically missing data increases. Similarly, the STI method showed bias and coverage issues in many scenarios due to inadequate modeling of the multilevel structure. In contrast, the multilevel imputation methods (MLMI and 2STG), which appropriately account for the multilevel structure, offer more reliable inferences, with smaller bias and confidence intervals closer to the nominal level. Nonetheless, 2STG exhibited slight underestimation in some scenarios, possibly due to treating categorical variables with sporadically missing values as systematically missing and imputing them accordingly. When the study size was small (
), the performance of all methods generally declined, highlighting their limitations in settings with limited data.
To demonstrate the approaches, we considered an IPD meta-analysis of 15 studies involving 4,774 patients with CKD stages III–V, defined by an estimated glomerular filtration rate below 60. The anonymized individual participant datasets for these studies were retrieved from the PLOS ONE website.
In this purely methodological exercise, we aimed to explore whether alcohol consumption and the number of existing health conditions such as hypertension and diabetes are associated with smoking status in this specific population. Focusing on the analysis model 3.2, we considered smoking status as the outcome y and alcohol consumption (yes vs no) and the number of comorbidities (count) as the covariates
and
, respectively. The number of comorbidities was defined as the sum of the presence of hypertension, diabetes, obesity, and cardiovascular disease. While the analysis model was restricted to these variables, the imputation models also included additional covariates—age and sex—to make the MAR assumption more plausible. ^29^ Table 4 summarizes the percentage of systematically and sporadically missing data for each variable included in this case study.Table 4Percentage of missing data by variable and study in the empirical exampleStudySizeSmokeGenderAgeAlcoholComorbidities18321001025330002603880000041,144000100055510001000657800010007116000100084439005509730005010550001000112700010001240500010001322400000145900010010015450001000
We employed four methods to investigate the strength of association and the degree of between-study heterogeneity for the covariates in model 3.2. First, we performed the CCA, excluding studies with systematically missing data and patients with sporadically missing data. This resulted in a subset of 1,615 patients across six studies.
Second, we applied the stratified MI method (STI), which accounts for heterogeneous intercepts by including study indicators in the imputation model. For STI, binary and count variables were imputed using the logreg and pmm imputation procedures in the R package
, respectively.
Finally, we applied the newly developed MLMI to account for heterogeneity across studies by incorporating joint random effects into the imputation models. Specifically, the imputation models for alcohol consumption (binary) and the number of comorbidities (count) included random effects on the intercept and the outcome variable smoking status, while the other covariates were treated as fixed effects. In addition, the two-stage MLMI (2STG) was included for comparison. For 2STG, the imputation models were identical to those used in MLMI.
For all imputation approaches (STI, 2STG, and MLMI), we allowed 20 iterations of the MICE algorithm to enhance convergence. Following White et al., ^15^ the number of imputed datasets should be at least as large as the percentage of missing data. Since approximately
of the data were missing in this study, we generated 100 imputed datasets from the original data.
Table 5 presents the estimates of the fixed-effects coefficients, their standard errors, and the between-study heterogeneity estimates (i.e., the random-effects parameters
and
) for this case study. Overall, the results align with the findings from the simulation study.Table 5Estimates of the fixed- and random-effects parameters in the empirical exampleFixed effectsRandom effects
EstseEstseEstseEstEstCCA
1.7940.5720.9130.4770.1020.0851.2170.677STI
1.2890.3150.6320.265
0.0140.0451.1290.543MLMI
1.3220.3530.6770.403
0.0120.0471.1370.6162STG
1.1250.4590.3720.600
0.0050.0481.2860.975 Note: Methods include complete case analysis (CCA), stratified multiple imputation (STI), multilevel multiple imputation (MLMI), and two-stage multilevel multiple imputation (2STG).
First, the estimates of the random-effects parameters indicate moderate between-study heterogeneity. The STI method tended to underestimate this heterogeneity, particularly for
, whereas other methods produced comparable estimates. Surprisingly, 2STG overestimated
, with its estimate increasing sharply compared to CCA and MLMI.
Next, CCA produced fixed-effects coefficient estimates that were further from the null compared to those obtained using the MI methods. This suggests that the missing data mechanism might not be completely at random. The MI methods yielded similar estimates for the fixed-effects parameters, but the estimate of the regression slope for alcohol consumption (
) appeared to differ with 2STG compared to the other imputation methods. Although this pattern was also observed in the simulation study (middle panel of Figure A1 in the Appendix), the corresponding coefficient is not statistically significant at the
level, suggesting that differences in point estimates are unlikely to be meaningful.
As expected, CCA overestimated the standard errors of the fixed-effects parameters due to substantial data reduction. In contrast, the STI method underestimated these standard errors relative to the multilevel imputation approaches, likely due to improper handling of between-study heterogeneity.
In summary, the multilevel imputation methods (particularly MLMI), which properly account for between-study heterogeneity, may provide more accurate results than other alternatives.
Meta-analysis of IPD should employ appropriate methodologies to deal with the simultaneous presence of systematically and sporadically missing data. Excluding incomplete studies or cases from the meta-analysis is undesirable, as full data potential is not optimally used, especially when the individual studies are too large or important to be excluded. This can lead to efficiency loss, reduced statistical power due to a smaller sample size, and potentially biased estimates.
Additionally, in the simulation study, CCA yielded conflicting results when comparing binary and continuous outcomes. Specifically, the estimates of the main parameters (
and
) were unbiased for the binary outcome but exhibited bias in the continuous case. This is not surprising because CCA leads to biased regression slope estimates in linear regression when missingness in a covariate depends on the outcome variable Y (i.e., the MAR mechanism). In contrast, in logistic regression with CCA, the regression slope estimates remain unbiased when missingness depends on Y due to the symmetry property of the odds ratio (see Bartlett et al. ^30^ and Carpenter and Kenward ^31^ ).
Hence, imputation strategies become an attractive solution in meta-analysis of IPD. Naïve imputation approaches, such as stratified imputation, which do not fully account for variability within and across studies typically lead to biased parameter estimates, and hence wrong conclusions may be drawn. ^18^ Multilevel (or hierarchical) imputation approaches are therefore a valuable alternative to preserve data integrity in IPD-MA and should be preferred in practice.
Building upon our previous work, we developed new multilevel imputation methodologies for clustered (or multilevel) data, which effectively handle the combined presence of systematically and sporadically missing values in dichotomous or count variables. This methodology utilizes a generalized linear mixed-effects model with random intercept terms and random slopes to complete the imputation task. The simulation and case study results demonstrate that the proposed multilevel imputation approach has desirable statistical properties in terms of bias and coverage rates and maintains appropriate levels of between-study heterogeneity. Furthermore, this methodology facilitates more complex post-imputation analyses by incorporating joint random effects during the imputation phase. Finally, the developed multilevel imputation routines are freely accessible via the popular R package
.
Within the MICE framework, missing values can be imputed per study, a process known as the two-stage imputation approach as opposed to the one-stage imputation approach (such as the proposed MLMI methodology) that imputes missing values within and across studies. The two-stage imputation approach does not impose any restrictions on sporadically missing data. However, it is no longer feasible to impute systematically missing data because the variable is entirely missing in certain studies, so no imputation can be generated for that variable within those studies.
Resche-Rigon and White ^24^ proposed a variant of the two-stage imputation approach (the 2STG method in the simulation study), where a two-stage estimator is used to handle both types of missing data. Compared to our one-stage imputation approach, the two-stage imputation approach ^24^ is computationally faster because estimation is performed in two steps. However, in our simulations, this approach produced mildly biased estimates for categorical covariates, particularly when the between-study heterogeneity was strong. A possible explanation is that it treats categorical variables with sporadically missing data as systematically missing, leading to inaccurate imputations. In addition, this method requires large studies to minimize the small-sample bias of the maximum likelihood estimator and to avoid separability issues, particularly with binary variables. Moreover, it is prone to overfitting when there are many covariates or high rates of missing data within each study. In contrast, the one-stage imputation approach is highly relevant and useful for IPD sets even with relatively few studies, as it can prevent overfitting issues. Audigier et al. ^19^ compared broadly properties of the one-stage and two-stage imputation approaches (together with other imputation methods) through an extensive simulation study.
One limitation of the one-stage imputation approaches within the MICE framework is the compatibility issue inherited from the standard (non-hierarchical) FCS algorithm. The conditional multilevel imputation models may not be consistent with a well-defined joint multilevel model, meaning the imputation models could be generally misspecified. ^24^
^,^
^32^ However, both papers showed via simulation studies that the impact of misspecified imputation models could be minor, particularly when the number of studies is large. Additionally, one-stage imputation approaches, such as MLMI, are computationally intensive and slow due to fitting several mixed-effects models during imputation. This is, however, less of an issue nowadays with fast processors and parallel computing options.
Multilevel imputation models often involve many random effects, which can cause convergence issues, particularly when the number of studies is limited. This can lead to inaccurate imputations, introducing bias in subsequent analyses. As suggested by Jolani et al., ^25^ simplifying the imputation model (e.g., by considering fewer or independent random effects) may help address this issue. Furthermore, in our simulation study, Rubin’s rule was applied to summarize the random-effects parameters (i.e., taking the arithmetic average of the estimated variances of the random effects from the imputed datasets). However, given the skewness of the distribution of these parameters, taking a simple mean might not be the optimal approach. Increasing the number of imputations and using the median of the estimated random parameters could be a better alternative, although further research is needed to systematically evaluate this proposal.
It is also worth mentioning that in IPD-MA, variables with systematically missing data are typically assumed to follow an MCAR or MAR mechanism. This assumption enables leveraging information from studies with observed data to impute plausible values for studies with systematically missing data. However, the possibility of an MNAR mechanism cannot be ruled out. A potential solution is to integrate the proposed multilevel imputation approach with the pattern-mixture approach and conduct a sensitivity analysis.
Although the development of the new imputation methodology was motivated by its usage in the meta-analysis of IPD, this methodology, by no means, is restricted to studies with both systematically and sporadically missing data. The proposed imputation methodology is general and can be directly applied to other hierarchical or multilevel settings, such as longitudinal studies where sporadically missing data are common.
While we focused on logistic and Poisson mixed-effects models in the current article, extending this methodology to other families of distributions, such as gamma or inverse Gaussian, would be relatively straightforward. Additionally, extending the methodology to handle other types of categorical variables (nominal and ordinal) with both systematically and sporadically missing data would be a valuable direction for future research. Another promising line of research is to broaden the methodology to handle treatment-covariate interactions (i.e., effect modifiers) in IPD-MA. Our current methodology does not account for effect modifiers during the imputation process, and addressing this challenge is not trivial. When a covariate has missing observations, its interaction with treatment will also be incomplete, complicating the imputation process.