rbmi: Statistical Specifications

2.2.1 Missing data prior to ICEs

Missing data may occur in subjects without an ICE or prior to the occurrence of an ICE. As such missing outcomes are not associated with an ICE, it is often plausible to impute them under a missing-at-random (MAR) assumption using a standard MMRM imputation model of the longitudinal outcomes. Informally, MAR occurs if the missing data can be fully accounted for by the baseline variables included in the model and the observed longitudinal outcomes, and if the model is correctly specified.

2.2.2 Implementation of the hypothetical strategy

The MAR imputation model described above is often also a good starting point for imputing data after an ICE handled using a hypothetical strategy (Mallinckrodt et al. (2020)). Informally, this assumes that unobserved values after the ICE would have been similar to the observed data from subjects who did not have the ICE and remained under follow-up. However, in some situations, it may be more reasonable to assume that missingness is “informative” and indicates a systematically better or worse outcome than in observed subjects. In such situations, MNAR imputation with a $\delta$-adjustment could be explored as a sensitivity analysis. $\delta$-adjustments add a fixed or random quantity to the imputations in order to make the imputed outcomes systematically worse or better than those observed as described in Cro et al. (2020). In rbmi only fixed $\delta$-adjustments are implemented.

2.2.3 Implementation of the treatment policy strategy

Ideally, data collection continues after an ICE handled with a treatment policy strategy and no missing data arises. Indeed, such post-ICE data are increasingly systematically collected in RCTs. However, despite best efforts, missing data after an ICE such as study treatment discontinuation may still occur because the subject drops out from the study after discontinuation. It is difficult to give definite recommendations regarding the implementation of the treatment policy strategy in the presence of missing data at this stage because the optimal method is highly context dependent and a topic of ongoing statistical research.

For ICEs which are thought to have a negligible effect on efficacy outcomes, standard MAR-based imputation which ignores whether an outcome is observed pre- or post-ICE may be appropriate. In contrast, an ICE such as treatment discontinuation may be expected to have a more substantial impact on efficacy outcomes. In such settings, the MAR assumption may still be plausible after conditioning on the subject’s time-varying treatment status (Guizzaro et al. (2021)). In this case, one option is to impute missing post-discontinuation data based on subjects who also discontinued treatment but continued to be followed up. Another option which may require somewhat less post-discontinuation data is to include all subjects in the imputation procedure but to model post-discontinuation data by using time-varying treatment status indicators (Guizzaro et al. (2021), Polverejan and Dragalin (2020), Noci et al. (2023), Drury et al. (2024), Bell et al. (2024)). In this approach, post-ICE outcomes are included in every step of the analysis, including in the fitting of the imputation model. It assumes that ICEs may impact post-ICE outcomes but that otherwise missingness is non-informative. The approach also assumes that the time-varying covariates do not contain missing values, deviations in outcomes after the ICE are correctly modeled by these time-varying covariates, and that sufficient post-ICE data are available to inform the regression coefficients of the time-varying covariates. The resulting imputation models are called “retrieved dropout models” in the statistical literature. These models tend to have less bias than alternative analysis approaches based on imputation under a basic MAR assumption or a reference-based missing data assumption. However, retrieved dropout models have been associated with inflated standard errors of associated treatment effect estimators which has a detrimental effect on study power. In particular, it has been observed that once the post-ICE observation percentages falls below 50%, the power loss can be quite dramatic (Bell et al. 2024). We illustrate the implementation of retrieved dropout models in the vignette “Implementation of retrieved-dropout models using rbmi” (vignette(topic = "retrieved_dropout", package = "rbmi")).

In some trial settings, only few subjects discontinue the randomized treatment. In other settings, treatment discontinuation rates are higher but it is difficult to retain subjects in the trial after treatment discontinuation leading to sparse data collection after treatment discontinuation. In both settings, the amount of available data after treatment discontinuation may be insufficient to inform an imputation model which explicitly models post-discontinuation data. Depending on the disease area and the anticipated mechanism of action of the intervention, it may be plausible to assume that subjects in the intervention group behave similarly to subjects in the control group after the ICE treatment discontinuation. In this case, reference-based imputation methods are an option (Mallinckrodt et al. (2020)). Reference-based imputation methods formalize the idea to impute missing data in the intervention group based on data from a control or reference group. For a general description and review of reference-based imputation methods, we refer to Carpenter, Roger, and Kenward (2013), Cro et al. (2020), I. White, Royes, and Best (2020) and Wolbers et al. (2022). For a technical description of the implemented statistical methodology for reference-based imputation, we refer to section 3 (in particular section 3.4).

2.2.4 Implementation of the composite strategy

The composite strategy is typically applied to binary or time-to-event outcomes but it can also be used for continuous outcomes by ascribing a suitably unfavorable value to patients who experience ICEs for which a composite strategy has been defined. One possibility to implement this is to use MI with a $\delta$-adjustment for post-ICE data as described in Darken et al. (2020).

3.1.1 Conventional MI

Conventional MI approaches include the following steps:

1. Base imputation model fitting step (Section 3.3)

• Fit a Bayesian multivariate normal mixed model for repeated measures (MMRM) to the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired (Section 3.3.3). Draw $M$ posterior samples of the estimated parameters (regression coefficients and covariance matrices) from this model.

• Alternatively, $M$ approximate posterior draws from the posterior distribution can be sampled by repeatedly applying conventional restricted maximum-likelihood (REML) parameter estimation of the MMRM model to nonparametric bootstrap samples from the original dataset (Section 3.3.4).

2. Imputation step (Section 3.4)

• Take a single sample $m$ ($m\in 1,\dots ,M)$ from the posterior distribution of the imputation model parameters.

• For each subject, use the sampled parameters and the defined imputation strategy to determine the mean and covariance matrix describing the subject’s marginal outcome distribution for all longitudinal outcome assessments (i.e. observed and missing outcomes).

• For each subjects, construct the conditional multivariate normal distribution of their missing outcomes given their observed outcomes (including observed outcomes after ICEs for which a reference-based assumption is desired).

• For each subject, draw a single sample from this conditional distribution to impute their missing outcomes leading to a complete imputed dataset.

• For sensitivity analyses, a pre-defined $\delta$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

3. Analysis step (Section 3.6)

• Analyze the imputed dataset using an analysis model (e.g. ANCOVA) resulting in a point estimate and a standard error (with corresponding degrees of freedom) of the treatment effect.

4. Pooling step for inference (Section 3.7)

• Repeat steps 2. and 3. for each posterior sample $m$, resulting in $M$ complete datasets, $M$ point estimates of the treatment effect, and $M$ standard errors (with corresponding degrees of freedom). Pool the $M$ treatment effect estimates, standard errors, and degrees of freedom using the rules by Barnard and Rubin to obtain the final pooled treatment effect estimator, standard error, and degrees of freedom.

3.1.2 Conditional mean imputation

The conditional mean imputation approach includes the following steps:

1. Base imputation model fitting step (Section 3.3)

• Fit a conventional multivariate normal/MMRM model using restricted maximum likelihood (REML) to the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired (Section 3.3.2).

2. Imputation step (Section 3.4)

• For each subject, use the fitted parameters from step 1. to construct the conditional distribution of missing outcomes given observed outcomes (including observed outcomes after ICEs for which reference-based missing data imputation is desired) as described above.

• For each subject, impute their missing data deterministically by the mean of this conditional distribution leading to a complete imputed dataset.

• For sensitivity analyses, a pre-defined $\delta$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

3. Analysis step (Section 3.6)

• Apply an analysis model (e.g. ANCOVA) to the completed dataset resulting in a point estimate of the treatment effect.

4. Jackknife or bootstrap inference step (Section 3.8)

• Inference for the treatment effect estimate from 3. is based on re-sampling techniques. Both the jackknife and the bootstrap are supported. Importantly, these methods require repeating all steps of the imputation procedure (i.e. imputation, conditional mean imputation, and analysis steps) on each of the resampled datasets.

3.1.3 Bootstrapped MI

The bootstrapped MI approach includes the following steps:

1. Base imputation model fitting step (Section 3.3)

• Apply conventional restricted maximum-likelihood (REML) parameter estimation of the MMRM model to $B$ nonparametric bootstrap samples from the original dataset using the observed longitudinal outcomes after exclusion of data after ICEs for which reference-based missing data imputation is desired.

2. Imputation step (Section 3.4)

• Take a bootstrapped dataset $b$ ($b\in 1,\dots ,B)$ and its corresponding imputation model parameter estimates.

• For each subject (from the bootstrapped dataset), use the parameter estimates and the defined strategy for dealing with their ICEs to determine the mean and covariance matrix describing the subject’s marginal outcome distribution for all longitudinal outcome assessments (i.e. observed and missing outcomes).

• For each subjects (from the bootstrapped dataset), construct the conditional multivariate normal distribution of their missing outcomes given their observed outcomes (including observed outcomes after ICEs for which reference-based missing data imputation is desired).

• For each subject (from the bootstrapped dataset), draw $D$ samples from this conditional distributions to impute their missing outcomes leading to $D$ complete imputed dataset for bootstrap sample $b$.

• For sensitivity analyses, a pre-defined $\delta$-adjustment may be applied to the imputed data prior to the analysis step. (Section 3.5).

3. Analysis step (Section 3.6)

• Analyze each of the $B\times D$ imputed datasets using an analysis model (e.g. ANCOVA) resulting in $B\times D$ point estimates of the treatment effect.

4. Pooling step for inference (Section 3.9)

• Pool the $B\times D$ treatment effect estimates as described in von Hippel and Bartlett (2021) to obtain the final pooled treatment effect estimate, standard error, and degrees of freedom.

3.3.1 Included data and model specification

The purpose of the imputation model is to estimate (covariate-dependent) mean trajectories and covariance matrices for each group in the absence of ICEs handled using reference-based imputation methods. Conventionally, publications on reference-based imputation methods have implicitly assumed that the corresponding post-ICE data is missing for all subjects (Carpenter, Roger, and Kenward (2013)). We also allow the situation where post-ICE data is available for some subjects but needs to be imputed using reference-based methods for others. However, any observed data after ICEs for which reference-based imputation methods are specified is not compatible with the imputation model described below and they are therefore removed and considered as missing for the purpose of estimating the imputation model, and for this purpose only. For example, if a patient has an ICE addressed with a reference-based method but outcomes after the ICE are collected, these post-ICE outcomes will be excluded when fitting the base imputation model (but they will be included again in the following steps). That is, the base imputation model is fitted to $Y_{i!}^{\prime }$ and not to $Y_{i!}$. If we did not exclude these data, then the imputation model would mistakenly estimate mean trajectories based on a mixture of observed pre- and post-ICE data which are not relevant for reference-based imputations.

Observed post-ICE outcomes in the control or reference group are also excluded from the base imputation model if the user specifies a reference-based imputation strategy for such ICEs. This ensures that an ICE has the same impact on the data included in the imputation model regardless whether the ICE occurred in the control or the intervention group. On the other hand, imputation in the reference group is based on a MAR assumption even for reference-based imputation methods and it may be preferable in some settings to include post-ICE data from the control group in the base imputation model. This can be implemented by specifying a MAR strategy for the ICE in the control group and a reference-based strategy for the same ICE in the intervention group.

The base imputation model of the longitudinal outcomes $Y_i^{\prime }$ assumes that the mean structure is a linear function of covariates. Full flexibility for the specification of the linear predictor of the model is supported. At a minimum the covariates should include the treatment group, the (categorical) visit, and treatment-by-visit interactions. Typically, other covariates including the baseline outcome are also included. External time-varying covariates (e.g. calendar time of the visit) as well as internal time-varying (e.g. time-varying indicators of treatment discontinuation or initiation of rescue treatment) may in principle also be included if indicated (Guizzaro et al. (2021)). Missing covariate values are not allowed. This means that the values of time-varying covariates must be non-missing at every visit regardless of whether the outcome is measured or missing.

Denote the $J\times p$ design matrix for subject $i$ corresponding to the mean structure model by $X_i$ and the same matrix after removal of rows corresponding to missing outcomes in $Y_{i!}^{\prime }$ by $X_{i!}^{\prime }$. Here $p$ is the number of parameters in the mean structure of the model for the elements of $Y_{i!}^{\prime }$. The base imputation model for the observed outcomes is defined as: $$Y_{i!}^{\prime }=X_{i!}^{\prime }\beta +{ϵ}_{i!}\text{ with }{ϵ}_{i!}\sim N(0,{\mathrm{\Sigma }}_{i!!})$$ where $\beta$ is the vector of regression coefficients and ${\mathrm{\Sigma }}_{i!!}$ is a covariance matrix which is obtained from the complete-data $J\times J$-covariance matrix $\mathrm{\Sigma }$ by omitting rows and columns corresponding to missing outcome assessments for subject $i$.

Typically, a common unstructured covariance matrix for all subjects is assumed for $\mathrm{\Sigma }$ but separate covariate matrices per treatment group are also supported. Indeed, the implementation also supports the specification of separate covariate matrices according to an arbitrarily defined categorical variable which groups the subjects into disjoint subset. For example, this could be useful if different covariance matrices are suspected in different subject strata. Finally, for all imputation methods described below that do not rely on Bayesian model fitting through MCMC, there is further flexibility in the choice of the covariance structure, i.e. unstructured (default), heterogeneous Toeplitz, heterogeneous compound symmetry, and AR(1) covariance structures are supported.

3.3.2 Restricted maximum likelihood estimation (REML)

Frequentist parameter estimation for the base imputation is based on REML. The use of REML as an improved alternative to maximum likelihood (ML) for covariance parameter estimation was originally proposed by Patterson and Thompson (1971). Since then, it has become the default method for parameter estimation in linear mixed effects models. rbmi allows to choose between ML and REML methods to estimate the model parameters, with REML being the default option.

3.3.3 Bayesian model fitting

The Bayesian imputation model is fitted with the R package rstan (Stan Development Team (2020)). rstan is the R interface of Stan. Stan is a powerful and flexible statistical software developed by a dedicated team and implements Bayesian inference with state-of-the-art MCMC sampling procedures. The multivariate normal model with missing data specified in section 3.3.1 can be considered a generalization of the models described in the Stan user’s guide (see Stan Development Team (2020, sec. 3.5)).

The same prior distributions as in the SAS implementation of the “five macros” are used (Roger (2021)), i.e. an improper flat priors for the regression coefficients and a weakly informative inverse Wishart prior for the covariance matrix (or matrices). Specifically, let $S\in {\mathbb{R}}^{J\times J}$ be a symmetric positive definite matrix and $\nu \in (J-1,\mathrm{\infty })$. Then the symmetric positive definite matrix $x\in {\mathbb{R}}^{J\times J}$ has density: $$\text{InvWish}(x|\nu ,S)=\frac{1}{2^{\nu J/2}}\frac{1}{{\mathrm{\Gamma }}_J(\frac{\nu }{2})}|S{|}^{\nu /2}|x{|}^{-(\nu +J+1)/2}\text{exp}(-\frac{1}{2}\text{tr}(S{x}^{-1})).$$ For $\nu >J+1$ the mean is given by: $$E[x]=\frac{S}{\nu -J-1}.$$ We choose $S$ equal to the estimated covariance matrix from the frequentist REML fit and $\nu =J+2$ as these are the lowest degrees of freedom that guarantee a finite mean. Setting the degrees of freedom with such a low $\nu$ ensures that the prior has little impact on the posterior. Moreover, this choice allows to interpret the parameter $S$ as the mean of the prior distribution.

As in the “five macros”, the MCMC algorithm is initialized at the parameters from a frequentist REML fit (see section 3.3.2). As described above, we are using only weakly informative priors for the parameters. Therefore, the Markov chain is essentially starting from the targeted stationary posterior distribution and only a minimal amount of burn-in of the chain is required.

The initial values and other expert options for the MCMC sampling can be defined via the control argument to method_bayes, which is simplified with the corresponding control_bayes() function call. With the default initial values, i.e. when using init = "mmrm" in the control list, then for obtaining reproducible results an external set.seed() call is required before running draws(). In particular, please note that it is not sufficient to just set the seed option in the control list. If more than one chain is being run, the chains are currently initialized at "random" values as per rstan default. Note that for consistency with rstan, the burn-in can be customized via the warmup argument, and the burn-between can be customized via the thin argument. Independent of the number of chains used, the total number of returned samples is always n_samples. In the case that more than one chain is used, these samples are distributed across the multiple chains appropriately.

3.3.4 Approximate Bayesian posterior draws via the bootstrap

Several authors have suggested that a stabler way to get Bayesian posterior draws from the imputation model is to bootstrap the incomplete data and to calculate REML estimates for each bootstrap sample (Little and Rubin (2002), Efron (1994), Honaker and King (2010), von Hippel and Bartlett (2021)). This method is proper in that the REML estimates from the bootstrap samples are asymptotically equivalent to a sample from the posterior distribution and may provide additional robustness to model misspecification (Little and Rubin (2002, sec. 10.2.3, part 6), Honaker and King (2010)). In order to retain balance between treatment groups and stratification factors across bootstrap samples, the user is able to provide stratification variables for the bootstrap in the rbmi implementation.

3.4.1 Marginal imputation distribution for a subject - MAR case

For each subject $i$, the marginal distribution of the complete $J$-dimensional outcome vector from all assessment visits according to the imputation model is a multivariate normal distribution. Its mean ${\tilde{\mu }}_i$ is given by the predicted mean from the imputation model conditional on the subject’s baseline characteristics, group, and, optionally, time-varying covariates. Its covariance matrix ${\tilde{\mathrm{\Sigma }}}_i$ is given by the overall estimated covariance matrix or, if different covariance matrices are assumed for different groups, the covariance matrix corresponding to subject $i$’s group.

3.4.2 Marginal imputation distribution for a subject - reference-based imputation methods

For each subject $i$, we calculate the mean and covariance matrix of the complete $J$-dimensional outcome vector from all assessment visits as for the MAR case and denote them by ${\mu }_i$ and ${\mathrm{\Sigma }}_i$. For reference-based imputation methods, a corresponding reference group is also required for each group. Typically, the reference group for the intervention group will be the control group. The reference mean ${\mu }_{ref,i}$ is defined as the predicted mean from the imputation model conditional on the reference group (rather than the actual group subject $i$ belongs to) and the subject’s baseline characteristics. The reference covariance matrix ${\mathrm{\Sigma }}_{ref,i}$ is the overall estimated covariance matrix or, if different covariance matrices are assumed for different groups, the estimated covariance matrix corresponding to the reference group. In principle, time-varying covariates could also be included in reference-based imputation methods. However, this is only sensible for external time-varying covariates (e.g. calendar time of the visit) and not for internal time-varying covariates (e.g. treatment discontinuation) because the latter likely depend on the actual treatment group and it is typically not sensible to assume the same trajectory of the time-varying covariate for the reference group.

Based on these means and covariance matrices, the subject’s marginal imputation distribution for the reference-based imputation methods is then calculated as detailed in Carpenter, Roger, and Kenward (2013, sec. 4.3). Denote the mean and covariance matrix of this marginal imputation distribution by ${\tilde{\mu }}_i$ and ${\tilde{\mathrm{\Sigma }}}_i$. Recall that the subject’s first visit which is affected by the ICE is denoted by ${\tilde{t}}_i\in \{1,\dots ,J\}$ (and visit ${\tilde{t}}_i-1$ is the last visit unaffected by the ICE). The marginal distribution for the patient $i$ is then built according to the specific assumption for the data up to and post the ICE as follows:

1. Jump to reference (JR): the patient’s outcome distribution is normally distributed with the following mean: $${\tilde{\mu }}_i=({\mu }_i[1],\dots ,{\mu }_i[{\tilde{t}}_i-1],{\mu }_{ref,i}[{\tilde{t}}_i],\dots ,{\mu }_{ref,i}[J]{)}^T.$$ The covariance matrix is constructed as follows. First, we partition the covariance matrices ${\mathrm{\Sigma }}_i$ and ${\mathrm{\Sigma }}_{ref,i}$ in blocks according to the time of the ICE ${\tilde{t}}_i$: $${\mathrm{\Sigma }}_i=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{i,11}& {\mathrm{\Sigma }}_{i,12}\\ {\mathrm{\Sigma }}_{i,21}& {\mathrm{\Sigma }}_{i,22}\end{array}\right]$$ $${\mathrm{\Sigma }}_{ref,i}=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{ref,i,11}& {\mathrm{\Sigma }}_{ref,i,12}\\ {\mathrm{\Sigma }}_{ref,i,21}& {\mathrm{\Sigma }}_{ref,i,22}\end{array}\right].$$ We want the covariance matrix ${\tilde{\mathrm{\Sigma }}}_i$ to match ${\mathrm{\Sigma }}_i$ for the pre-deviation measurements, and ${\mathrm{\Sigma }}_{ref,i}$ for the conditional components for the post-deviation given the pre-deviation measurements. The solution is derived in Carpenter, Roger, and Kenward (2013, sec. 4.3) and is given by: $$\begin{array}{c}{\tilde{\mathrm{\Sigma }}}_{i,11}={\mathrm{\Sigma }}_{i,11}\\ {\tilde{\mathrm{\Sigma }}}_{i,21}={\mathrm{\Sigma }}_{ref,i,21}{\mathrm{\Sigma }}_{ref,i,11}^{-1}{\mathrm{\Sigma }}_{i,11}\\ {\tilde{\mathrm{\Sigma }}}_{i,22}={\mathrm{\Sigma }}_{ref,i,22}-{\mathrm{\Sigma }}_{ref,i,21}{\mathrm{\Sigma }}_{ref,i,11}^{-1}({\mathrm{\Sigma }}_{ref,i,11}-{\mathrm{\Sigma }}_{i,11}){\mathrm{\Sigma }}_{ref,i,11}^{-1}{\mathrm{\Sigma }}_{ref,i,12}.\end{array}$$

2. Copy increments in reference (CIR): the patient’s outcome distribution is normally distributed with the following mean: $$\begin{array}{rl}{\tilde{\mu }}_i=& ({\mu }_i[1],\dots ,{\mu }_i[{\tilde{t}}_i-1],{\mu }_i[{\tilde{t}}_i-1]+({\mu }_{ref,i}[{\tilde{t}}_i]-{\mu }_{ref,i}[{\tilde{t}}_i-1]),\dots ,\\ & {\mu }_i[{\tilde{t}}_i-1]+({\mu }_{ref,i}[J]-{\mu }_{ref,i}[{\tilde{t}}_i-1]){)}^T.\end{array}$$ The covariance matrix is derived as for the JR method.

3. Copy reference (CR): the patient’s outcome distribution is normally distributed with mean and covariance matrix taken from the reference group: $${\tilde{\mu }}_i={\mu }_{ref,i}$$ $${\tilde{\mathrm{\Sigma }}}_i={\mathrm{\Sigma }}_{ref,i}.$$

4. Last mean carried forward (LMCF): the patient’s outcome distribution is normally distributed with the following mean: $${\tilde{\mu }}_i=({\mu }_i[1],\dots ,{\mu }_i[{\tilde{t}}_i-1],{\mu }_i[{\tilde{t}}_i-1],\dots ,{\mu }_i[{\tilde{t}}_i-1]{)}^{\prime }$$ and covariance matrix: $${\tilde{\mathrm{\Sigma }}}_i={\mathrm{\Sigma }}_i.$$

3.4.3 Imputation of missing outcome data

The joint marginal multivariate normal imputation distribution of subject $i$’s observed and missing outcome data has mean ${\tilde{\mu }}_i$ and covariance matrix ${\tilde{\mathrm{\Sigma }}}_i$ as defined above. The actual imputation of the missing outcome data is obtained by conditioning this marginal distribution on the subject’s observed outcome data. Of note, this approach is valid regardless whether the subject has intermittent or terminal missing data.

The conditional distribution used for the imputation is again a multivariate normal distribution and explicit formulas for the conditional mean and covariance are readily available. For completeness, we report them here with the notation and terminology of our setting. The marginal distribution for the outcome of patient $i$ is $Y_i\sim N({\tilde{\mu }}_i,{\tilde{\mathrm{\Sigma }}}_i)$ and the outcome $Y_i$ can be decomposed in the observed ($Y_{i,!}$) and the unobserved ($Y_{i,?}$) components. Analogously the mean ${\tilde{\mu }}_i$ can be decomposed as $({\tilde{\mu }}_{i,!},{\tilde{\mu }}_{i,?})$ and the covariance ${\tilde{\mathrm{\Sigma }}}_i$ as: $${\tilde{\mathrm{\Sigma }}}_i=\left[\begin{array}{cc}{\tilde{\mathrm{\Sigma }}}_{i,!!}& {\tilde{\mathrm{\Sigma }}}_{i,!?}\\ {\tilde{\mathrm{\Sigma }}}_{i,?!}& {\tilde{\mathrm{\Sigma }}}_{i,??}\end{array}\right].$$ The conditional distribution of $Y_{i,?}$ conditional on $Y_{i,!}$ is then a multivariate normal distribution with expectation $$E(Y_{i,?}|Y_{i,!})={\tilde{\mu }}_{i,?}+{\tilde{\mathrm{\Sigma }}}_{i,?!}{\tilde{\mathrm{\Sigma }}}_{i,!!}^{-1}(Y_{i,!}-{\tilde{\mu }}_{i,!})$$ and covariance matrix $$Cov(Y_{i,?}|Y_{i,!})={\tilde{\mathrm{\Sigma }}}_{i,??}-{\tilde{\mathrm{\Sigma }}}_{i,?!}{\tilde{\mathrm{\Sigma }}}_{i,!!}^{-1}{\tilde{\mathrm{\Sigma }}}_{i,!?}.$$

Conventional random imputation consists in sampling from this conditional multivariate normal distribution. Conditional mean imputation imputes missing values with the deterministic conditional expectation $E(Y_{i,?}|Y_{i,!})$.

3.8.1 Point estimate of the treatment effect

The point estimator is obtained by applying the analysis model (Section 3.6) to a single conditional mean imputation of the missing data (see Section 3.4.3) based on the REML estimator of the parameters of the imputation model (see Section 3.3.2). We denote this treatment effect estimator by $\hat{\theta }$.

As demonstrated in Wolbers et al. (2022) (Section 2.4), this treatment effect estimator is valid if the analysis model is an ANCOVA model or, more generally, if the treatment effect estimator is a linear function of the imputed outcome vector. Indeed, if this is the case, then the estimator is identical to the pooled treatment effect across multiple random REML imputation with an infinite number of imputations and corresponds to a computationally efficient implementation of a proposal by von Hippel and Bartlett (2021). We expect that the conditional mean imputation method is also applicable to some other analysis models (e.g. for general MMRM analysis models) but this has not been formally justified.

3.8.2 Jackknife standard errors, confidence intervals (CI) and tests for the treatment effect

For a dataset containing $n$ subjects, the jackknife standard error depends on treatment effect estimates ${\hat{\theta }}_{(-b)}$ ($b=1,\dots ,n$) from samples of the original dataset which leave out the observation from subject $b$. As described previously, to obtain treatment effect estimates for leave-one-subject-out datasets, all steps of the imputation procedure (i.e. imputation, conditional mean imputation, and analysis steps) need to be repeated on this new dataset.

Then, the jackknife standard error is defined as $${\widehat{se}}_{jack}=[\frac{(n-1)}{n}\cdot \sum _{b=1}^n({\hat{\theta }}_{(-b)}-{\overline{\theta }}_{(.)}{)}^2{]}^{1/2}$$ where ${\overline{\theta }}_{(.)}$ denotes the mean of all jackknife estimates (Efron and Tibshirani (1994), chapter 10). The corresponding two-sided normal approximation $1-\alpha$ CI is defined as $\hat{\theta }\pm z^{1-\alpha /2}\cdot {\widehat{se}}_{jack}$ where $\hat{\theta }$ is the treatment effect estimate from the original dataset. Tests of the null hypothesis $H_0:\theta ={\theta }_0$ are then based on the $Z$-score $Z=(\hat{\theta }-{\theta }_0)/{\widehat{se}}_{jack}$ using a standard normal approximation.

A simulation study reported in Wolbers et al. (2022) demonstrated exact protection of the type I error for jackknife-based inference with a relatively low sample size (n = 100 per group) and a substantial amount of missing data (>25% of subjects with an ICE).

3.8.3 Bootstrap standard errors, confidence intervals (CI) and tests for the treatment effect

As an alternative to the jackknife, the bootstrap has also been implemented in rbmi (Efron and Tibshirani (1994), Davison and Hinkley (1997)).

Two different bootstrap methods are implemented in rbmi: Methods based on the bootstrap standard error and the normal approximation and percentile bootstrap methods. Denote the treatment effect estimates from $B$ bootstrap samples by ${\hat{\theta }}_b^{\ast }$ ($b=1,\dots ,B$). The bootstrap standard error ${\widehat{se}}_{boot}$ is defined as the empirical standard deviation of the bootstrapped treatment effect estimates. Confidence intervals and tests based on the bootstrap standard error can then be constructed in the same way as for the jackknife. Confidence intervals using the percentile bootstrap are based on empirical quantiles of the bootstrap distribution and corresponding statistical tests are implemented in rbmi via inversion of the confidence interval. Explicit formulas for bootstrap inference as implemented in the rbmi package and some considerations regarding the required number of bootstrap samples are included in the Appendix of Wolbers et al. (2022).

A simulation study reported in Wolbers et al. (2022) demonstrated a small inflation of the type I error rate for inference based on the bootstrap standard error (up to $5.3\mathrm{\%}$ for a nominal type I error rate of $5\mathrm{\%}$) for a sample size of n = 100 per group and a substantial amount of missing data (>25% of subjects with an ICE). Based on this simulations, we recommend the jackknife over the bootstrap for inference because it performed better in our simulation study and is typically much faster to compute than the bootstrap.

3.10.1 Treatment effect estimation

All approaches provide consistent treatment effect estimates for standard and reference-based imputation methods in case the analysis model of the completed datasets is a general linear model such as ANCOVA. Methods other than conditional mean imputation should also be valid for other analysis models. The validity of conditional mean imputation has only been formally demonstrated for analyses using the general linear model (Wolbers et al. (2022, sec. 2.4)) though it may also be applicable more widely (e.g. for general MMRM analysis models).

Treatment effects based on conditional mean imputation are deterministic. All other methods are affected by Monte Carlo sampling error and the precision of estimates depends on the number of imputations or bootstrap samples, respectively.

3.10.2 Standard errors of the treatment effect

All approaches for imputation under a MAR assumption provide consistent estimates of the frequentist standard error.

For reference-based imputation methods, the situation is more complicated and two different types of variance estimators have been proposed in the statistical literature (Bartlett (2023)). The first is the frequentist variance which describes the actual repeated sampling variability of the estimator. If the reference-based missing data assumption is correctly specified, then the resulting inference based on this variance is correct in the frequentist sense, i.e. hypothesis tests have asymptotically correct type I error control and confidence intervals have correct coverage probabilities under repeated sampling (Bartlett (2023), Wolbers et al. (2022)). Reference-based missing data assumptions are strong and borrow information from the reference arm for imputation in the active arm. As a consequence, the size of frequentist standard errors for treatment effects may decrease with increasing amounts of missing data. The second proposal is the so-called “information-anchored” variance which was originally proposed in the context of sensitivity analyses (Cro, Carpenter, and Kenward (2019)). This variance estimator is based on disentangling point estimation and variance estimation altogether. The information-anchoring principle described in Cro, Carpenter, and Kenward (2019) states that the relative increase in the variance of the treatment effect estimator under MAR imputation with increasing amounts of missing data should be preserved for reference-based imputation methods. The resulting information-anchored variance is typically very similar to the variance under MAR imputation and typically increases with increasing amounts of missing data. However, the information-anchored variance does not reflect the actual variability of the reference-based estimator under repeated sampling and the resulting inference is highly conservative resulting in a substantial power loss (Wolbers et al. (2022)). Moreover, to date, no Bayesian or frequentist framework has been developed under which the information-anchored variance provides correct inference for reference-based missingness assumptions, nor is it clear whether such a framework can even be developed.

Reference-based conditional mean imputation (method_condmean()) and bootstrapped likelihood-based multiple methods (method = method_bmlmi()) obtain standard errors via resampling and hence target the frequentist variance (Wolbers et al. (2022), von Hippel and Bartlett (2021)). For finite samples, simulations for a sample size of $n=100$ per group reported in Wolbers et al. (2022) demonstrated that conditional mean imputation combined with the jackknife (method_condmean(type = "jackknife")) provided exact protection of the type one error rate whereas the bootstrap (method_condmean(type = "bootstrap")) was associated with a small type I error inflation (between 5.1% to 5.3% for a nominal level of 5%). For reference-based conditional mean imputation, an alternative information-anchored variance can be obtained by following a proposal by Lu (2021). The basic idea of Lu (2021) is to obtain the information-anchored variance via a MAR imputation combined with a delta-adjustment where delta is selected in a data-driven way to match the reference-based estimator. For conditional mean imputation, the proposal by Lu (2021) can be implemented by choosing the delta-adjustment as the difference between the conditional mean imputation under the chosen reference-based assumption and MAR on the original dataset. An illustration of how the different variances can be obtained for conditional mean imputation in rbmi is provided in the vignette “Frequentist and information-anchored inference for reference-based conditional mean imputation” (vignette(topic = "CondMean_Inference", package = "rbmi")).

Reference-based Bayesian (or approximate Bayesian) multiple imputation methods combined with Rubin’s rules (method_bayes() and method_approxbayes()) target the information-anchored variance (Cro, Carpenter, and Kenward (2019)). A frequentist variance for these methods could in principle be obtained via bootstrap or jackknife re-sampling of the treatment effect estimates but this would be very computationally intensive and is not directly supported by rbmi.

Our view is that for primary analyses, accurate type I error control (which can be obtained by using the frequentist variance) is more important than adherence to the information anchoring principle which, to us, is not fully compatible with the strong reference-based assumptions. In any case, if reference-based imputation is used for the primary analysis, it is critical that the chosen reference-based assumption can be clinically justified, and that suitable sensitivity analyses are conducted to stress-test these assumptions.

Conditional mean imputation combined with the jackknife is the only method which leads to deterministic standard error estimates and, consequently, confidence intervals and $p$-values are also deterministic. This is particularly important in a regulatory setting where it is important 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.

3.10.3 Computational complexity

Bayesian MI methods rely on the specification of prior distributions and the usage of Markov chain Monte Carlo (MCMC) methods. All other methods based on multiple imputation or bootstrapping require no other tuning parameters than the specification of the number of imputations $M$ or bootstrap samples $B$ and rely on numerical optimization for fitting the MMRM imputation models via REML. Conditional mean imputation combined with the jackknife has no tuning parameters.

In our rbmi implementation, the fitting of the MMRM imputation model via REML is computationally most expensive. MCMC sampling using rstan (Stan Development Team (2020)) is typically relatively fast in our setting and requires only a small warmup phase and thinning of the chains. In addition, the number of random imputations for reliable inference using Rubin’s rules is often smaller than the number of resamples required for the jackknife or the bootstrap (see e.g. the discussions in I. R. White, Royston, and Wood (2011, sec. 7) for Bayesian MI and the Appendix of Wolbers et al. (2022) for the bootstrap). Thus, for many applications, we expect that conventional MI based on Bayesian posterior draws will be fastest, followed by conventional MI using approximate Bayesian posterior draws and conditional mean imputation combined with the jackknife. Conditional mean imputation combined with the bootstrap and bootstrapped MI methods will typically be most computationally demanding. Of note, all implemented methods are conceptually straightforward to parallelise and some parallelisation support is provided by rbmi.