Standard and reference-based conditional mean imputation

The ICH E9(R1) addendum on estimands stresses the importance of a precise description of the targeted treatment effect reflecting the clinical question posed by the trial objectives [1]. One important attribute of an estimand is a list of possible intercurrent events (ICEs), i.e. events occurring after treatment initiation that affect either the interpretation or the existence of the measurements associated with the clinical question of interest, and the definition of appropriate strategies to deal with ICEs. The two most relevant strategies for the purpose of this article are the hypothetical strategy and the treatment policy strategy. For the hypothetical strategy, a scenario is envisaged in which the ICE would not occur. For the treatment policy strategy, the treatment effect in the presence of the ICEs is targeted and analyses under a treatment policy strategy include all observed outcomes regardless whether the subject had an ICE or not.

The ICH E9(R1) addendum distinguishes between ICEs and missing data. However, there are many connections between them. Under the hypothetical scenario, outcome values after the ICE are not directly observable and treated using models for missing data. Under a treatment policy scenario, the outcome assessments after the ICE are observable but may be incomplete because study drop out cannot be prevented. In addition, intermittently missing outcome assessments may also occur due to missed visits.

For more than a decade, a likelihood-based mixed model for repeated measures (MMRM) approach has been the de facto standard for the primary analysis of continuous longitudinal endpoints in pharmaceutical clinical trials [2]. MMRM approaches typically model the longitudinal outcomes at scheduled visits using a multivariate normal model with the visit treated as a categorical variable in each treatment group, adjustment for baseline covariates, and an unstructured covariance matrix. This provides valid treatment effect estimates if the missing data are missing at random (MAR) after accounting for baseline covariates, the assigned treatment group, and observed outcomes and if the complete data model is correctly specified. Such a MAR assumption is often a good starting point for implementing a hypothetical strategy [3]. However, it is typically less plausible when implementing a treatment policy strategy for an ICE such as study treatment discontinuation. For this latter setting as well as for settings which include multiple ICEs handled with different strategies, methods based on imputation of missing data may be more suitable than MMRM [4]. For example, it may be feasible to impute missing data after study treatment discontinuation based on subjects from the same treatment group who also discontinued treatment but continued to be followed up [5, 6]. This is compatible with a more complex MAR assumption which also accounts for the subject’s discontinuation status. An alternative which is especially attractive if post-discontinuation data are sparse is to use reference-based imputation methods [7, 4, 3]. These methods formalize the idea to impute missing data in the intervention group based on data from a control or reference group. A typical approach to implement these imputation methods is to use multiple imputation based on Bayesian posterior draws of model parameters combined with Rubin’s rules to make inferences.

In this article, conditional mean imputation for MAR and reference-based imputation of missing data is introduced, justified, and explored in a clinical trial application and a simulation study. The proposed approach differs from the conventional multiple imputation methods in several ways. First, it estimates the parameters of the imputation model using maximum likelihood rather than via Bayesian posterior draws [8]. Second, it uses a single conditional mean imputation rather than random imputations. Third, inference is based on resampling methods, i.e. the jackknife or the bootstrap, and targets the frequentist repeated sampling variance. In contrast, the variance estimated by Rubin’s rules is approximately information-anchored, i.e. the proportion of information lost due to missing data under MAR is approximately preserved in reference-based imputation methods, but tends to over-estimate the frequentist variance for reference-based imputation methods [9, 10, 11, 12, 13].

Unlike the conventional approach, it avoids the specification of prior distributions and the complexity of Markov chain Monte Carlo (MCMC) sampling and leads to deterministic treatment effect estimates (i.e. they are free from Monte Carlo sampling error) and, in the case of the jackknife, also to deterministic confidence intervals and $p$-values.

In brief, the procedure is implemented in four steps. First, a MMRM imputation model is fitted to the observed outcome data. Second, conditional mean imputation of missing data is performed based on the parameter estimates from the MMRM imputation model and the chosen imputation method, e.g. reference-based imputation. Third, the completed data is analysed using a simple ANCOVA model. Finally, inference is performed based on resampling techniques. All four steps are described and justified below.

We applied the proposed methods to a publicly available example data set from an antidepressant clinical trial of an active drug versus placebo which was also used in Mallinckrodt et al. [23] and Tang [12]. The relevant endpoint is the Hamilton 17-item rating scale for depression (HAMD17) which was assessed at baseline and at weeks 1, 2, 4, and 6. Study drug discontinuation occurred in 24% (20/84) for the active drug and 26% (23/88) for placebo. All data after study drug discontinuation were missing and there was a single additional intermittent missing observation. The imputation model had the mean change from baseline in the HAMD17 score as the outcome, included the treatment group, the (categorical) visit, treatment-by-visit interactions, the baseline HAMD17 score, and baseline HAMD17-by-visit interactions as covariates, and assumed a common unstructured covariance matrix in both groups. The data set included one intermittent missing value which was imputed under a MAR assumption. For imputing missing data after study drug discontinuation, we explored MAR and reference-based imputation methods. The ANCOVA analysis model also adjusted for the baseline HAMD17 value. The results for conventional Bayesian multiple imputation combined with Rubin’s rules (with $M$=1,000 imputations) and conditional mean imputation using jackknife or bootstrap (with $B$=10,000 bootstrap samples) standard errors, respectively, are reported in Table 1.

As shown in Table 1, conventional Bayesian multiple imputation methods and conditional mean imputation result in nearly identical estimates of the treatment effect for all imputation approaches. However, for reference-based imputation methods, the resampling-based standard errors of the treatment effect from conditional mean imputations are considerably smaller than those from Rubin’s rules. The reported treatment effects estimates and standard errors are very similar to those obtained in Tang [12] using a different implementation. The treatment effect estimate and the standard error for conditional mean imputation were smallest for J2R and the bootstrap produced slightly smaller standard error estimates than the jackknife for all methods. The estimated computing time to obtain all four reported treatment effect estimates and corresponding standard errors with our R software implementation were 14.48 minutes for Bayesian multiple imputation ($M=1,000$), 5.16 minutes for conditional mean imputation plus jackknifing, and 312 minutes (5.2 hours) for conditional mean imputation plus bootstrapping ($B=10,000$, i.e. on average 1.87 seconds per bootstrap sample) on a single core of a low-power laptop with an Intel Core i5-8365U processor (1.60 GHz), 16 GB RAM, and Windows 10 Enterprise.

To further explore the statistical properties of the methods, we conducted a simulation study of a trial of an active drug versus placebo with 100 subjects per group and bi-monthly visits from baseline until 12 months. The other simulation parameters were chosen as follows:

• •

  The mean outcome trajectory in the placebo group (and in the active group under the null hypothesis) increased linearly from 50 at baseline to 60 at 12 months.

• •

  The mean outcome trajectory in the active group under the alternative hypothesis was identical to the placebo group up to month 4. From month 4 onward, the slope decreased by 50% to 5 points/year.

• •

  The covariance structure of the baseline and follow-up values in both groups was implied by a random intercept and slope model with a standard deviation of 5 for both the intercept and the slope, and a correlation of 0.25. In addition, an independent residual error with standard deviation 2.5 was added to each assessment. This implies marginal standard deviations ranging from 5.59 at baseline to 8.29 at the last visit in both groups.

• •

  The probability of study drug discontinuation after each visit was calculated according to a logistic model which also depended on the observed outcome at that visit. Specifically, a visit-wise discontinuation probability of 1.5% and 2.5% in the placebo and active group, respectively, was specified if the observed outcome was 50 or lower. For outcomes above 50, the odds of discontinuation further increase by 50% for each 10 point increase.

• •

  Study drug discontinuation had no effect on the mean trajectory in the placebo group. In the active group, subjects who discontinued followed the slope of the mean trajectory from the placebo group from that time point onward (CIR).

• •

  Study drop-out at the study drug discontinuation visit occurred with a probability of 75% leading to missing outcome data from that time point onward.

Under the alternative, the simulation parameters imply a true treatment effect of -2.59 (7.41 vs 10.00) under a treatment policy strategy which is the target of this simulation study. The treatment discontinuation probabilities were 34% and 24% in the active and placebo groups, respectively, under the null hypothesis and 31% and 24% under the alternative hypothesis.

Reference-based imputation methods of missing data were compared to MAR imputation. Observed post-discontinuation data was included in the analysis model but excluded from the imputation model for both reference-based and MAR-based imputation methods. Imputation and analysis models were as described for the antidepressant clinical trial application. Results for conventional Bayesian multiple imputation combined with Rubin’s rules were based on $M=100$ random imputation and bootstrap standard errors were calculated for $B$=1,000 bootstrap samples per simulation scenario.

Results of the simulations under the null hypothesis based on 100,000 simulated data sets are presented in Table 2. This choice of the number of simulations provides a Monte Carlo standard error for type I error estimates of approximately 0.07%. As expected, the mean treatment effect estimate was essentially 0 for all imputation methods and scenarios. For Bayesian multiple imputation, the observed standard deviation of treatment effect estimates across simulation runs ($SD$ $\hat{\theta}$) was substantially smaller than the mean treatment effect standard errors from Rubin’s rules (mean $\hat{se}_{\hat{\theta}}$) for reference-based imputation methods. Consequently, the simulated type I error ranged from 0.92%-2.54% and was substantially below the nominal type I error of 5%. For conditional mean imputation approaches, there was close agreement between the repeated sampling standard deviation $SD$ $\hat{\theta}$ and the average standard error. The simulated type I error for the jackknife ranged from 4.84% to 4.96% which was slightly below but very close to the targeted type I error of 5%. Inference based on bootstrap standard errors was associated with a small type I error inflation with simulated type I error ranging from 5.10% to 5.26%.

It is sometimes reported that MMRM models such as our imputation model using an unstructured covariance matrix result in convergence problems [2]. In this case, a simplified correlation structure, e.g. a Toeplitz or autoregressive structure, may be applied. In our simulations with a relatively small sample size and a substantial proportion of subjects with treatment discontinuations and missing data, the imputation model failed to converge in only 8/100,000 (0.008%) simulated data sets and these were excluded from our summaries in Table 1. Jackknifing did not result in any additional convergence problems and convergence problems occurred only for 1/((100,000-8)*1,000) bootstrap samples. In case convergence problems occurred for a bootstrap sample, we replaced it by a new bootstrap sample which seems unproblematic due to the minimal number of observed convergence problems.

Results of the simulations under the alternative hypothesis based on 10,000 simulated data sets (providing a Monte Carlo standard error for power estimates of <0.5%) are reported in Table 3. The mean treatment effect estimate under a CIR assumption was identical to the simulation truth of -2.59 for all three methods. MAR imputation gave inflated treatment effect estimates whereas CR and J2R were conservative. Rubin’s rules overestimated the frequentist standard error for all reference-based methods leading to a substantial reduction in the estimated power compared to methods based on conditional mean imputation plus re-sampling. The frequentist standard error was lowest for J2R. Surprisingly, this led to a larger simulated power for J2R than for CIR or CR. Inference based on bootstrapping led to small power gains of 0.3-0.4% across the four scenarios compared to the jackknife.

We introduced and justified conditional mean imputation based on the REML estimates of the parameters of the imputation model for MAR and reference-based imputation methods with longitudinal multivariate normal outcomes. For inference, we explored the use of the jackknife and the bootstrap. The data example and the simulations illustrated that this leads to nearly identical treatment effect estimates compared to established methods based on multiple imputation using Bayesian posterior draws form the imputation model parameters and Rubin’s rules. The resampling methods produced standard errors which were consistent with the observed repeated sampling variation. The jackknife demonstrated exact protection of the type I error in simulations with a relatively low sample size ($n=100$ per group) and a substantial amount of missing data (>25% of subjects with treatment discontinuations) whereas the bootstrap based on a normal approximation showed a slightly increased type I error rate compared to the nominal value (up to 5.26% at a nominal 5% significance level). The method based on Bayesian posterior draws and Rubin’s rules had inflated standard errors relative to the true repeated sampling variance, very conservative observed type I error rates, and a substantially decreased statistical power for reference-based imputation methods.

Based on our simulations, we recommend the jackknife over the bootstrap for inference because it performed better in our simulation study, is typically much faster to compute than the bootstrap, and leads to deterministic inference. In both the data example and the simulation study, jackknife standard errors tended to be slightly larger than bootstrap standard errors. This is in line with theoretical investigations which suggest that the jackknife variance tends to be biased upwards for finite sample sizes [24] whereas bootstrap variances are more susceptible to downwards bias [20, Section 3.4.3]. As a consequence, the jackknife demonstrated exact type I error protection whereas the bootstrap led to slight type I error inflation. While we haven’t formally investigated this, it could be that bootstrap methods which use more accurate methods for inference than the normal approximation have an advantage over the jackknife for very small sample sizes or other settings where the normal approximation to the test statistics is in doubt.

There are several differences between the proposed conditional mean imputation approach combined with the jackknife and conventional multiple imputation approaches based on Bayesian posterior draws combined with Rubin’s rules. First, our approach leads to deterministic treatment effect estimates, standard errors, and inference. This is particularly important in a regulatory setting where it is critical to ascertain whether a calculated p-value which is close to the critical boundary of 5% is truly below or above that threshold rather than being uncertain about this because of Monte Carlo error. Second, our approach does not rely on the specification of prior distributions, other tuning parameters (such as the number of imputations), or MCMC sampling. Third, previous publications on reference-based imputation methods have implicitly assumed that the corresponding post-ICE data is missing for all subjects whereas our implementation also covers the setting where post-ICE data is available for some subjects but requires imputation using reference-based methods for others. Fourth, our approach provides frequentist consistent estimates of the standard error for reference-based imputation methods whereas it is well known that Rubin’s rules do not provide frequentist consistent estimates of the standard error for reference-based imputation methods [9, 11, 12, 10, 13]. Intuitively, this occurs because reference-based imputation methods borrow information from the reference group for imputations in the intervention group. This creates a positive correlation between outcomes in the two treatment groups and leads to a reduction in the frequentist variance of the resulting treatment effect contrast which is not captured by Rubin’s variance estimator. Formally, the discrepancy is due to uncongeniality between the imputation and analysis models in reference-based imputation methods [25, 13]. Cro et al. [10] argued that Rubin’s variance (which cannot be obtained from our approach) is nevertheless valid for reference-based imputation methods because it is approximately information-anchored, i.e. the proportion of information lost due to missing data under MAR is approximately preserved in reference-based imputation methods. In contrast, several other authors have implicitly or explicitly favored inference that is correct from a frequentist repeated-sampling perspective [9, 11, 12, 13]. To us, information anchoring is a sensible concept for sensitivity analyses, whereas for a primary analyses, we feel that that it is more important to adhere to the principles of frequentist inference.

The proposed approach requires resampling for inference which is generally applicable but also computationally intensive. However, the time-consuming resampling procedure can easily be parallelized and computing time reduced if computations are distributed over multiple cores. A further limitation is that the validity of our approach has only been demonstrated for cases where parameter estimates from the analysis model are a linear transformation of the outcome such as general linear models including ANCOVA. This hinders the generalizability of the approach to binary or count outcomes. In addition, we have not covered the case of missing baseline covariates.

Under a standard MAR-based missing data assumption (accounting for baseline covariates, the treatment group, and observed outcomes), direct likelihood approaches such as MMRM are the most established and computationally efficient approach for treatment effect estimation. However, our approach and other approaches based on missing data imputation still play an important role for implementing more complex MAR assumptions. As an example, Noci et al. [26] implemented a treatment policy strategy by including interaction terms between the treatment group and two time-varying covariates in the imputation model to account for the impact of treatment discontinuations and initiations of symptomatic treatment, respectively, on outcomes. Completed data sets based on this more complex imputation model were then analyzed using a standard ANCOVA analysis model. It may be possible to derive an equivalent treatment effect estimator (and a corresponding standard error) directly from an MMRM model with time-varying covariates but this would be conceptually less straightforward and require additional algebraic manipulations. Providing formal justification and reliable implementations of such alternative MMRM-based approaches is an area of ongoing research.

We believe that methods based on missing data imputation are the most flexible and established approach to implement reference-based missing data assumptions but we acknowledge that related approaches based on MMRM model fitting have also been proposed by Mehrotra et al. [27] and Liu and Pang [11]. Instead of using subject‐level imputation, Mehrotra et al. [27] propose a simpler approach for the primary analysis in which the entire missing mean for the dropouts from the intervention group is explicitly replaced with the MAR‐based estimated overall control group mean. The resulting treatment effect estimator is the difference between the mean outcome of completers from the intervention group and the MAR‐based estimated mean outcome from the full control group, which is then multiplied by the proportion of completers from the intervention group. The mean outcomes in these two patient groups are obtained by separate MMRM models evaluated at the observed overall mean(s) of the baseline covariate(s) based on all randomized patients in the trial. Of note, Mehrotra et al. [27] compare two non-randomized patient groups (completers from the intervention group versus all subjects from the control group) and observed outcomes from subjects in the intervention group who drop out have no impact on the resulting treatment effect estimator. In case the probability of dropout depends only on baseline covariates but not on observed outcomes, this method leads to similar estimates as methods based on a jump-to-reference missing data assumption. However, in many clinical settings, it is plausible that subjects with worse observed outcome values are more likely to drop out. Consequently, completers in the intervention group would on average have better outcomes than the full intervention group. Under the null hypothesis of equal mean outcomes in the intervention and control group, the approach by Mehrotra et al. [27] would therefore target a treatment effect different from zero which could also lead to type I error inflation. In contrast, none of the reference-based imputation methods included in our simulation study showed evidence of bias under the null hypothesis or type I error inflation even though the probability of study drug discontinuation depended on observed outcomes. Liu and Pang [11] derived treatment effect estimators compatible with a reference-based imputation assumption based on parameter estimates from a standard MMRM fit to the data. However, it would be challenging to extend this approach to more general settings (e.g. allowing for different covariance structures in the two treatment groups or accounting for the inclusion of post-ICE data as proposed here).

Analyses of data with missing observations generally rely on unverifiable missing data assumptions and the assumptions for reference-based imputation methods are relatively strong [16, 13]. Informally, they assume that treatment discontinuations in the control or placebo group do not affect outcomes, i.e. that imputations under a standard MAR assumption are valid, and that subjects in the intervention group who discontinue treatment subsequently immediately follow the mean trajectories implied by the chosen reference-based imputation assumption. However, these assumptions may not be plausible if subjects who discontinue randomized treatment are likely to also initiate rescue medication or if the effects of discontinuing the intervention treatment are not immediate. Therefore, these assumptions need to be clinically justified as appropriate or at least conservative for the considered disease area and the anticipated mechanism of action of the intervention treatment. Moreover, if sufficient post-discontinuation data is available, it is possible to check whether reference-based imputations of missing post-discontinuation data are consistent with the observed post-discontinuation data. For example, one could graphically examine whether there are any systematic deviations between observed outcomes and corresponding conditional mean imputations. Furthermore, appropriate sensitivity analyses should always be planned for missing data methods. While we have not explicitly discussed it, $\delta$-based imputation approaches which adjust imputed data sets prior to the analysis for sensitivity analysis purposes could also be incorporated in our framework in a straightforward manner [4].

In summary, we believe that conditional mean imputation for MAR and reference-based imputation of missing data combined with the jackknife is a useful complementary approach to Bayesian multiple imputation for longitudinal multivariate normal outcomes. It provides replicable parameter estimates and standard errors and correct frequentist inference.

All presented methods have been implemented in an R package ‘rbmi‘ which which is available from CRAN at https://cran.r-project.org/web/packages/rbmi/.

The data analysis of the antidepressant clinical trial is described in the ‘Quick start‘ vignette of the package. All simulation code has been uploaded to github and can be found at https://github.com/nociale/Simulation_condmean_manuscript.

The authors thank Rachid Abbas for helpful remarks on a draft version of the manuscript.

Bartlett’s contribution was supported by the UK Medical Research Council (Grant MR/T023953/1).

Bartlett’s institution has received consultancy fees for the author’s advice on statistical methodology from AstraZeneca, Bayer, Novartis, and Roche. Bartlett has received consultancy and course fees from Bayer and Roche.