Multiple imputation for systematically missing effect modifiers in individual participant data meta-analysis

Abstract

Individual participant data (IPD) meta-analysis of randomised trials is a crucial method for detecting and investigating effect modifications in medical research. However, few studies have explored scenarios involving systematically missing data on discrete effect modifiers (EMs) in IPD meta-analyses with a limited number of trials. This simulation study examines the impact of systematic missing values in IPD meta-analysis using a two-stage imputation method. We simulated IPD meta-analyses of randomised trials with multiple studies that had systematically missing data on the EM. A multivariable Weibull survival model was specified to assess beneficial (Hazard Ratio (HR) $=$ 0.8), null (HR $=$ 1.0), and harmful (HR $=$ 1.2) treatment effects for low, medium, and high levels of an EM, respectively. Bias and coverage were evaluated using Monte-Carlo simulations. The absolute bias for common and heterogeneous effect IPD meta-analyses was less than 0.016 and 0.007, respectively, with coverage close to its nominal value across all EM levels. An uncongenial imputation model resulted in larger bias, even when the proportion of studies with systematically missing data on the EM was small. Overall, the proposed two-stage imputation approach provided unbiased estimates with improved precision. The assumptions and limitations of this approach are discussed.

Keywords

Individual participant data meta-analysis systematically missing data two-stage meta-analysis conditional quantile imputation Monte-Carlo simulation treatment effect modification

1. Introduction

Individual participant data (IPD) meta-analysis of randomised trials is a key method to identify and investigate differential treatment effects (effect modification) in medical research.^1–4 Single trials may lack the statistical power to detect detailed subgroup differences in treatment effects.^4–6 For this reason, IPD meta-analyses provide an important opportunity to increase power to detect genuine effect modification.^2,7 A two-stage IPD approach for estimating treatment effect modification mitigates aggregation bias by first evaluating the effect modification within each individual study and then combining the results in the second stage.⁸ However, missing data in individual studies are common and pose challenges when performing an IPD meta-analysis. While sporadically missing data can be managed within a single study, systematically missing data might lead to the exclusion of entire studies.

Systematically missing data is present when one or more variables are not available in one or more studies.^9–11 Variables might be systematically missing for different reasons (e.g., different survey instruments or measurement devices, or lack of information) and often pose practical as well as methodological challenges such as the risk of decreasing the ability to evaluate a broader range of effect modifications. Multiple imputation (MI) is a popular method that can be used to retain studies with systematically missing values in IPD meta-analyses.^10,12–14 In brief, MI ensures that missing data are replaced by their corresponding imputation samples, resulting in $M$ completed datasets where $M$ is the number of imputations. The substantive (outcome) model is then fit in $M$ imputed datasets. Rubin’s rules are applied to pool the $M$ estimates and take into account the variability in the results between the imputed datasets. An important feature of MI is its flexibility, since the imputation model and the substantive model do not have to be the same. However, the two models have to be congenial.¹⁵ In other words, the imputation and substantive model can be derived from a joint distribution by appropriate conditioning. This implies that for the imputation model to be congenial, it must include all variables that are also used in the substantive model (including the outcome and transformations of variables).^16–18

While some studies have looked at systematically missing covariates including confounders,⁹ few have assessed scenarios on when a key effect modifier (EM) might be missing. Therefore, in this analysis we consider the use and evaluation of a two-stage imputation approach to impute systematically missing values of an EM in an IPD meta-analysis with a small number of trials. Identifying, validating, and successfully analysing treatment effect modifications in clinical practice is challenging. Therefore, it is crucial to consider all available studies in an IPD meta-analysis to avoid losing vital information.¹⁹ The rationale for using MI in this context is primarily to retain all studies in the analysis, thereby assessing treatment effect heterogeneity across the entire population.

The remainder of this paper is structured as follows. First, we describe the method employed for imputing systematically missing values. Subsequently, we outline the structure of the simulation study, encompassing the data generating mechanism (DGM), analytical model, and performance criteria.^20,21 The results of the study are then presented, followed by an application of the imputation method to an IPD meta-analysis of 10 randomised controlled trials assessing the efficacy of postoperative radiotherapy (PORT) in patients with completely resected non-small cell lung cancer. We conclude with a discussion of the strengths and limitations of the proposed approach.

2. Conditional quantile imputation

In this section, conditional quantile imputation (CQI) is introduced for discrete systematically missing data in IPD meta-analysis. The approach extends that of the standard approach of imputing categorical data in single studies. An overview of CQI for continuous data can be found elsewhere.^22,23 An explanation of the notations used throughout the paper follows. Let the index $i$ denote the $i$ th individual in each study. Further, let the index $j = 1, 2, \dots, J$ denote the $j$ th study included in the IPD meta-analysis. We consider a discrete variable, $z_{i j}$ , with $k = 1, 2, \dots, K$ levels, and, a $p$ -dimensional set of predictors of $z_{i j}$ , $w_{i j} = (w_{1 i j}, w_{2 i j}, \dots, w_{p i j})$ . Suppose the vector $w_{i j}$ is completely observed for all $i \in 1, \dots, n_{j}$ and all $j \in 1, \dots, J$ , while $z_{i j}$ can have missing values. Let $M_{j} \subset {1, \dots, n_{j}}$ denote the set of indexes corresponding to the individuals with missing values for $z_{i j}$ . The cardinality of $M_{j}$ divided by the sample size $n_{j}$ denotes the proportion of missing data for the $j$ th study, that is $F M_{j} = | M_{j} | / n_{j}$ , ranging from 0 (no missing) to 1 (all missing, i.e., systematically missing). Systematically missing data in a single study is therefore present if $F M_{j} = 1$ . We identify $J = A \cup B$ , with $A = j \in J : F M_{j} = 1$ as the set of studies with systematically missing data and $B = j \in J : F M_{j} < 1$ as the set of studies that can be used to fit the imputation model. The imputation for a discrete systematically missing variable in an IPD meta-analysis can be achieved according to three steps: (1) Specification of the imputation model; (2) prediction of conditional probabilities; (3) imputations of the systematically missing variable.

2.1. The imputation model

An imputation model is specified in studies $J \in B$ . Let $θ_{i k j}$ denote the probability that $z_{i j}$ is equal to level $k$ conditionally on a set of observed variables, the vector of $w_{i j}$ , $θ_{i k j} | w_{i j} = P (z_{i j} = k | w_{i j})$ with $\sum_{k = 1}^{K} θ_{k} = 1$ . A study-specific multinomial logistic regression model is considered for the discrete EM, $z_{i j}$

\ln (\frac{θ_{i k j}}{θ_{i r j}}) = f_{i k j} (w_{i j})

(1)

by specifying a baseline level

r

. Estimates of the regression coefficients,

{\hat{γ}}_{i k j}

, and their variances,

\hat{Var} ({\hat{γ}}_{i k j})

, in the linear predictor

f_{i k j} (w_{i j})

are obtained using the maximum likelihood method separately for each study with complete or partial information.

2.2. Prediction of conditional probabilities

For all studies in $J \in B$ , the estimated regression coefficients ${\hat{γ}}_{i k j}$ are combined using the inverse variance method. Specifically, a weighted average of the imputation regression coefficients from all studies in $J \in B$ is taken, with the weights equal to the inverse of the variance, assigning greater weight to coefficients with smaller variances.²⁴ The combined average regression coefficients are denoted by ${\bar{γ}}_{i k}$ . An estimate of the variance of ${\bar{γ}}_{i k}$ , denoted by $\hat{Var} ({\bar{γ}}_{i k})$ , is estimated using the inverse of the sum of weights

\hat{Var} ({\bar{γ}}_{i k}) = {(\sum_{J \in B} w_{j})}^{- 1}

(2)

where

w_{j} = \hat{Var} ({\hat{γ}}_{i k j})^{- 1}

is the weight assigned to study

j

. The predicted conditional probabilities of falling into the levels of the systematic missing predictor are denoted by

{\hat{θ}}_{i k} = \hat{P} (z_{i} = k | w_{i})

. Next, for each regression coefficient

{\bar{γ}}_{i k}

in studies

J \in B

, including the intercept, and for all linear predictors in the multinomial logistic regression model, a value is drawn from a normal distribution with the mean equal to the average regression coefficients

{\bar{γ}}_{i k}

and standard deviation equal to the square root of the estimated variance of the pooled coefficients.²⁵ In total, for every imputation, there are

(p + 1) \times (K - 1)

independent draws from univariate normal distributions. Specifically, each draw is given by

γ_{i k}^{*} \sim N ({\bar{γ}}_{i k}, \sqrt{\hat{Var} ({\bar{γ}}_{i k})})

(3)

The resulting draws,

γ_{i k}^{*}

, are used to define linear predictors,

{\bar{f}}_{i k} (w_{i})

, which are then utilised to compute conditional predicted probabilities of the missing EM in studies

J \in A

, with

r

as the lowest level

{\hat{θ}}_{i k} = {\begin{cases} \frac{1}{1 + \sum_{k = 2}^{K} e^{{\bar{f}}_{i k} (w_{i})}} & if k = r \\ \frac{e^{{\bar{f}}_{i k} (w_{i})}}{1 + \sum_{k = 2}^{K} e^{{\bar{f}}_{i k} (w_{i})}} & if k > r \end{cases}

(4)

2.3. Imputations

For all studies in $J \in A$ , the inverse of the cumulative distribution function (CDF) – quantile function – is used to assign a value to the missing EM, $z_{i}$ , given the vector of observed predictors $w_{i}$ . The conditional predicted CDF, ${\hat{Θ}}_{i k}$ , is equal to the sum of the probabilities of the EM being less than or equal to $k$ :

{\hat{Θ}}_{i k} = \hat{P} (z_{i} \leq k | w_{i}) = \sum_{k = 1}^{K} {\hat{θ}}_{i k}

(5)

An imputation of the missing values of the EM in any study with systematic missing

z_{i}

is obtained by inverting the CDF of a continuous random uniform distribution

U \sim Uniform (0, 1)

. The imputed value of the EM is obtained by mapping such random draw from a uniform distribution to the predicted cumulative probabilities

{\hat{Θ}}_{i k}

3. Simulation

3.1. Data generating mechanism

We first describe the DGM for the simulation of a single trial followed by a description of the DGM of multiple trials forming the IPD meta-analysis. We then outline the different scenarios considered for the simulation study.

3.1.1. Single trial

We considered a randomised controlled trial ( $n = 500$ ) with a time-to-death outcome, one binary treatment, and one discrete EM. The following random variables defined a single trial:

$x$ is a binary treatment (randomly allocated) generated from a Bernoulli distribution with probability of 0.5.

$z$ is a discrete EM, modifying the effect of $x$ on the outcome at different levels of $z$ . The EM, $z$ , was generated by categorising a random uniform distribution over a (0,1) interval into 3 values (1, 2, and 3) using 0.4, 0.8, and 1.0 as cumulative probabilities.

$t$ is the time from baseline to death (in years) or end of follow-up (10 years), whichever came first.

The outcome random variable,

T

, denoting individual time-to-death conditional on the treatment

x

, covariate

z

(modelled with

3 - 1 = 2

indicator variables) as well as their product terms was described by a multivariable Weibull survival model as follows:

\begin{matrix} P (T > t) = S (t) = e^{- λ t^{γ}} \end{matrix}

(6)

\begin{matrix} t_{i} = {[- \frac{\ln (S_{i})}{λ_{i}}]}^{1 / γ} \end{matrix}

(7)

\begin{matrix} λ_{i} = e^{β_{0} + β_{1} x_{i} + β_{2} I (z_{i} = 2) + β_{3} I (z_{i} = 3) + β_{4} x \cdot I (z_{i} = 2) + β_{5} x_{i} \cdot I (z_{i} = 3)} \end{matrix}

(8)

The shape parameter was defined as

γ = 1.4

signifying an increase of the baseline mortality rate over the follow-up period. Drawing a random value from a continuous uniform distribution for the survival probability,

S

, over a 0 to 1 interval, a random value of the individual time-to-death outcome conditionally on

x

and

z

was generated as follows:

t_{i} = {[- \frac{\ln (S_{i})}{e^{β_{0} + β_{1} x_{i} + β_{2} I (z_{i} = 2) + β_{3} I (z_{i} = 3) + β_{4} x_{i} \cdot I (z_{i} = 2) + β_{5} x_{i} \cdot I (z_{i} = 3)}}]}^{1 / 1.4}

(9)

Any death beyond the follow-up time was set to 10 years and considered censored (i.e., still alive at the end of the follow-up). An indicator variable

d = I (t < 10)

was created to take into account censoring in the estimation of the parameters of the survival model.

3.1.2. Multiple trials

A common and heterogeneous DGM for the generation of multiple trials was considered. We now describe the parameters underlying the Weibull survival model for a common and heterogeneous treatment effect modification.

Common treatment effect modification: The primary characteristic of common treatment effect modification is that the parameter values across the studies do not differ from one another. In other words, all trials in the IPD meta-analysis share a single fixed parameter value for the effect modification. The following values of the parameters were defined:

$β_{0}$ : The mortality rate when all predictors were set to 0 – the intercept – was 10 deaths per 1000 person-years resulting in the regression coefficient $\ln (\frac{10}{1000}) = - 4.605$ .

The hypothesised treatment effect modification underlying the individual time-to-event outcomes was a beneficial

(H R_{x | z = 1} = 0.8)

, null

(H R_{x | z = 2} = 1.0)

, and harmful

(H R_{x | z = 3} = 1.2)

treatment effect for low, medium, and high levels of

z

, respectively. As such, we set the following parameters:

$β_{1}$ : The beneficial effect of treatment $x$ when $z = 0$ resulted in the regression coefficient $β_{1} = \ln (0.80) = - 0.223$ .

The EM,

z

, increased the mortality rate in the control group:

$β_{2}$ : The log mortality hazard ratio(HR) was set to $β_{2} = \ln (1.2) = 0.182$ for the indicator variable $I (z = 2)$ .

$β_{3}$ : The log mortality HR for the indicator variable $I (z = 3)$ was set to $β_{3} = \ln (1.5) = 0.405$ .

Next, the regression coefficients of the two product terms between the treatment,

x

, and the EM,

z

, were derived from the following assumed interaction mechanism:

$β_{4}$ : Given $H R_{x | z = 2} = 1.0$ , it follows that the regression coefficient of the first product term $x \cdot I (z = 2)$ was set to $β_{4} = \ln (1) - \ln (0.8) = 0.223$ .

$β_{5}$ : Given $H R_{x | z = 3} = 1.2$ , we have that the regression coefficient of the second product term $x \cdot I (z = 3)$ was set to $β_{5} = \ln (1.2) - \ln (0.8) = 0.405$ .

Heterogeneous treatment effect modification: In addition to a common treatment effect modification as outlined above, we considered a heterogeneous (random) treatment effect modification scenario. A distinctive feature of this scenario is that there is a distribution, as opposed to fixed values, of regression coefficients that govern the magnitude and direction of the effect modification underlying all trials included in the IPD meta-analysis. The set of regression coefficients underlying the data for the

j

th study, compactly denoted as

β_{j}^{*} = (β_{1}^{*}, \dots, β_{5}^{*})

, followed a multivariate normal distribution centred around the average of

β

and an identity variance/covariance matrix with a common variability across studies

\sqrt{τ^{2}} \cdot I_{5}

. The regression coefficients underlying the

j

th study denoted by

β_{j}^{*}

were randomly picked from

β_{j}^{*} \sim M V N (β, \sqrt{τ^{2}} \cdot I_{5})

, where the between-study variability was set to

τ = 0.05

for all coefficients. To facilitate a comparison with the common effect scenario, the vector

β

was defined with the identical average parameter values of

β = (- 0.223, 0.182, 0.405, 0.223, 0.405)

3.2. Scenarios

We generated IPD meta-analyses under a common and heterogeneous treatment effect modification with $N = 6$ and $N = 12$ number of trials. The proportion of trials with systematic missing EM was set to 1/6 and 1/3 following the approach presented in Jolani et al.,¹⁰ Audigier et al.¹⁴ The EM was set to missing at a single site independently of all observed covariates. We can reasonably assume that data on the EM might be missing due to administrative, financial, and/or study design reasons that are unrelated to other observed information but related to the study site itself. For all scenarios, the sample size for the studies with complete information was fixed at $n = 500$ participants (with an average of 127 events per study). In total, we simulated eight main scenarios each with two different imputation models according to the number of product terms between the treatment and survival time/censoring indicator included in the imputation model. Each scenario was simulated 1000 times with each 30 imputations. As a reference scenario we considered a complete case analysis.

3.3. Imputation models

In all scenarios, we imputed the systematically missing values for $z$ with two different specifications of the imputation model in the same IPD meta-analysis. Given the operating mechanism between treatment and EM in predicting mortality, a congenial imputation model for $z$ should include main effects and the product terms between the treatment and mortality (both survival time and censoring indicator) in predicting $z$ . We used the Nelson–Aalen estimator of the cumulative hazard rate, $\hat{H} (t)$ , in each study as explained by White et al.²⁶ The imputation model was specified according to equation 1 in the description of CQI with the probability that $θ_{i k}$ is equal to level $k$ of the EM conditionally on the set of observed variables $w_{i}$ . Here, $w_{i}$ is a vector of variables with $w_{i} = (x_{i}, {\hat{H}}_{i} (t), d_{i})$ . The imputation model was fit in studies with complete data on the EM. The first specification of the imputation model with all product terms (Model I: two product terms) imputing the systematic missing variable was defined as:

\ln (\frac{θ_{i k}}{θ_{i 0}} | x_{i}, {\hat{H} (t)}_{i}, d_{i}) = γ_{0 k} + γ_{1 k} x_{i} + γ_{2 k} {\hat{H} (t)}_{i} + γ_{3 k} d_{i} + γ_{4 k} x_{i} \cdot {\hat{H} (t)}_{i} + γ_{5 k} x_{i} \cdot d_{i}

(10)

The imputation model is congenial with the substantive outcome model and should be derived after careful consideration as outlined in White et al.¹⁷

The following imputation model differs from the initial one in a way that it is incomplete. We specified the imputation model omitting all product terms, $x_{i} \cdot {\hat{H} (t)}_{i}$ and $x_{i} \cdot d_{i}$ (Model II: zero product terms).

\ln (\frac{θ_{i k}}{θ_{i 0}} | x_{i}, {\hat{H} (t)}_{i}, d_{i}) = γ_{0 k} + γ_{1 k} x_{i} + γ_{2 k} {\hat{H} (t)}_{i} + γ_{3 k} d_{i}

(11)

3.4. Estimation and statistical inference

Based on the observed and imputed data of the trials included in the IPD meta-analysis, a two-stage approach was used to obtain an estimate of the regression coefficients of the multivariable Weibull survival model. In a two-stage approach, study-specific regression coefficients are first estimated within each study based on imputed data and then combined across studies with a multivariate meta-regression model according to a common or heterogeneous treatment effect modification.²⁷ The inverse variance method and the restricted maximum likelihood were used for the common and heterogeneous (random) treatment effect modification, respectively.²⁴ Rubin’s rules were used to aggregate the estimates of the regression coefficients across imputations. Following standard notations,²⁵ we briefly outline Rubin’s rules for a single parameter. Starting from the $M$ different ${\hat{β}}_{m}$ , the MI point estimate is their average: ${\hat{β}}_{M I} = 1 / M \sum_{m = 1}^{M} {\hat{β}}_{m}$ . The variance of ${\hat{β}}_{M I}$ , $σ_{M I}$ , is estimated by the sum of a within imputation variance, $W_{M I}$ , and a between imputation variance, $B_{M I}$ . A $(1 + 1 / M)$ term is added to account for the finite number of imputations:

\begin{aligned} W_{M I} & = \frac{1}{M} \sum_{m = 1}^{M} {\hat{σ}}_{m}^{2} \\ B_{M I} & = \frac{1}{M - 1} \sum_{m = 1}^{M} ({\hat{β}}_{k} - β_{M I})^{2} \\ σ_{M I} & = W_{M I} + (1 + \frac{1}{M}) B_{M I} \end{aligned}

(12)

Furthermore, a multivariate Wald type test, conducted with a Type I error of 5%, for the hypothesis of no effect modification translated into testing that the two regression coefficients $β_{4}$ and $β_{5}$ of the treatment EM product terms, were jointly equal to zero.

3.5. Performance measures

The principle measures to assess the performance of the two-stage imputation method were bias and coverage.²⁰ Bias was computed as the distance between the average estimated regression coefficients, ( $β_{1}, β_{4}, β_{5}, β_{1} + β_{4}, β_{1} + β_{5}$ ), and the set parameters ( $β_{1} = - 0.223, β_{4} = 0.223, β_{5} = 0.405, β_{1} + β_{4} = 0.000, β_{1} + β_{5} = 0.182$ ). Coverage was computed as the fraction of studies in which the estimated 0.025 and 0.975 quantiles of confidence covered the parameter values. In addition, empirical standard error (ESE) and mean squared error (MSE) were reported.

3.6. Results

The results of the performance of CQI in all scenarios for the common and heterogeneous treatment effect modification are shown in Tables 1 and 2, respectively. The average estimates of the regression coefficients describing the effect of the treatment, $x$ , at different levels of the EM, $z$ , are presented in addition to bias and coverage. Figures 1 and 2 show the different distributions of the parameters under the two specifications of the imputation model with 6 and 12 studies, respectively. Bias was comparably low for both common and heterogeneous treatment effect modification and no stark differences could be pointed out between the number of studies with systematically missing data in the IPD meta-analysis. The ESEs were lower in scenarios with 12 studies whereas more uncertainty was observed in smaller IPD meta-analyses. The complete case analysis was unbiased and coverage of the effect estimates was close to 95%. The following two sections describe in more detail the results for the common and heterogeneous treatment effect modification.

Figure 1.

Sampling distribution of the estimated mortality hazard ratios conferred by the treatment at different levels of the low ( $z = 1$ ), null ( $z = 2$ ), and high ( $z = 3$ ) effect modifiers (EMs) under different complexities of the imputation model for common and heterogeneous treatment effect modification. Both scenarios, common and heterogeneous treatment effect modification, were complied of 6 studies with 1/6 and 1/3 of the studies with systematically missing information on the EM. The sample size of all studies was set to 500. Model 1 (two product terms): $x_{i} \cdot \hat{H} (t)_{i}$ and $x \cdot d_{i}$ . Model 2 (zero product terms): Omitting $x \cdot \hat{H} (t)$ and $x \cdot d_{i}$ .

Figure 2.

Sampling distribution of the estimated mortality hazard ratios conferred by the treatment at different levels of the low ( $z = 1$ ), null ( $z = 2$ ), and high ( $z = 3$ ) effect modifiers(EMs) under different complexities of the imputation model for common and heterogeneous treatment effect modification. Both scenarios, common and heterogeneous treatment effect modification, were complied of 12 studies with 1/6 and 1/3 of the studies with systematically missing information on the EM. The sample size of all studies was set to 500. Model 1 (two product terms): $x_{i} \cdot \hat{H} (t)_{i}$ and $x \cdot d_{i}$ . Model 2 (zero product terms): Omitting $x \cdot \hat{H} (t)$ and $x \cdot d_{i}$ .

Table 1.

Log hazard ratios (average and performance measure) for the treatment effect at different levels of the EM by the complexity of the imputation models and by the size of the studies with systematic missing data on the EM.

	Complete case analysis					Model I: Two product terms					Model II: Zero product terms
Parameters	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$
6 studies with 1 systematically missing (1/6 missing)
Average	$- 0.217$	0.221	0.396	0.003	0.178	$- 0.217$	0.219	0.394	0.002	0.176	$- 0.186$	0.184	0.330	$- 0.002$	0.143
Bias	0.006	$- 0.002$	$- 0.010$	0.003	$- 0.004$	0.006	$- 0.004$	$- 0.012$	0.002	$- 0.006$	0.037	$- 0.039$	$- 0.076$	$- 0.002$	$- 0.039$
Coverage¹	94.7	95.3	95.0	94.7	93.6	94.5	95.1	94.9	94.9	93.6	96.7	97.8	96.8	96.8	95.7
ESE	0.143	0.188	0.210	0.127	0.158	0.139	0.187	0.209	0.121	0.154	0.124	0.157	0.175	0.110	0.135
MSE	0.021	0.035	0.044	0.016	0.025	0.019	0.035	0.044	0.015	0.024	0.017	0.026	0.036	0.012	0.020
6 studies with 2 systematically missing (1/3 missing)
Average	$- 0.225$	0.230	0.408	0.005	0.184	$- 0.221$	0.227	0.402	0.006	0.181	$- 0.157$	0.153	0.272	$- 0.004$	0.115
Bias	$- 0.001$	0.006	0.003	0.005	0.001	0.002	0.004	$- 0.003$	0.006	$- 0.001$	0.066	$- 0.070$	$- 0.134$	$- 0.004$	$- 0.068$
Coverage²	95.1	95.6	95.2	94.9	95.3	93.9	94.4	94.8	93.8	95.5	96.3	99.0	97.1	98.2	97.5
ESE	0.158	0.212	0.231	0.140	0.168	0.149	0.209	0.228	0.128	0.157	0.115	0.142	0.154	0.103	0.116
MSE	0.025	0.045	0.053	0.020	0.028	0.022	0.044	0.052	0.016	0.024	0.018	0.025	0.042	0.011	0.018
12 studies with 2 systematically missing (1/6 missing)
Average	$- 0.216$	0.214	0.393	$- 0.002$	0.176	$- 0.215$	0.212	0.391	$- 0.003$	0.175	$- 0.185$	0.178	0.327	$- 0.007$	0.142
Bias	0.007	$- 0.009$	$- 0.013$	0.002	$- 0.006$	0.008	$- 0.011$	$- 0.015$	$- 0.003$	$- 0.007$	0.038	$- 0.045$	$- 0.078$	$- 0.007$	$- 0.040$
Coverage³	96.5	94.6	96.1	94.2	96.5	95.7	94.0	95.9	93.7	95.7	97.0	96.2	96.7	96.3	95.8
ESE	0.097	0.130	0.138	0.089	0.104	0.094	0.130	0.137	0.085	0.102	0.084	0.109	0.115	0.077	0.090
MSE	0.009	0.017	0.019	0.008	0.011	0.009	0.017	0.019	0.007	0.010	0.008	0.014	0.019	0.006	0.010
12 studies with 4 systematically missing (1/3 missing)
Average	$- 0.220$	0.218	0.394	$- 0.003$	0.173	$- 0.217$	0.215	0.389	$- 0.002$	0.172	$- 0.156$	0.145	0.262	$- 0.011$	0.106
Bias	0.003	$- 0.006$	$- 0.012$	$- 0.003$	$- 0.009$	0.006	$- 0.008$	$- 0.016$	$- 0.002$	$- 0.011$	0.067	$- 0.078$	$- 0.143$	$- 0.011$	$- 0.076$
Coverage⁴	95.4	95.2	95.8	95.2	95.4	95.4	94.5	95.1	93.6	94.0	95.0	97.9	94.0	97.8	95.2
ESE	0.110	0.146	0.162	0.098	0.123	0.102	0.144	0.158	0.089	0.115	0.079	0.097	0.107	0.071	0.086
MSE	0.012	0.021	0.026	0.010	0.015	0.010	0.021	0.025	0.008	0.013	0.011	0.016	0.032	0.005	0.013

Results are presented for 1/6 and 1/3 of the studies with systematically missing data. All simulated data was generated under a common treatment effect modification with parameter values equal to: $β_{1} = - 0.223, β_{4} = 0.223, β_{5} = 0.405, β_{1} + β_{4} = 0.000, β_{1} + β_{5} = 0.182$ . Each simulated IPD meta-analysis was analysed with a two-stage multivariate Weibull survival model. Each scenario was simulated 1000 times with each with 30 imputations using conditional quantile imputation, ESE, MSE. Model I (two product terms): $x_{i} \cdot \hat{H} (t)_{i}$ and $x \cdot d_{i}$ . Model II (zero product terms): Omitting $x \cdot \hat{H} (t)$ and $x \cdot d_{i}$ . MCE was below a maximum of 0.007 for all parameters except coverage; ¹:MCE < 0.891; ²:MCE < 1.042; ³:MCE < 1.024; ⁴:MCE < 1.550. IPD: individual participant data; EM: effect modifier; ESE: empirical standard error; MSE: mean squared error; MCE: Monte Carlo error.

Table 2.

	Complete case analysis					Model I: Two product terms					Model II: Zero product terms
Parameters	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$	${\hat{β}}_{1}$	${\hat{β}}_{4}$	${\hat{β}}_{5}$	${\hat{β}}_{1} + {\hat{β}}_{4}$	${\hat{β}}_{1} + {\hat{β}}_{5}$
6 studies with 1 systematically missing (1/6 missing)
Average	$- 0.221$	0.220	0.399	$- 0.001$	0.179	$- 0.221$	0.219	0.397	$- 0.002$	0.176	$- 0.189$	0.183	0.333	$- 0.007$	0.143
Bias	0.003	$- 0.003$	$- 0.006$	$- 0.001$	$- 0.003$	0.002	$- 0.004$	$- 0.009$	$- 0.002$	$- 0.006$	0.034	$- 0.041$	$- 0.073$	$- 0.007$	$- 0.039$
Coverage¹	94.8	95.6	94.9	94.1	94.9	94.6	95.3	94.7	93.7	95.2	96.5	98.0	96.7	96.2	97.3
ESE	0.143	0.186	0.208	0.130	0.159	0.138	0.186	0.207	0.125	0.153	0.123	0.155	0.174	0.113	0.135
MSE	0.020	0.035	0.043	0.017	0.025	0.019	0.034	0.043	0.016	0.023	0.016	0.026	0.035	0.013	0.020
6 studies with 2 systematically missing (1/3 missing)
Average	$- 0.219$	0.222	0.405	0.003	0.186	$- 0.217$	0.219	0.400	0.002	0.183	$- 0.152$	0.146	0.267	$- 0.007$	0.114
Bias	0.004	$- 0.001$	$- 0.000$	0.003	0.004	0.006	$- 0.004$	$- 0.006$	0.002	0.000	0.070	$- 0.077$	$- 0.138$	$- 0.007$	$- 0.068$
Coverage²	94.9	94.3	94.8	93.1	95.6	94.9	93.0	92.5	91.5	94.5	97.5	98.7	96.9	97.6	97.7
ESE	0.162	0.220	0.237	0.150	0.174	0.150	0.219	0.235	0.140	0.166	0.116	0.147	0.158	0.111	0.122
MSE	0.026	0.048	0.056	0.022	0.030	0.023	0.048	0.055	0.020	0.028	0.018	0.028	0.044	0.012	0.019
12 studies with 2 systematically missing (1/6 missing)
Average	$- 0.219$	0.220	0.400	0.001	0.181	$- 0.218$	0.219	0.398	0.001	0.180	$- 0.187$	0.184	0.334	$- 0.003$	0.147
Bias	0.004	$- 0.003$	$- 0.005$	0.001	$- 0.001$	0.005	$- 0.004$	$- 0.008$	0.001	$- 0.002$	0.036	$- 0.040$	$- 0.072$	$- 0.003$	$- 0.035$
Coverage³	96.6	96.2	96.4	95.2	95.4	96.2	96.6	95.9	95.0	94.6	96.6	97.8	96.5	96.8	96.4
ESE	0.093	0.126	0.019	0.088	0.107	0.091	0.126	0.138	0.085	0.105	0.082	0.106	0.116	0.077	0.092
MSE	0.009	0.016	0.140	0.008	0.011	0.008	0.016	0.019	0.007	0.011	0.008	0.013	0.019	0.006	0.010
12 studies with 4 systematically missing (1/3 missing)
Average	$- 0.225$	0.223	0.407	$- 0.002$	0.181	$- 0.220$	0.220	0.401	$- 0.000$	0.180	$- 0.158$	0.148	0.271	$- 0.009$	0.113
Bias	$- 0.002$	$- 0.000$	0.001	$- 0.002$	$- 0.001$	0.002	$- 0.003$	$- 0.004$	$- 0.000$	$- 0.002$	0.065	$- 0.075$	$- 0.135$	$- 0.009$	$- 0.069$
Coverage⁴	95.6	94.9	96.0	94.3	94.9	94.5	93.4	94.9	93.1	93.2	95.9	97.6	93.8	97.7	95.7
ESE	0.111	0.152	0.166	0.103	0.126	0.102	0.151	0.164	0.095	0.119	0.078	0.102	0.111	0.076	0.089
MSE	0.012	0.023	0.028	0.011	0.016	0.010	0.023	0.027	0.009	0.014	0.010	0.016	0.030	0.006	0.013

Results are presented for 1/6 and 1/3 of the studies with systematically missing data. All simulated data was generated under a heterogeneous treatment effect modification with the average parameter values equal to: $β_{1} = - 0.223, β_{4} = 0.223, β_{5} = 0.405 β_{1} + β_{4} = 0.000, β_{1} + β_{5} = 0.182$ . Each simulated IPD meta-analysis was analysed with a two-stage multivariate Weibull survival model. Each scenario was simulated 1000 times with each with 30 imputations using conditional quantile imputation, ESE, MSE. Model I (two product terms): $x_{i} \cdot \hat{H} (t)_{i}$ and $x \cdot d_{i}$ . Model II (zero product terms): Omitting $x \cdot \hat{H} (t)$ and $x \cdot d_{i}$ . MCE was below a maximum of 0.007 for all parameters except coverage; ¹:MCE < 0.828; ²:MCE < 1.063; ³:MCE < 1.024; ⁴:MCE < 1.556. IPD: individual participant data; EM: effect modifier; ESE: empirical standard error; MSE: mean squared error.

3.6.1. Common treatment effect modification

The absolute value of the bias after using CQI to impute the systematically missing values for the EM was less than 0.016 for all effect estimates for the specification of the imputation model with two product terms. Bias was comparably low for scenarios with 6 and 12 studies, regardless of the number of studies with systematically missing data on the EM. The largest bias with the congenial imputation model across all scenarios was found for the scenario with 12 studies when 1/3 of the studies had systematically missing data on the EM with 0.006, $- 0.002$ , and $- 0.011$ for the beneficial, null, and harmful treatment effect of the EM, respectively. Specifying the imputation with zero product terms led to larger biases in all effect estimates throughout the different scenarios. Here, the largest bias was found for 12 studies when 1/3 of the studies had systematically missing values on the EM, at around 0.067, $- 0.011$ , and $- 0.076$ for the beneficial, null, and harmful treatment effect of the EM, respectively. Figures 1 and 2 show the difference in sampling distribution of the different levels of the EM between the two specifications of the imputation model. The ability of the IPD meta-analysis to differentiate between levels of the EM decreased drastically when the imputation model included only main effects of treatment, time, and censoring indicator. Here, the magnitude of the effect modification was close to 0, making it difficult for the statistical test to reject the null hypothesis of no effect modification.

The coverage of the treatment effect estimates were close to 95% in all scenarios. Coverage was at 95.7%, 93.7%, and 95.7% for the beneficial, null, and harmful treatment effect of the EM, respectively, for 12 studies when a sixth of studies had systematic missing values on the EM. With a third of the studies having systematically missing data, the coverage for the same scenario remained very similar at 95.4%, 93.6%, and 94.0% for the beneficial, null, and harmful treatment effect of the EM, respectively. The small differences in coverage are explained by random variability in the estimates. When the imputation model included no product terms, coverage was 95.0%, 97.8%, and 95.2% for the three levels of the EM, respectively.

Notably, the precision in all effect estimates increased after using MI to retain trials with systematically missing values on the EM. The largest gains in precision were made when a higher proportion of studies with missing data was imputed. This is evident in the consistently smaller ESE and corresponding MSE between the complete case analysis and MI with the congenial specification of the imputation model. For example, the ESE for the complete case analysis in the scenario of 12 trials with 4 trials having systematically missing data was 0.110, 0.098, and 0.123 for the low, null, and harmful treatment effect, respectively. After retaining the four trials in the analysis, the ESE was reduced to 0.102, 0.089, and 0.115, respectively. In addition, the MSE was lower (about 50%) for the scenarios with 12 studies compared to the scenario with only six studies, regardless of the number of studies with systematically missing data on the EM. Table 1 contains all values for the scenarios of the common treatment effect modification.

3.6.2. Heterogeneous treatment effect modification

Overall, bias was lower compared to the scenarios of common treatment effect modification. In all scenarios simulated under a heterogeneous treatment EM, the absolute value of the bias remained less than 0.009 for all effect estimates under the congenial specification of the imputation model. There were no stark differences in the magnitude of bias between the scenarios. With 12 studies, when 1/3 of the studies were imputed, the bias for all levels of the EM remained below 0.002 under the specification of the imputation model with two product terms. Under the incomplete (zero product terms) specification of the imputation model, large bias was introduced in all parameters in all scenarios. The absolute value of the bias increased up to 0.069 for the estimates of the EM. Here, with 1/3 of the studies having systematically missing data on the EM, bias was highest for the additional effect of the harmful EM ( $β_{5}$ ) with $- 0.138$ and $- 0.135$ for 6 and 12 studies, respectively.

Coverage remained close to 95% in most scenarios of the heterogeneous treatment effect modification. For 6 studies with 1/3 of the studies being imputed, the coverage was at 94.9%, 91.5%, and 94.5% for the beneficial, null, and harmful treatment effect of the EM. For the same scenario with 12 studies, coverage levels were at 94.5%, 93.1%, and 93.2%, respectively. Here, the coverage of the interaction terms itself was marginally closer to the nominal level at 94.5% and 95.1%.

Similar to the common treatment effect modification, precision improved when using MI to retain trials with systematically missing values on the EM, most notable in the smaller ESE and MSE values compared to the complete case analysis. Additionally, the MSE was 50% lower in the scenarios with 12 studies, compared to those with 6.

In summary, in all scenarios under the specification of the congenial imputation model, no substantial bias was found in estimates at all levels of the EM. The magnitude of the bias was not affected with an increase in the number of studies with systematic missing values on the EM. The specification of the imputation model with zero product terms always resulted in an increase in bias of the estimates, regardless of the number of studies with systematically missing data on the EM. Using a two-stage imputation approach to impute systematically missing EM values, we observed an improvement in the precision of effect estimates in all scenarios under the congenial specification of the imputation model compared to the complete case analysis.

4. The effect of PORT on survival at different stages of the disease

In this section, we illustrate the use of CQI for a systematically missing EM in an IPD meta-analysis of trials. The main question was about the effectiveness of PORT in patients with completely resected non-small cell lung cancer at different stages of the disease. A multivariate IPD meta-analysis of PORT versus surgery alone on survival in patients with resected non-small-cell lung cancer at different stages of the disease was performed.²⁸ The outcome in all trials was time from randomisation to death from all causes or censoring, whichever came first. The disease stage of the patients was measured in three different levels ranging from stage I to III. The 10 trials included in this study had a varying sample size ranging from 69 to 316 participants per trial with a total of 1642 participants and 1082 deaths across all 10 trials. A common treatment effect modification underlying the trials was estimated with a two-stage multivariable Cox proportional hazard regression model with $x$ being the treatment and $z$ the EM (disease stage):

h (t ∣ x_{i}, z_{i}) = h_{0} (t) \cdot e^{β_{1} x_{i}} + β_{2} I (z_{i} = 2) + β_{3} I (z_{i} = 3) + β_{4} x_{i} \cdot I (z_{i} = 2) + β_{5} x_{i} \cdot I (z_{i} = 3)

(13)

where

h_{0} (t)

is the baseline hazard function. The function

h_{0} (t)

indicates some unspecified trend of the hazard over time when all predictors are zero. Seven out of the 10 trials had complete information on the disease stage of patients (the EM). However, three trials with a sample size of 110, 169, and 176 had no information on the disease stage, i.e., systematically missing information. Different from the settings in the simulation study, age and sex of the patients were measured and could also be used for the imputation model. CQI was used to impute the systematically missing data on the disease stage in three trials based on the trials with complete information on the stage of the disease of patients (30 imputations). This allowed us to retain all 10 trials in the analysis and not lose valuable treatment-outcome information. In the seven studies with complete data on disease stage, the imputation model was specified in accordance to the final outcome model, incorporating the effect modification between the time-varying and indicator variable for death:

\begin{aligned} \ln (\frac{θ_{i k}}{θ_{i 0}} | x_{i}, {\hat{H} (t)}_{i}, d_{i}, a g e_{i}, s e x_{i}) & = γ_{0 k} + γ_{1 k} x_{i} + γ_{2 k} {\hat{H} (t)}_{i} + γ_{3 k} d_{i} \\ + γ_{4 k} x_{i} \cdot {\hat{H} (t)}_{i} + γ_{5 k} x_{i} \cdot d_{i} + γ_{6 k} a g e_{i} + γ_{7 k} s e x_{i} \end{aligned}

(14)

The data used in this paper is a publicly available, real-data based simulated example with similar characteristics to the PORT data and is only used for the illustration of the presented imputation method.²⁹

4.1. Results

In this sample of data, only using the seven studies with complete information on disease stage, the statistical test did not provide a strong indication against the hypothesis of a homogeneous effect of PORT on mortality across stages of the disease ( $χ^{2} = 5.00$ , $p = 0.082$ ). The estimated mortality HRs of PORT were 1.00, 1.02, and 1.43, for the first, second, and third stage of the disease, respectively. At the latest stage of the disease, the effect of PORT on survival was estimated with a 43% increase in the mortality hazard rate (HR $=$ 1.43, 95% CI: 1.11, 1.85, $z = 2.73$ , $p = 0.006$ ).

Similar to the complete case analysis, the statistical test when using all 10 studies after using CQI to impute the systematically missing values for disease stage in three trials did not provide a clear indication against the hypothesis of a homogeneous effect of PORT on mortality across stages of the disease ( $χ^{2} = 5.53$ , $p = 0.063$ ). For the first, second, and third stage of the disease the estimated mortality HRs of PORT were 1.03, 1.03, and 1.45, respectively. The mortality HR increased by 45% ( $z = 3.14$ , $p = 0.002$ ) when using PORT compared to surgery alone at disease stage III (HR $=$ 1.45, 95% CI: 1.15, 1.83) (see Table 3).

Table 3.
Mortality HRs, 95% CIs, and SEs conferred by postoperative radiotherapy at different stages of the disease in complete case and MI datasets.

7 studies (Complete case) 10 studies (Using MI)

N $=$ 1187, deaths $=$ 785 N $=$ 1642, deaths $=$ 1082

Mortality HR (95% CI) HR 95% CI SE HR 95% CI SE

Stage I 1.00 (0.78, 1.27) 0.13 1.03 (0.82, 1.28) 0.12

Stage II 1.02 (0.80, 1.29) 0.12 1.03 (0.83, 1.28) 0.11

Stage III 1.43 (1.11, 1.85) 0.19 1.45 (1.15, 1.83) 0.17

	7 studies (Complete case)	10 studies (Using MI)
Stage I	1.00	(0.78, 1.27)	0.13	1.03	(0.82, 1.28)	0.12
Stage II	1.02	(0.80, 1.29)	0.12	1.03	(0.83, 1.28)	0.11
Stage III	1.43	(1.11, 1.85)	0.19	1.45	(1.15, 1.83)	0.17

The individual participant data meta-analysis included 10 trials, 7 with complete information on the disease stage, 3 with systematically missing data on disease stage. Conditional quantile imputation was used to impute the systematically missing data based on the trials with complete information. Data were analysed with a two-stage multivariable Cox-regression model. MI: multiple imputation; HR: hazard ratio; CI: confidence interval; SE: standard error.

4.2. Interpretation

The described example of imputing the three trials with systemically missing values on disease stage suggests a beneficial use of MI in this scenario. While the empirical results of the analysis did not change substantially compared to the complete case analysis, all studies are retained in the analysis without introducing substantial bias and pointing towards a harmful effect of PORT on survival in the later stage of the disease. By using MI to retain the three trials with systematically missing values on disease stage, we were able to include an additional 455 participants including 297 mortality cases in the analysis. While no substantial differences in the mortality HRs between the complete case analysis and the analysis with MI were estimated, the precision in the effect estimates at all levels of disease stage increased, without distorting the association that was estimated in the complete case analysis. This is in line with the results presented from the simulation study. Overall, in the example of the effect of PORT on survival at different stages of the disease, using MI to retain three trials with systematically missing values on disease stage indicates a worthwhile use to increase the generalisbility and precision of clinical effect measures without introducing bias.

5. Discussion

This simulation study evaluated a two-stage imputation method based on conditional quantiles to assess its performance on retaining studies with systematically missing EMs in IPD meta-analyses. We evaluated the feasibility of imputing systematically missing EM in IPD meta-analysis with a limited number of trials (6 and 12) under a common and heterogeneous treatment effect modification. The results demonstrated that the bias for all levels of the EM was considerably low for the common and heterogeneous treatment effect modification. Compared to the complete case analysis, using the two-stage imputation approach to retain all trials with missing data in the analysis improved the precision of effect estimates in all scenarios under a congenial specification of the imputation model. The bias increased for all scenarios when the imputation model was incompletely specified, i.e., missing important product terms.

5.1. Performance of CQI

First, the proposed approach, CQI, indicated no substantial difference in bias between smaller and larger IPD meta-analysis nor substantial differences when the number of studies with systematically missing data on the EM doubled (i.e., from one sixths to one third). In fact, the performance showed a negligible bias for both common and heterogeneous treatment effect modification IPD meta-analysis that is comparable with biases reported for other MI methods.¹⁴ The average effect estimates after using CQI were comparable with those of the complete case analysis. Even though a complete case analysis can be justified when assuming missing completely at random (MCAR), it is suboptimal when the fraction of studies with missing EM is large. As we demonstrated in this simulation study, this can lead to a reduction in the precision of effect estimates and substantial loss of data. Bias increased substantially when the imputation model was incompletely specified with zero product terms (i.e., not congenial with the substantive outcome model). This is in line with previous work that has demonstrated the importance of correctly specifying the imputation model, whereby failing to do so leads to biased inference.¹⁶ Sensitivity analyses are often used in practice to explore the implications of different imputation models.^30–32

Second, coverage was close to the nominal level for all scenarios that were assessed. Some small deviations from the desired 95% were observed in scenarios with a larger number of studies and when a third of the studies had systematically missing data on the EM. In the same scenarios, similar deviations from 95% were seen for the complete case analysis. This is in line with simulation results presented by Resche-Rigon and White¹¹ indicating slight undercoverage in the presence of only systematically missing data.

Third, we showed that using a two-stage imputation approach in the context of systematically missing data in IPD meta-analysis improved the precision of effect estimates compared to the complete case analysis. That is, lower ESE and MSE across all scenarios were estimated reflecting the better use made of the data. In addition, by including all trials in the analysis, the generalisbility of study findings potentially increases, although it is difficult to quantify this numerically. In summary, CQI resulted in an analysis with (1) low bias, (2) coverage close to the nominal 95% level, and (3) small variance.

5.2. Assumptions

Despite satisfactory performance benefits of the proposed method, certain assumptions need to be discussed. First, we assumed that the underlying treatment effect at different levels of the EM is consistent across all trials. This assumption supports the interpretation of observed variations in treatment effects across levels of the EM as true interactions, as opposed to artefacts from trial-level heterogeneity, such as differences in study design, populations, or implementation. Second, we assumed that all trials collected similar information on other covariates used to inform the imputation model. This assumption is more likely to be satisfied in prospective IPD meta-analyses, where common data collection, analysis plans, and harmonisation strategies are implemented, thereby reducing the risk of systematic missing values.³³ In retrospective IPD meta-analyses, however, data collection may have occurred at different time points and protocols may vary across studies, increasing the risk of systematic missingness.³⁴ Related to this, we also assumed that the EM has the same number of categories (or levels) across trials and was measured consistently across studies. In IPD meta-analyses, harmonisation of covariates is a crucial step prior to implementing any MI approach. Therefore, EMs with different distributions across trials must be harmonised before applying the two-stage imputation approach. Last, we addressed systematic missing data under a MCAR assumption. As previously mentioned, a complete case analysis is theoretically justifiable, albeit suboptimal due to the significant loss of data. Under a MAR assumption, the approach to impute systematically missing values of the EM would not differ. While the results may become more sensitive to the choice of predictors included in the imputation model assuming MAR, simulation studies have shown small performance differences between MCAR and MAR in multilevel data settings for sporadically and systematically missing data.^11,35

5.3. Limitations

Despite demonstrating an approach for imputing a systematically missing EM in IPD meta-analyses with a limited number of trials, this study is subject to a number of limitations. First, we aimed at simulating simplified, yet complex enough, scenarios that are related to realistic challenges in IPD meta-analysis, whilst remaining accessible for researchers that are faced with such problems. While there is an abundance of scenarios that were not considered, we are confident that the chosen scenarios can give an intuition about the performance of the CQI method and observation of a general trend of when MI is worth considering in IPD meta-analysis to impute systematically missing EM. Second, we solely focussed on the impact of systematically missing data and did not consider scenarios with sporadically and systematically missing data at the same time. Third, in all our simulations we used 30 imputations. We did not test whether the performance, in particular the coverage, changes with an increasing number of imputations due to reasonable computation time. Last, we did not consider larger degrees of heterogeneity across studies in our simulations due to the focus on IPD meta-analyses with a limited number of studies.

5.4. Phases of methodological development

Based on the phases of methodological research in biostatistics according to Heinze et al.,³⁶ this work can be categorised as part of phase I to II developments of a two-stage imputation procedure for systematic missing covariates in IPD meta-analysis in line with previous work.¹¹ Facilitating software is available in Stata and provided in the resources linked to this study. Further analysis is needed to test the approach in a variety of settings such as different outcomes, varying degrees of heterogeneity and DGMs. In particular, future continuation on the refinement of this approach can be directed towards:

(1)
Extending the approach to multivariate imputation as proposed in Resche-Rigon and White:¹¹ Further software implementation of CQI into MI with chained equation is a logical step to make it more widely applicable in scenarios with sporadically and systematically missing values. An extension of CQI to continuous systematic missing covariates using quantile regression has been recently presented.³⁷
(2)
Integrating random effects into the two-stage imputation procedure to increase heterogeneity between imputed datasets: Given the small number of studies in this analysis, a random-effects model would provide a poor estimate of the variability between studies. However, using a common-effects model to derive average regression coefficients for assigning values to the missing EM can result in more homogeneous imputed datasets.

5.5. Conclusion

This simulation study presented and evaluated the use of a two-stage imputation procedure - CQI - to impute systematically missing EMs in IPD meta-analyses with a limited number of trials. The absolute bias for common and heterogeneous effect IPD meta-analyses was less than 0.016 and 0.007, respectively, with coverage close to its nominal value across all levels of the EM. In addition, CQI improved the precision of pooled effect estimates compared to a complete case analysis excluding trials with systematically missing values on the EM. An incomplete specification of the imputation model resulted in biased inference even if the proportion of studies with systematically missing data was small.

Footnotes

Acknowledgments

The authors would like to thank Stephanie G. Pitt (Karolinska Institutet) for proofreading the manuscript.

Data availability

The current version of CQI for discrete variables can be downloaded as a Stata 18 package. In addition to the package, the DGM used for this manuscript and the code using CQI for the applied example can be accessed here: .

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethical approval

Not applicable. This study only used simulated data.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the National Infrastructure NEAR, supported by the Swedish Research Council (Dnrs 2017-00639 and 2021-00178).

ORCID iD

Robert Thiesmeier

References

Riley

Debray

Fisher

. Individual participant data meta-analysis to examine interactions between treatment effect and participant-level covariates: Statistical recommendations for conduct and planning. Stat Med 2020; 39: 2115–2137.

Tierney

Riley

Smith

. Rationale for embarking on an ipd meta-analysis project. In: Individual participant data meta-analysis: A handbook for healthcare research, 2021. pp.7–19. John Wiley & Sons Ltd. ISBN 9781119333784.

Godolphin

Marlin

Cornett

. Use of multiple covariates in assessing treatment-effect modifiers: A methodological review of individual participant data meta-analyses. Res Syn Meth 2023; 15: 107–116.

Lambert

Sutton

Abrams

. A comparison of summary patient-level covariates in meta-regression with individual patient data meta-analysis. J Clin Epidemiol 2002; 55: 86–94.

Brookes

Whitely

Egger

. Subgroup analyses in randomized trials: Risks of subgroup-specific analyses; power and sample size for the interaction test. J Clin Epidemiol 2004; 57: 229–236.

Leon

Heo

. Sample sizes required to detect interactions between two binary fixed-effects in a mixed-effects linear regression model. Comput Stat Data Anal 2009; 53: 603–608.

Stewart

Tierney

. To ipd or not to ipd? Stat Med 2002; 25: 76–97.

Riley

Fisher

. Using ipd meta-analysis to examine interactions between treatment effect and participant-level covariates. In: Individual participant data meta-analysis: A handbook for healthcare research, 2021, pp.163–198. John Wiley & Sons Ltd. ISBN 9781119333784.

Resche-Rigon

White

Bartlett

. Multiple imputation for handling systematically missing confounders in meta-analysis of individual participant data. Stat Med 2013; 32: 4890–4905.

10.

Jolani

Debray

Koffijberg

. Imputation of systematically missing predictors in an individual participant data meta-analysis: A generalized approach using mice. Stat Med 2015; 34: 1841–1863.

11.

Resche-Rigon

White

. Multiple imputation by chained equations for systematically and sporadically missing multilevel data. Stat Methods Med Res 2018; 27: 1634.

12.

Quartagno

Carpenter

. Multiple imputation for ipd meta-analysis: Allowing for heterogeneity and studies with missing covariates. Stat Med 2016; 35: 2938–2954.

13.

Kunkel

Kaizar

. A comparison of existing methods for multiple imputation in individual participant data meta-analysis. Stat Med 2017; 36: 3507–3532.

14.

Audigier

White

Jolani

. Multiple imputation for multilevel data with continuous and binary variables. J Stat Softw 2018; 33: 160–183.

15.

Meng

. Multiple imputation with uncongenial sources of input (with discussion). Stat Sci 1994; 9: 538–573.

16.

Wood

White

Royston

. How should variable selection be performed with multiply imputed data? Stat Med 2008; 27: 3227–3246.

17.

White

Royston

Wood

. Multiple imputation using chained equations: Issues and guidance for practice. Stat Med 2011; 30: 377–399.

18.

D’Agostino McGowan

Lotspeich

SSAH

. The “why” behind including “y” in your imputation model. Stat Methods Med Res 2024; 33: 996–1020.

19.

Steyerberg

Nieboer

Debray

, et al. Assessment of heterogeneity in an individual participant data meta-analysis of prediction models: An overview and illustration. Stat Med 2918; 38: 4290–4309.

20.

Morris

White

Crowther

. Using simulation studies to evaluate statistical methods. Stat Med 2019; 38: 2074–2102.

21.

White

Pham

Quartagno

, et al. How to check a simulation study. Int J Epidemiol 2024; 53: dyad134.

22.

Muñoz

Rueda

. New imputation methods for missing data using quantiles. J Comput Appl Math 2009; 232: 305–317.

23.

Bottai

Zhen

. Multiple imputation based on conditional quantile estimation. Epidemiol Biostat Pub Health 2013; 1: 1–18.

24.

Hartung

Knapp

Sinha

. Statistical meta-analysis with applications. Hoboken, NJ: John Wiley & Sons, 2008.

25.

Rubin

. Multiple Imputation for Nonresponse in Surveys. Hoboken, NJ, USA: John Wiley & Sons, Inc., 1987.

26.

White

Royston

. Imputing missing covariate values for the cox model. Stat Med 2009; 28: 1982–1998.

27.

Riley

Debray

Morris

, et al. The two-stage approach to ipd meta-analysis. In: Individual participant data meta-analysis: A handbook for healthcare research, 2021, pp.87–125. John Wiley & Sons Ltd. ISBN 9781119333784.

28.

Burdett

Stewart

. PORT meta-analysis group postoperative radiotherapy in non-small-cell lung cancer: Update of an individual patient data meta-analysis. Lung Cancer 2005; 47: 81–83.

29.

Fisher

. Two-stage individual participant data meta-analysis and generalized forest plots. Stata J 2015; 15: 369–396.

30.

Noghrehchi

Stoklosa

Penev

Sea

. Selecting the model for multiple imputation of missing data: Just use an ic! Stat Med 2021; 40: 2467–2497.

31.

Nguyen

Carlin

Lee

. Model checking in multiple imputation: An overview and case study. Emerg Themes Epidemiol 2017; 14: 1–12.

32.

Carpenter

Kenward

White

. Sensitivity analysis after multiple imputation under missing at random: A weighting approach. Stat Methods Med Res 2007; 16: 259–275.

33.

Seidler

Hunter

Cheyne

, et al. A guide to prospective meta-analysis. BMJ 2019; 367: l5342.

34.

Debray

Snell

Quartagno

, et al. Dealing with missing data in an ipd meta-analysis. Indiv Partic Data Meta-Anal: Handb Health Res 2021; 499–524. doi:10.1002/9781119333784.ch18

35.

Andridge

. Quantifying the impact of fixed effects modeling of clusters in multiple imputation for cluster randomized trials. Biomet J 2011; 53: 57–74.

36.

Heinze

Boulesteix

Kammer

, et al. Phases of methodological research in biostatistics–building the evidence base for new methods. Biomet J 2024; 66: 2200222.

37.

Thiesmeier

Bottai

Orsini

. Systematically missing data in distributed data networks: Multiple imputation when data cannot be pooled. J Stat Comput Simul 2024; 94: 3807–3825.

	7 studies (Complete case)			10 studies (Using MI)
	N $=$ 1187, deaths $=$ 785			N $=$ 1642, deaths $=$ 1082
Mortality HR (95% CI)	HR	95% CI	SE	HR	95% CI	SE
Stage I	1.00	(0.78, 1.27)	0.13	1.03	(0.82, 1.28)	0.12
Stage II	1.02	(0.80, 1.29)	0.12	1.03	(0.83, 1.28)	0.11
Stage III	1.43	(1.11, 1.85)	0.19	1.45	(1.15, 1.83)	0.17