A causal modelling framework for reference-based imputation and tipping point analysis

TL;DR

This paper introduces a causal modeling framework for reference-based imputation in clinical trials, providing formal justification, flexibility for sensitivity analysis, and application to real trial data.

Contribution

It offers a formal causal model that justifies and generalizes existing reference-based imputation methods, enabling transparent sensitivity analysis.

Findings

Simulation studies demonstrate the model's validity.

The framework allows flexible sensitivity analysis.

Application to clinical trial data illustrates practical utility.

Abstract

We consider estimating the "de facto" or effectiveness estimand in a randomised placebo-controlled or standard-of-care-controlled drug trial with quantitative outcome, where participants who discontinue an investigational treatment are not followed up thereafter. Carpenter et al (2013) proposed reference-based imputation methods which use a reference arm to inform the distribution of post-discontinuation outcomes and hence to inform an imputation model. However, the reference-based imputation methods were not formally justified. We present a causal model which makes an explicit assumption in a potential outcomes framework about the maintained causal effect of treatment after discontinuation. We show that the "jump to reference", "copy reference" and "copy increments in reference" reference-based imputation methods, with the control arm as the reference arm, are special cases of the causal model with specific assumptions about the causal treatment effect. Results from simulation studies are presented. We also show that the causal model provides a flexible and transparent framework for a tipping point sensitivity analysis in which we vary the assumptions made about the causal effect of discontinued treatment. We illustrate the approach with data from two longitudinal clinical trials.

Tables

Table 1. Table 1: Notation.

Notation	Explanation
$Z$	Randomised treatment arm: $Z=a$ for the active treatment and $Z=c$ for the control treatment.
$D_a$ or $D_c$	Participant’s last observation time prior to stopping treatment in arm $a$ or $c$.
$Y_t$	Response at time $t$, where $t=0,\dots ,t_{max}$.
$𝒀=({𝒀}_{\le t},{𝒀}_{>t})$	The vector of all outcomes from $t=0$ to $t=t_{max}$, with subvectors for outcomes until time $t$, and for outcomes from time $t+1$.
$Y_t\left(s\right)$	Potential outcome at time $t$ if $s$ periods of the active treatment were received, where $0\le s\le t_{max}$.
$𝒀\left(s\right)=({𝒀}_{\le t}\left(s\right),{𝒀}_{>t}\left(s\right))$	The vector of all potential outcomes if $s$ periods of the active treatment were received, with subvectors for outcomes until time $t$, and for outcomes from time $t+1$.
$𝝁\left(s\right)$	Vector of overall means of counterfactuals if participants would receive only $s$ periods of the active treatment.
$𝚺\left(s\right)$	Variance-covariance matrix of $\mathrm{var}\left(𝒀\left(s\right)\right)$ for $s=0,t_{max}$.
${𝜷}_t\left(s\right)$	Regression coefficient of ${𝒀}_{>t}\left(s\right)$ on ${𝒀}_{\le t}\left(s\right)$ for $s=0,t_{max}$.
${𝚺}_{>t>t|\le t}\left(s\right)$	Residual variance of ${𝒀}_{>t}\left(s\right)$ given ${𝒀}_{\le t}\left(s\right)$ for $s=0,t_{max}$.

Table 2. Table 2: Mean of the imputation distribution of 𝒀 > t ​ ( t ) subscript 𝒀 absent 𝑡 𝑡 \boldsymbol{Y}_{>t}(t) for t < t m ​ a ​ x 𝑡 subscript 𝑡 𝑚 𝑎 𝑥 t<t_{max} given randomisation Z = a 𝑍 𝑎 Z=a , past 𝒀 ≤ t subscript 𝒀 absent 𝑡 \boldsymbol{Y}_{\leq{t}} and treatment discontinuation time D a = t subscript 𝐷 𝑎 𝑡 D_{a}=t , under various reference-based imputation methods (Carpenter et al. , 2013 ) and the proposed causal model.

Method	Mean of the imputation distribution
Reference-based imputation methods
LMCF	${𝜷}_t\left(t_{max}\right)\left\{{𝒀}_{\le t}-{𝝁}_{\le t}\left(t\right)\right\}$	$+{𝑪}_t{𝝁}_{\le t}\left(t\right)$
J2R*	${𝜷}_t\left(0\right)\left\{{𝒀}_{\le t}-{𝝁}_{\le t}\left(t\right)\right\}$		$+{𝝁}_{>t}\left(0\right)$
CIR*	${𝜷}_t\left(0\right)\left\{{𝒀}_{\le t}-{𝝁}_{\le t}\left(t\right)\right\}$	$+{𝑪}_t\left\{{𝝁}_{\le t}\left(t\right)-{𝝁}_{\le t}\left(0\right)\right\}$	$+{𝝁}_{>t}\left(0\right)$
CR*	${𝜷}_t\left(0\right)\left\{{𝒀}_{\le t}-{𝝁}_{\le t}\left(t\right)\right\}$	$+{𝜷}_t\left(0\right)\left\{{𝝁}_{\le t}\left(t\right)-{𝝁}_{\le t}\left(0\right)\right\}$	$+{𝝁}_{>t}\left(0\right)$
Causal model	${𝜷}_t\left(t\right)\left\{{𝒀}_{\le t}-{𝝁}_{\le t}\left(t\right)\right\}$	$+{𝑲}_t\left\{{𝝁}_{\le t}\left(t\right)-{𝝁}_{\le t}\left(0\right)\right\}$	$+{𝝁}_{>t}\left(0\right)$
*with control arm as the reference and variance-covariance from the control arm

Table 3. Table 3: Simulation study: estimates of treatment effect at the final visit using complete data, causal model imputation and RBI imputation. 𝜷 1 ​ ( 0 ) = ( 0 , 0.5 ) ′ subscript 𝜷 1 0 superscript 0 0.5 ′ \boldsymbol{\beta}_{1}(0)=(0,0.5)^{\prime} in all cases. 𝜷 1 ​ ( 2 ) ≠ 𝜷 1 ​ ( 0 ) subscript 𝜷 1 2 subscript 𝜷 1 0 \boldsymbol{\beta}_{1}(2)\neq\boldsymbol{\beta}_{1}(0) means 𝜷 1 ​ ( 2 ) = ( − 0.12 , 0.74 ) ′ subscript 𝜷 1 2 superscript 0.12 0.74 ′ \boldsymbol{\beta}_{1}(2)=(-0.12,0.74)^{\prime} .

	Data generating mechanisms for observed data
	MCAR		MAR
	${𝜷}_1\left(2\right)={𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)\ne {𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)={𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)\ne {𝜷}_1\left(0\right)$
A. Complete data generated with:
${𝜷}_1\left(1\right)={𝜷}_1\left(0\right)$
$k=0.00$	0.99	0.99	1.00	0.70
$k=0.50$	1.24	1.24	1.25	0.95
$k=0.74$	1.36	1.36	1.37	1.07
$k=1.00$	1.49	1.49	1.50	1.20
${𝜷}_1\left(1\right)={𝜷}_1\left(2\right)$
$k=0.00$	0.99	1.00	1.00	1.00
$k=0.50$	1.24	1.25	1.25	1.25
$k=0.74$	1.36	1.37	1.37	1.37
$k=1.00$	1.49	1.50	1.50	1.50
B. Causal model imputation with assumed ${\widetilde{𝜷}}_1\left(1\right)$ and $\tilde{k}$
${\widetilde{𝜷}}_1\left(1\right)={𝜷}_1\left(0\right)$
$\tilde{k}=0.00$	1.00	1.00	1.00	0.71
$\tilde{k}=0.50$	1.24	1.25	1.25	0.96
$\tilde{k}=0.74$	1.36	1.37	1.37	1.08
$\tilde{k}=1.00$	1.49	1.50	1.50	1.21
${\widetilde{𝜷}}_1\left(1\right)={𝜷}_1\left(2\right)$
$\tilde{k}=0.00$	1.00	1.00	1.00	1.00
$\tilde{k}=0.50$	1.24	1.25	1.25	1.25
$\tilde{k}=0.74$	1.36	1.37	1.37	1.37
$\tilde{k}=1.00$	1.49	1.50	1.50	1.50
C. RBI imputation with assumed variance-covariance and method
Variance-covariance matrix from control arm
J2R	1.00	1.00	1.00	0.71
CR	1.24	1.25	1.25	0.96
CIR	1.49	1.50	1.50	1.21
Variance-covariance matrix from active arm
J2R	1.00	1.00	1.00	1.00
CR	1.24	1.37	1.25	1.38
CIR	1.49	1.50	1.50	1.50
Note: Maximum Monte Carlo standard error

Table 4. Table 4: Simulation study: average standard error (empirical standard error) for the treatment difference at the final visit using complete data, causal model imputation and RBI imputation. 𝜷 1 ​ ( 0 ) subscript 𝜷 1 0 \boldsymbol{\beta}_{1}(0) and 𝜷 1 ​ ( 2 ) subscript 𝜷 1 2 \boldsymbol{\beta}_{1}(2) as in Table 3 .

	Data generating mechanisms for observed data
	MCAR		MAR
	${𝜷}_1\left(2\right)={𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)\ne {𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)={𝜷}_1\left(0\right)$	${𝜷}_1\left(2\right)\ne {𝜷}_1\left(0\right)$
A. Complete data generated with:
${𝜷}_1\left(1\right)={𝜷}_1\left(0\right)$
$k=0.00$	0.276 (0.273)	0.318 (0.311)	0.260 (0.255)	0.295 (0.288)
$k=0.50$	0.272 (0.269)	0.315 (0.308)	0.261 (0.256)	0.298 (0.291)
$k=0.74$	0.302 (0.206)	0.342 (0.247)	0.311 (0.219)	0.351 (0.258)
$k=1.00$	0.270 (0.266)	0.313 (0.306)	0.262 (0.257)	0.301 (0.295)
${𝜷}_1\left(1\right)={𝜷}_1\left(2\right)$
$k=0.00$	0.276 (0.273)	0.318 (0.316)	0.260 (0.255)	0.292 (0.288)
$k=0.50$	0.272 (0.269)	0.315 (0.312)	0.261 (0.256)	0.296 (0.291)
$k=0.74$	0.271 (0.268)	0.314 (0.311)	0.262 (0.256)	0.298 (0.294)
$k=1.00$	0.270 (0.266)	0.313 (0.310)	0.262 (0.257)	0.300 (0.296)
B. Causal model imputation with assumed ${\widetilde{𝜷}}_1\left(1\right)$ and $\tilde{k}$
${\widetilde{𝜷}}_1\left(1\right)={𝜷}_1\left(0\right)$
$\tilde{k}=0.00$	0.310 (0.168)	0.337 (0.189)	0.305 (0.171)	0.328 (0.181)
$\tilde{k}=0.50$	0.301 (0.190)	0.327 (0.226)	0.299 (0.189)	0.322 (0.214)
$\tilde{k}=0.74$	0.302 (0.206)	0.327 (0.249)	0.301 (0.206)	0.325 (0.237)
$\tilde{k}=1.00$	0.305 (0.226)	0.332 (0.277)	0.306 (0.227)	0.332 (0.267)
${\widetilde{𝜷}}_1\left(1\right)={𝜷}_1\left(2\right)$
$\tilde{k}=0.00$	0.310 (0.168)	0.359 (0.187)	0.315 (0.188)	0.359 (0.206)
$\tilde{k}=0.50$	0.301 (0.190)	0.344 (0.224)	0.310 (0.204)	0.350 (0.236)
$\tilde{k}=0.74$	0.302 (0.206)	0.342 (0.247)	0.311 (0.219)	0.351 (0.258)
$\tilde{k}=1.00$	0.305 (0.226)	0.345 (0.275)	0.316 (0.239)	0.356 (0.285)
C. RBI imputation with assumed variance-covariance and method
Variance-covariance matrix from control arm
J2R	0.310 (0.168)	0.337 (0.189)	0.305 (0.171)	0.328 (0.181)
CR	0.303 (0.192)	0.328 (0.229)	0.306 (0.200)	0.331 (0.226)
CIR	0.305 (0.226)	0.332 (0.277)	0.306 (0.227)	0.332 (0.267)
Variance-covariance matrix from active arm
J2R	0.310 (0.168)	0.359 (0.187)	0.315 (0.188)	0.359 (0.206)
CR	0.305 (0.192)	0.344 (0.249)	0.332 (0.236)	0.370 (0.283)
CIR	0.305 (0.226)	0.345 (0.275)	0.316 (0.239)	0.356 (0.285)
Note: Maximum Monte Carlo standard error

Table 5. Table 5: HAMD17 and pain score data: estimated treatment effect at the final visit using standard multiple imputation with 100 imputations, mixed model for repeated measures (MMRM) and RBI methods.

Estimand &	HAMD17				Pain score
Method	Estimate1	Std. error	p-value		Estimate2	Std. error	p-value
De jure
Standard MI	-2.62	0.99	0.01		-0.88	0.39	0.03
MMRM	-2.58	1.03	0.01		-0.88	0.39	0.03
De facto
RBI1: variance-covariance from the placebo arm
J2R	-2.01	1.01	0.05		-0.64	0.40	0.11
CR	-2.22	0.99	0.03		-0.75	0.39	0.06
CIR	-2.30	0.99	0.02		-0.77	0.39	0.05
RBI2: variance-covariance from the drug arm
J2R	-1.99	1.01	0.05		-0.60	0.39	0.13
CR	-2.20	0.99	0.03		-0.71	0.39	0.07
CIR	-2.28	0.99	0.02		-0.73	0.39	0.06

Equations

$$Y_t(s)=Y_t(t)\text{mbox}foralls>t.$$

$$p(D_c=t|Z=c,\mathbf{Y},D_c\ge t)=p(D_c=t|Z=c,{\mathbf{Y}}_{\le t},D_c\ge t)$$

$$p(D_a=t|Z=a,\mathbf{Y}(t_{max}),D_a\ge t)=p(D_a=t|Z=a,{\mathbf{Y}}_{\le t}(t_{max}),D_a\ge t)$$

$$E[{\mathbf{Y}}_{>t}(t)|Z=a,{\mathbf{Y}}_{\le t},D_a=t]$$

$$p(D_a=t|Z=a,\mathbf{Y}(s),D_a\ge t)=p(D_a=t|Z=a,{\mathbf{Y}}_{\le t}(s),D_a\ge t)$$

$$\begin{array}[b]{llll}{\mbox{E}}\left[\boldsymbol{Y}_{>t}(t)|Z=a,\boldsymbol{Y}_{\leq{t}}(t),D_{a}=t\right]&-&\boldsymbol{\mu}_{>t}(t)&\\ \hskip 20.00003pt={\mbox{E}}\left[\boldsymbol{Y}_{>t}(t)|Z=a,\boldsymbol{Y}_{\leq{t}}(t)\right]&-&\boldsymbol{\mu}_{>t}(t)&\mbox{ (by \ref{ass:actMAR2})}\\ \hskip 20.00003pt={\mbox{E}}\left[\boldsymbol{Y}_{>t}(t)|\boldsymbol{Y}_{\leq{t}}(t)\right]&-&\boldsymbol{\mu}_{>t}(t)&\mbox{ (by randomisation)}\\ \hskip 20.00003pt=\boldsymbol{\beta}_{t}(t)\left\{\boldsymbol{Y}_{\leq{t}}(t)-\boldsymbol{\mu}_{\leq{t}}(t)\right\}&&&\mbox{ (by \ref{ass:linear})}\end{array}$$

$${\mathbf{\mu }}_{>t}(t)-{\mathbf{\mu }}_{>t}(0)={\mathbf{K}}_t\left\{{\mathbf{\mu }}_{\le t}(t)-{\mathbf{\mu }}_{\le t}(0)\right\}$$

$$\text{mbox}E\left[Y_s(t)-Y_s(0)\right]$$

$$E[Y_t{}_{max}(D_a)-Y_t{}_{max}(0)]={\alpha }_t{}_{max}{\delta }_t{}_{max}+\sum _{t<t_{max}}{\alpha }_t{\mathbf{e}}_t^{\prime }{\mathbf{K}}_t{\mathbf{\delta }}_t$$

$$+\sum _{t<t_{max}}{\alpha }_t{\mathbf{e}}_t^{\prime }({\mathbf{\beta }}_t(t)-{\mathbf{\beta }}_t(t_{max}))\text{mbox}E\left[{\mathbf{Y}}_{\le t}(t|D_a=t)-{\mathbf{Y}}_{\le t}(t)\right]$$

$${\mathbf{K}}_t=k{\mathbf{C}}_t$$

$$E[{\mathbf{Y}}_{t_{max}}(D_a)-{\mathbf{Y}}_{t_{max}}(0)]={\delta }_{J2R}+\sum _{t<t_{max}}{\alpha }_t{\mathbf{e}}_t^{\prime }{\mathbf{K}}_t{\mathbf{\delta }}_t$$

$$E[{\mathbf{Y}}_{t_{max}}(D_a)-{\mathbf{Y}}_{t_{max}}(0)]$$

$$E[Y_t{}_{max}(D_a)-Y_t{}_{max}(0)]=\sum _{t=0}^{t_{max}}{\alpha }_t\left\{{\mu }_t{}_{max}(t|D_a=t)-{\mu }_t{}_{max}(0)\right\}$$

$${\mu }_t{}_{max}(t|D_a=t)-{\mu }_t{}_{max}(0)$$

$${\mathbf{\beta }}_t(t_{max})\left\{{\mathbf{\mu }}_{\le t}(t|D_a=t)-{\mathbf{\mu }}_{\le t}(t)\right\}$$

$$\sum _{t<t_{max}}{\alpha }_t{\mathbf{e}}_t^{\prime }{\mathbf{\beta }}_t(t_{max})\left\{{\mathbf{\mu }}_{\le t}(t|D_a=t)-{\mathbf{\mu }}_{\le t}(t)\right\}$$

$$+{\alpha }_t{}_{max}\left\{{\mu }_t{}_{max}(t_{max}|D_a=t_{max})-{\mu }_t{}_{max}(t_{max})\right\}=0$$

$${\alpha }_t{}_{max}{\mu }_t{}_{max}(t_{max}|D_a=t_{max})={\alpha }_t{}_{max}{\mu }_t{}_{max}(t_{max})$$

$$-\sum _{t<t_{max}}{\alpha }_t{\mathbf{e}}_t^{\prime }{\mathbf{\beta }}_t(t_{max})\left\{{\mathbf{\mu }}_{\le t}(t|D_a=t)-{\mathbf{\mu }}_{\le t}(t)\right\}.$$

$$E[Y_t{}_{max}(D)-Y_t{}_{max}(0)]$$

Full text

1 Introduction

Missing outcome data represent a major threat to the validity of randomised clinical trials (RCTs), and appropriate analysis methods have been much discussed. An influential report showed that different analysis methods may target different estimands (different measures of treatment effect) and argued that specification of the estimand is an important part of trial design and should inform trial analysis and reporting (National Research Council, 2010). Regulators have joined the call for estimands to be defined clearly (ICH steering committee, 2014).

Estimands have long been discussed in the causal inference literature. For example, Angrist et al. (1996) discussed all-or-nothing non-compliance in RCTs and gave conditions under which instrumental variable methods estimate the complier-average causal effect (CACE), that is, the mean effect of treatment in those who will receive treatment if and only if they are assigned to treatment. This differs from the intention-to-treat estimand which is the mean effect of allocation to treatment. In general, the estimand can be designed to take into account or not take into account potential confounding of the data by post-randomisation events such as non-compliance, discontinuation of intervention, treatment switching, or use of rescue medication.

We consider two types of estimand considered by National Research Council (2010): (E1) difference in outcome improvement at the planned endpoint if all participants had tolerated or adhered to trial protocol; (E2) difference in outcome improvement at the planned endpoint for all randomised participants. The former measures how treatment works in an ideal setting (efficacy), while the latter measures how treatment might work in practice (effectiveness). To encompass outcomes that measure harms of treatment, Carpenter et al. (2013) proposed the broader terms de jure and de facto estimand for (E1) and (E2) respectively.

Sometimes investigators continue to collect data after treatment discontinuation, which we call “off-treatment data”. The use of such data depends on the estimand (Permutt, 2016). For the estimation of a de jure estimand, off-treatment data for participants who discontinued randomised treatment could be used in a CACE analysis (Dunn et al., 2003). In practice, however, off-treatment data are typically either not collected or are excluded from the primary analysis, and the missing data are assumed to be missing at random (MAR): that is, it is assumed that participants who discontinued treatment would have benefited from continued treatment in the same way as those who remained on treatment. However, estimation of a de facto estimand ideally makes use of the off-treatment data, and guidance recommends its collection where possible (National Research Council, 2010). When all discontinuers are followed up and outcome data obtained, the de facto estimand can be directly estimated by comparing observed means (Little and Kang, 2015).

This paper considers estimation of a de facto estimand for a quantitative outcome in the absence of off-treatment data. This requires assumptions about how participants who have discontinued treatment benefit from their previous treatment. For example, participants who received only a partial period of treatment may not have benefited after treatment discontinuation in a similar manner to those who completed the treatment in the trial. We also focus on the situation in which rescue treatment (over and above the per protocol treatment regime for the control arm) for those who discontinue randomised treatment is not available, or interest is in the effect attributable to the initially randomised treatment without the confounding effects of rescue medications (where the latter corresponds to estimand 6 in Mallinckrodt et al. (2012)).

In previous work on this topic, Little and Yau (1996) presented a multiple imputation (MI) approach that can incorporate a variety of alternative assumptions about treatment dose after treatment discontinuation for the estimation of de facto estimands in RCTs. Carpenter et al. (2013) presented a generic algorithm for MI of post-discontinuation outcome data. They derived the conditional distribution of post-discontinuation outcome data given pre-discontinuation outcome data and time of discontinuation from a joint Normal distribution which could be formed with reference to other arms in several ways. For the “copy reference” (CR) method, the mean of the joint distribution is taken as the mean for the reference arm. For other methods, the mean of the joint distribution in a specific arm follows the mean of that arm until treatment discontinuation; after treatment discontinuation, it is formed either by directly using the means over time for the reference arm (the “jump to reference” (J2R) method); by using the mean in the specific arm at the time of treatment discontinuation (the “last mean carried forward” (LMCF) method); or by using the mean for the specific arm at the time of treatment discontinuation plus the mean increments over time in the reference arm (the “copy increments in reference” (CIR) method). The variance-covariance matrix of the joint distribution is set equal to the variance-covariance matrix from either the active (investigational) treatment arm or the reference arm. However, Carpenter et al. (2013) did not provide a theoretical justification of these “reference-based imputation” (RBI) methods.

Specification of estimands is clarified by using counterfactual outcomes – outcomes that have not or could not have been observed. Such counterfactuals are best described using potential outcomes notation (Little and Rubin, 2000; Angrist et al., 1996). Our aims in this paper are to provide mathematical derivation of the RBI methods, including their identifying assumptions, in a potential outcomes framework, and to propose a post-discontinuation imputation model using explicit causal assumptions to enable generalization of the RBI methods. We do this by assuming that after discontinuing their randomised treatment, participants receive treatment that is similar to that allocated to the control arm: thus in the case of reference-based imputation, we take the reference to be the control treatment, where this is typically either placebo or standard of care (SoC). Further, while we allow both a delayed response to the active treatment and a maintained response after discontinuing the active treatment, we implicitly assume that there is no delayed response to the control treatment: thus when a patient discontinues randomised treatment, we assume the effects of any treatments they switch to are similar to the effects they would have experienced had they received the control treatment from the start of the trial.

Since the imputation models in our proposal rely on untestable assumptions, sensitivity analyses are desirable for understanding the impact of the assumptions on inferences and conclusions from the primary analysis. Kenward et al. (2001) described a principled approach to sensitivity analyses which varies a sensitivity parameter that quantifies deviations from the missing data assumption. A tipping point sensitivity analysis (e.g. Yan et al. (2009); Liublinska and Rubin (2014)) extends this approach by varying the sensitivity parameter until the conclusion from the primary analysis is overturned. In this paper, we propose a tipping point sensitivity analysis using the causal model.

In Section 2, we investigate the underlying assumptions for RBI methods and present a causal model that requires an explicit assumption on post-discontinuation behaviour in terms of the average treatment effect which is maintained after treatment discontinuation. In Section 3, we demonstrate the equivalence of RBI and causal model estimates in a simulation study. In Section 4, we illustrate the proposed approach and demonstrate the tipping point analysis with two example data sets. We conclude with summary remarks in Section 5.

2 RBI methods and causal model

Our aim here is to propose a causal model which can be used to impute missing values, and to relate the mean imputations from the causal model to the mean imputations derived under RBI. We start by presenting notation and estimands, and then present the assumptions used in imputation, a comparison of the methods, and details of implementation.

2.1 Notation and simplifying assumptions

We consider a two-arm longitudinal RCT with $t_{max}$ scheduled observation times after randomisation. Let $Z$ denote the randomised treatment arm with $Z=a$ for the active treatment and $Z=c$ for the control treatment. Let $Y_{t}$ denote the outcome at the $t$th observation time, where $t=0,1,...,t_{max}$ and $t=0$ denotes the baseline observation. Each $Y_{t}$ may be observed or missing; if missing, it is the value of the actual outcome that is experienced by the participant despite being unobserved by the investigators. Let $\boldsymbol{Y}=(Y_{0},...,Y_{t}{}_{max})^{^{\prime}}$ denote the vector of all outcomes, with subvectors $\boldsymbol{Y}_{\leq{t}}=(Y_{0},...,Y_{t})^{^{\prime}}$ for outcomes until time $t$ (including the baseline), and $\boldsymbol{Y}_{>t}=(Y_{t+1},...,Y_{t}{}_{max})^{^{\prime}}$ for outcomes after time $t$. We assume for clarity that the only baseline covariate is the baseline value of the outcome, $Y_{0}$, but other baseline covariates can be accommodated by replacing $Y_{0}$ with a baseline covariate vector $\boldsymbol{Y}_{0}$.

We assume that individuals in both arms may discontinue their randomised treatment during the course of follow-up, but they do not restart treatment. We let $D_{a}$, $0\leq D_{a}\leq t_{max}$, denote the last observation time for an active arm participant prior to stopping treatment, and similarly for $D_{c}$ in the control arm. $D_{a}$ or $D_{c}=t_{max}$ denotes a participant who continues treatment throughout the trial. We assume that, by the design of the trial, there is no follow-up after treatment discontinuation, and that data may also be intermittently missing before treatment discontinuation.

We define the potential outcome $Y_{t}(s)$, for $0\leq s,t\leq t_{max}$, as the true outcome at time $t$ if the active treatment were received for $s$ periods and then control treatment were received. We assume that outcomes are not causally affected by future treatment, so

[TABLE]

In particular, (1) means that the baseline value is unaffected by treatment, $Y_{0}(s)=Y_{0}$ for all $s$. All counterfactual outcomes $Y_{t}(s)$ are assumed to be independent of $Z$ (by randomisation), but not of $D_{a}$ or $D_{c}$ (treatment discontinuation may be related to counterfactual outcomes). Let $\boldsymbol{Y}(s)=(Y_{0}(s),...,Y_{t}{}_{max}(s))^{^{\prime}}$ with subvectors $\boldsymbol{Y}_{\leq{t}}(s)=(Y_{0}(s),...,Y_{t}(s))^{^{\prime}}$ and $\boldsymbol{Y}_{>t}(s)=(Y_{t+1}(s),...,Y_{t}{}_{max}(s))^{^{\prime}}$. The overall means of counterfactuals are defined as $\boldsymbol{\mu}(s)={\mbox{E}}\left[\boldsymbol{Y}(s)\right]$, $\boldsymbol{\mu}_{\leq{t}}(s)={\mbox{E}}\left[\boldsymbol{Y}_{\leq{t}}(s)\right]$ and $\boldsymbol{\mu}_{>t}(s)={\mbox{E}}\left[\boldsymbol{Y}_{>t}(s)\right]$.

Finally, let the variance-covariance matrices be $\boldsymbol{\Sigma}(s)={\rm var}\left(\boldsymbol{Y}(s)\right)$ for $s=0$ or $t_{max}$, with submatrices $\boldsymbol{\Sigma}_{\leq{t}\leq{t}}(s)$, $\boldsymbol{\Sigma}_{>t\leq{t}}(s)$, $\boldsymbol{\Sigma}_{\leq{t}>t}(s)$, $\boldsymbol{\Sigma}_{>t>t}(s)$. For the linear regression of $\boldsymbol{Y}_{>t}(s)$ on $\boldsymbol{Y}_{\leq{t}}(s)$, we write the regression coefficients as $\boldsymbol{\beta}_{t}(s)=\boldsymbol{\Sigma}_{>t\leq{t}}(s)\boldsymbol{\Sigma}_{\leq{t}\leq{t}}(s)^{-1}$ and the residual variance as $\boldsymbol{\Sigma}_{>t>t|\leq{t}}(s)=\boldsymbol{\Sigma}_{>t>t}(s)-\boldsymbol{\Sigma}_{>t\leq{t}}(s)\boldsymbol{\Sigma}_{\leq{t}\leq{t}}(s)^{-1}\boldsymbol{\Sigma}_{\leq{t}>t}(s)$. The notation is summarised in Table 1.

[TABLE 1 ABOUT HERE.]

2.2 Estimands

Given this notation, we can write the de jure estimand (estimand E1) at time $t>0$ as $E[Y_{t}(t)-Y_{t}(0)]$. The de facto estimand (estimand E2) is the estimand of interest in this paper and is defined as $E[{Y}_{t}(D_{a})-{Y}_{t}(0)]$. Often, primary interest is in the last observation time, $t=t_{max}$.

2.3 Imputation

Our aim is to impute $\boldsymbol{Y}_{>}t(t)$ in the $D_{a}=t$ subgroup where $t<t_{max}$. Because our main concern is bias, we focus on the mean of the imputation distribution. Variances are considered in Section 2.5.

2.3.1 Initial assumptions and modelling

Assumption A1**.**

$(\boldsymbol{Y}|Z=c)$* is missing at random (MAR).*

Assumption A1 states that the observed outcomes in the control arm are MAR. If there are no missing data before treatment discontinuation then we can also write this

[TABLE]

for all $t$. In other words, treatment discontinuation in the control arm does not depend on future (untreated) outcomes, given the past and present.

Assumption A2**.**

$(\boldsymbol{Y}(t_{max})|Z=a)$* is MAR.*

Assumption A2 states that the counterfactual fully treated outcomes up to time $t$ in the active arm are MAR. If there are no missing data before treatment discontinuation then we can also write this

[TABLE]

for all $t$. In other words, treatment discontinuation in the active arm does not relate to future fully-treated outcomes, given the past and present.

We do not assume that the actual outcomes in the active arm, $(\boldsymbol{Y}|Z=a)$, are MAR. Indeed, this is unlikely to be true, since (if treatment works) treatment discontinuation causally affects actual future outcomes. Thus treatment discontinuation is allowed to relate to future actual outcomes, given the past and present.

Assumption A1 allows us to estimate $\boldsymbol{\mu}(0)$ and $\boldsymbol{\Sigma}(0)$ by applying likelihood-based methods to the observed outcomes in arm $Z=c$. Assumption A2 allows us to estimate $\boldsymbol{\mu}(t_{max})$ and $\boldsymbol{\Sigma}(t_{max})$ by applying likelihood-based methods to the observed outcomes in arm $Z=a$. These two assumptions are therefore enough to identify the de jure estimand. They also allow us to impute intermittent missing data before treatment discontinuation in the active treatment arm, and all missing data in the control arm, so henceforth we assume that the only missing data are actual outcomes after treatment discontinuation in the active arm.

2.3.2 Post-discontinuation imputation using explicit causal assumptions

For the de facto estimand, we need to impute actual outcomes after treatment discontinuation. We consider an individual in the active arm who stops treatment at time $t<t_{max}$ with history $\boldsymbol{Y}_{\leq{t}}$. We use as imputation model a multivariate linear regression model with variance-covariance matrix described in Section 2.5. We write the mean in this model as the sum of three terms:

[TABLE]

where the first term represents the difference between the subgroup who discontinue at time $t$ and the whole group (“selection term”), the second term represents the treatment effect in the whole group (“maintained treatment effect”), and the third term is the mean in the control arm. The third term is identified by assumption A1, and we make the identifying assumptions below for the first two terms on the right hand side of (2).

2.3.3 Selection term

To identify the first term in (2), we assume

Assumption A3**.**

$\boldsymbol{Y}_{>t}(t)|\boldsymbol{Y}_{\leq{t}}(t)$* follows a linear regression for each $t$.*

Assumption A4**.**

$p(D_{a}=t|Z=a,\boldsymbol{Y}(t))=p(D_{a}=t|Z=a,\boldsymbol{Y}_{\leq{t}}(t))$.

A4 appears similar to A2, but whereas A2 states that treatment discontinuation at time $t$ is unaffected by future fully-treated potential outcomes, A4 refers instead to future partly-treated potential outcomes. If there are no missing data before treatment discontinuation then a stronger assumption which implies both A2 and A4 is

[TABLE]

for all $t$ and all $s>t$.

We can now write the selection term as

[TABLE]

where $\boldsymbol{\beta}_{t}(t)$ is the matrix of regression coefficients of $\boldsymbol{Y}_{>t}(t)$ on $\boldsymbol{Y}_{\leq{t}}(t)$ and is not yet identified. However $\boldsymbol{\beta}_{t}(0)$, the regression of $\boldsymbol{Y}_{>t}(0)$ on $\boldsymbol{Y}_{\leq{t}}(0)$, and $\boldsymbol{\beta}_{t}(t_{max})$, the regression of $\boldsymbol{Y}_{>t}(t_{max})$ on $\boldsymbol{Y}_{\leq{t}}(t)$, are identified by assumptions A1 and A2 respectively. If all treatment effects are homogeneous (i.e. if $\boldsymbol{Y}(t)-\boldsymbol{Y}(0)$ does not vary between individuals for any $t$) then $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(t_{max})=\boldsymbol{\beta}_{t}(0)$. Otherwise, they typically differ. Often, we may be willing to assume that $\boldsymbol{\Sigma}(t_{max})=\boldsymbol{\Sigma}(0)$, in which case $\boldsymbol{\beta}_{t}(t_{max})=\boldsymbol{\beta}_{t}(0)$ and this is clearly a sensible choice for $\boldsymbol{\beta}_{t}(t)$. More generally we propose assuming either $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(t_{max})$ or $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(0)$.

2.3.4 Maintained treatment effect

To identify the second term in (2), we make an explicit assumption about the maintained effect of treatment after discontinuation:

Assumption A5**.**

${\mbox{E}}\left[\boldsymbol{Y}_{>t}(t)-\boldsymbol{Y}_{>t}(0)\right]=\boldsymbol{K}_{t}{\mbox{E}}\left[\boldsymbol{Y}_{\leq{t}}(t)-\boldsymbol{Y}_{\leq{t}}(0)\right]$**

or equivalently

[TABLE]

where $\boldsymbol{K}_{t}$ is a $(t_{max}-t)\times(t+1)$ matrix of sensitivity parameters: it is not identified by the data and must be specified by the user. $\boldsymbol{\mu}_{\leq{t}}(0)$ and $\boldsymbol{\mu}_{\leq{t}}(t)$ are identified by assumptions A1 and A2, respectively.

Equation (1) implies $\mu_{0}(t)-\mu_{0}(0)=0$ so the first column of $\boldsymbol{K}_{t}$ does not affect equation (4) and need not be defined. In most cases the post-discontinuation treatment effect depends only on the treatment effect at the time of discontinuation and not on treatment effects at earlier times, and therefore $\boldsymbol{K}_{t}$ needs to have non-zero elements only in the final column: our software implementation below relies on this assumption. For tipping point sensitivity analyses, we consider two causal models with just one parameter in $\boldsymbol{K}_{t}$:

[TABLE]

for $s>t$, where $v_{s},v_{t}$ are the times (on a suitable scale) of observations $s,t$. The maintained treatment effect after treatment discontinuation is constant in model (5) but decays exponentially in model (6), being multiplied by $k_{1}$ in every period, where $0\leq k_{1}\leq 1$. A combined model is

[TABLE]

2.4 Comparison with reference-based imputation in active arm

For imputing $\boldsymbol{Y}_{>t}(t)$, Carpenter et al. (2013, 2014) specified four conditional distributions of post-discontinuation outcomes given pre-discontinuation outcomes under assumptions A1, A2 and A3. The means of these imputation distributions in potential outcomes notation are given in Table 2 for the $D_{a}=t$ subgroup for any $t<t_{max}$. In Table 2, $\boldsymbol{C}_{t}$ is a ‘carry-forward’ $(t_{max}-t)\times(t+1)$ matrix containing $t$ columns of zeroes and a final column of ones, so that $\boldsymbol{C}_{t}\boldsymbol{\mu}_{\leq{t}}(t)$ is a column vector containing $t_{max}-t$ copies of $\mu_{t}(t)$. The expression for CR would usually be written $\boldsymbol{\beta}_{t}(0)\left\{\boldsymbol{Y}_{\leq{t}}-\boldsymbol{\mu}_{\leq{t}}(0)\right\}+\boldsymbol{\mu}_{>t}(0)$ but has been expanded here to demonstrate similarity with the other expressions. The corresponding expression for the causal model is derived by substituting equations (3) and (4) into (2) and is given in the last row of Table 2.

[TABLE 2 ABOUT HERE]

In these expressions, the first term again represents a selection effect describing how an individual with $D_{a}=t$ differs systematically from the whole group. The second term represents the effect of treatment up to time $t$ (compared with not receiving the active arm treatment at all) on the outcome, except for LMCF. From Table 2, it is clear that RBI methods other than LMCF correspond to particular choices of $\boldsymbol{K}_{t}$ with $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(0)$: $\boldsymbol{K}_{t}=\boldsymbol{0}$ for J2R, $\boldsymbol{K}_{t}=\boldsymbol{C}_{t}$ for CIR, and $\boldsymbol{K}_{t}=\boldsymbol{\beta}_{t}(0)$ for CR. This makes precise the statement of Mallinckrodt et al. (2012) that, under CIR, CR and J2R, the $D_{a}=t$ subgroup have the treatment effect at time $t$ maintained, diminished and eliminated respectively at time $t_{max}$. The LMCF method does not correspond to the causal model. Table 2 assumes that the variance-covariance matrix is taken from the control arm; it could also be taken from the active treatment arm, in which case the same equivalences apply, but with $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(t_{max})$.

A useful expression (derived in Appendix A) for the de facto estimand at time $t_{max}$ is

[TABLE]

where $\alpha_{t}=p(D_{a}=t)$ expresses the discontinuation distribution in the active treatment arm; $\delta_{t}={\mbox{E}}\left[Y_{t}(t)-Y_{t}(0)\right]$ is the de jure estimand, with $\delta_{0}=0$; $\boldsymbol{\delta}_{t}=(\delta_{0},\delta_{1},\ldots,\delta_{t})^{\prime}$; and $\boldsymbol{e}_{t}=(0,...,0,1)^{\prime}$ of length $t_{max}-t$. Similar expressions apply at times $t<t_{max}$. From (8), the bias in the de facto estimate associated with wrongly assuming $\boldsymbol{\beta}_{t}(t)$ to equal $\tilde{\boldsymbol{\beta}}_{t}(t)$ is $\sum_{t<t_{max}}\alpha_{t}\boldsymbol{e}_{t}^{\prime}(\tilde{\boldsymbol{\beta}}_{t}(t)-\boldsymbol{\beta}_{t}(t)){\mbox{E}}\left[\boldsymbol{Y}_{\leq{t}}(t|D_{a}=t)-\boldsymbol{Y}_{\leq{t}}(t)\right]$, and the bias associated with wrongly assuming $\boldsymbol{K}_{t}$ to equal $\tilde{\boldsymbol{K}_{t}}$ is $\sum_{t<t_{max}}\alpha_{t}\boldsymbol{e}_{t}^{\prime}(\tilde{\boldsymbol{K}_{t}}-\boldsymbol{K}_{t})\boldsymbol{\delta}_{t}$.

2.5 Implementation of the causal model

2.5.1 Option 1: directly using the SAS 5 macros

The SAS 5 macros were developed by James Roger to perform multiple imputation under the RBI methods and are available at the web page (on www.missingdata.org.uk) of the DIA working group for missing data. We modified the Part2A macro to impute under the causal model with

[TABLE]

where $k$ is a scalar that may vary between participants. This enables causal model (5) to be implemented by setting $k=k_{0}$, the same for all participants. When interest is in the outcome at time $t_{max}$, causal model (6) can be implemented by setting $k=k_{1}^{t_{max}-D_{a}}$, which varies across participants with different values of $D_{a}$. The modified macro is available on the DIA working group web page and sample code is provided in Appendix B. By default the variance-covariance matrices in the two arms, $\boldsymbol{\Sigma}(t_{max})$ and $\boldsymbol{\Sigma}(0)$, are assumed equal, but the user can specify them to be unequal, which is the case we consider.

2.5.2 Option 2: indirectly using the SAS 5 macros (based on J2R)

Equation (8) for the de facto estimand at time $t_{max}$ can be written as

[TABLE]

where $\delta_{J2R}$ is the J2R estimand. Thus we can run the J2R method to estimate the first term, fit a mixed model for repeated measures to estimate the $\delta_{t}$, directly estimate the $\alpha_{t}$, and assume values for $\boldsymbol{K}_{t}$.

To approximate the standard error, we propose using the Rubin’s rule (RR) estimate of standard error from J2R, as we will see in Section 3 that the standard error does not vary much with $\boldsymbol{K}_{t}$.

2.5.3 Option 3: using Mixed Model Repeated Measures (MMRM) estimates or standard Multiple Imputation (MI)-based estimates

If we are willing to assume either that $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(t_{max})$ or that there is no selection effect, then the third term in expression (8) vanishes and the de facto estimand at time $t_{max}$ can be written as

[TABLE]

That is, the de facto estimate at time $t_{max}$ can be estimated as a linear combination of de jure MAR-based estimates at times $t\leq t_{max}$. To find the standard error, we could ignore uncertainty in the $\hat{\alpha}_{t}$ and use the standard error of the linear combination, or we could do better by applying the delta approximation method (Oehlert, 1992) as done by Liu and Pang (2016).

3 Simulation

We performed a simulation study to demonstrate equivalence of the RBI methods and the proposed causal model for estimating the treatment effect at the final time point, and to assess the impact of mis-specification of $\boldsymbol{K}_{t}$ and $\boldsymbol{\beta}_{t}(t)$.

3.1 Design

We consider a RCT with one baseline observation and two post-baseline follow-ups during the treatment period (that is, $t_{max}=2$). Potential untreated outcomes $\boldsymbol{Y}(0)$ in both arms were generated from a multivariate normal (MVN) distribution with mean $\boldsymbol{\mu}(0)=(10,12,14)^{^{\prime}}$ and variance-covariance matrix $AR(1)$ with standard deviation 3 and lag-one correlation $0.5$. Potential fully-treated outcomes were derived using the causal models $Y_{1}(1)-Y_{1}(0)=\delta_{1}+u_{1}$ and $Y_{2}(2)-Y_{2}(0)=\delta_{2}+u_{1}$, where $\delta_{1}=1$ and $\delta_{2}=2$ are average treatment effects and $u_{1}$ may vary between participants to give treatment effect heterogeneity. In the control arm, we assumed that all participants had complete data. Observed outcomes in the control arm were set equal to untreated potential outcomes.

In the active arm, we assumed that participants either completed the entire period of treatment and provided complete measurements ($D_{a}=2$), or completed only the first period of treatment and failed to provide measurements at visit 2 ($D_{a}=1$). The treatment discontinuation mechanism was ${\rm logit}\ p(D_{a}=1|Y_{0},Y_{1}(1))=\tau_{0}+\tau_{1}Y_{1}(1)$. Observed outcomes were set equal to the potential fully-treated outcomes. We used two values of $\tau_{1}$: [math] for missing completely at random (MCAR) and $1$ for MAR. $\tau_{0}$ was chosen to give $p(D_{a}=1)=0.5$.

We explored the implications of differences between $\boldsymbol{\beta}_{t}(t_{max})$ and $\boldsymbol{\beta}_{t}(0)$. Since missing data occur only in the $D_{a}=1$ group, we needed only consider differences for $t=1$, that is, between $\boldsymbol{\beta}_{1}(2)$ and $\boldsymbol{\beta}_{1}(0)$. The distribution above for $\boldsymbol{Y}(0)$ implies $\boldsymbol{\beta}_{1}(0)=(0,0.50)^{\prime}$. We used two treatment effect models: homogeneity, $u_{1}=0$ implying $\boldsymbol{\beta}_{1}(2)=(0,0.50)^{\prime}=\boldsymbol{\beta}_{1}(0)$; and heterogeneity, $u_{1}\sim N(0,2.5^{2})$ implying $\boldsymbol{\beta}_{1}(2)=(-0.12,0.74)^{\prime}\neq\boldsymbol{\beta}_{1}(0)$. The latter model resulted in slightly more heterogeneity in the active arm than in the control arm. We considered 250 participants per arm and 1000 repetitions for the simulation study.

This completes the mechanism for simulating the observed data. However, the de facto estimand depends on the unobserved data including the potential partly-treated outcomes $Y_{2}(1)$ whose distribution has not yet been specified. We generated $Y_{2}(1)$ from the causal model $Y_{2}(1)-Y_{2}(0)=k\delta_{1}+u_{2}$ where $k$ took values 0, 0.5, 0.74 or 1, and where the marginal distribution of $u_{2}$ was the same as that of $u_{1}$ with either (a) $corr(u_{1},u_{2})=0.5$, so that $\boldsymbol{\beta}_{1}(1)=\boldsymbol{\beta}_{1}(0)$, or (b) $corr(u_{1},u_{2})=1$, so that $\boldsymbol{\beta}_{1}(1)=\boldsymbol{\beta}_{1}(2)$. For the treatment effect homogeneity model, options (a) and (b) are equivalent.

For each mechanism for simulating the observed data, we analysed the data in three ways. Firstly (reported in panel A of Tables 3 and 4), we generated the complete data by setting $Y_{2}=Y_{2}(D_{a})$ under the various values of $k$ and the two choices of $\boldsymbol{\beta}_{1}(1)$ described above, and analysed the complete data. Secondly (reported in panel B of Tables 3 and 4) we imputed the missing data using the causal model. Here $\boldsymbol{K}_{1}$ was a $1\times 2$ vector whose first element was irrelevant and whose second element was the scalar $k$. Using $\tilde{k}$ and $\tilde{\boldsymbol{\beta}}_{1}(1)$ to denote assumed values of $k$ and $\boldsymbol{\beta}_{1}(1)$, we took the same range of values for $\tilde{k}$ as proposed for $k$ above, and the same choices of $\tilde{\boldsymbol{\beta}}_{1}(1)$. Thirdly (reported in panel C of Tables 3 and 4) we imputed the missing data using the reference-based imputation methods CR, CIR and J2R with the two choices of variance-covariance matrix. In all cases, we estimated the treatment effect from a linear regression of $Y_{2}$ on randomised arm and baseline $Y_{0}$. With imputed data, standard errors were computed using Rubin’s rules.

3.2 Results

Table 3 displays the average estimated treatment difference at the final visit for each data generating mechanism.

Comparing panels A and B shows that the causal model imputation methods result in unbiased estimates when the assumed values of $\tilde{\boldsymbol{\beta}}_{1}(1)$ and $\tilde{k}$ agree with the true values of ${\boldsymbol{\beta}}_{1}(1)$ and ${k}$.

Comparing panels B and C shows that the RBI estimates with variance-covariance matrix drawn from the control arm are special cases of the causal model estimates with $\boldsymbol{\beta}_{1}(1)=\boldsymbol{\beta}_{1}(0)$, and the RBI estimates with variance-covariance matrix drawn from the active arm are special cases of the causal model estimates with $\boldsymbol{\beta}_{1}(1)=\boldsymbol{\beta}_{1}(2)$. Specifically, J2R corresponds to $\tilde{k}=0$, CIR corresponds to $\tilde{k}=1$, and CR corresponds to $\tilde{k}=$ the second element of $\boldsymbol{\beta}_{1}(1)$ which is 0.50 or 0.74 depending on the data generating mechanism.

Comparing different choices of $\tilde{\boldsymbol{\beta}}_{1}(1)$ when the observed data had either no selection effect (i.e. under MCAR) or $\boldsymbol{\beta}_{1}(0)=\boldsymbol{\beta}_{1}(2)$ (that is, in the first three data generating mechanisms), we see that choice of $\tilde{\boldsymbol{\beta}}_{1}(1)$ does not affect estimates, as expected from Section 2.4. Sensitivity to choice of $\tilde{\boldsymbol{\beta}}_{1}(1)$ was observed in the fourth data generating mechanism (MAR with $\boldsymbol{\beta}_{1}(0)\neq\boldsymbol{\beta}_{1}(2)$): mean causal model estimates were reduced by 0.29 by assuming $\boldsymbol{\beta}_{1}=\boldsymbol{\beta}_{1}(0)$ instead of $\boldsymbol{\beta}_{1}=\boldsymbol{\beta}_{1}(2)$. Sensitivity to choice of $\tilde{k}$ was the same for all values of $\boldsymbol{\beta}_{1}$ (see Section 2.4): for example, assuming $\tilde{k}=0$ instead of $\tilde{k}=1$ reduced mean causal model estimates by 0.50 irrespective of the value of $\boldsymbol{\beta}_{1}$.

[TABLE 3 ABOUT HERE]

Table 4 displays the average standard error (SE) (the average of the 1000 SEs) and the empirical SE (the sample standard deviation of the 1000 point estimates) for the treatment difference at the final visit. Empirical and average SEs for J2R and CIR are similar to those for corresponding causal model estimates. The SEs for CR are slightly larger than that for causal model with $\tilde{k}=$ the second element of $\boldsymbol{\beta}_{1}$ because $\tilde{k}$ is an assumed value while $\boldsymbol{\beta}_{1}$ is estimated. With MAR data, the larger average and empirical SEs due to using the variance-covariance matrix from the active arm rather than from the control arm is mainly because of no missing data in the control arm and larger heterogeneity in the active arm than in the control arm. More importantly, as shown in Seaman et al. (2014), the results confirm that both RBI and causal model methods give (1) smaller empirical SEs than the estimator based on the complete data, and (2) larger average SEs (estimated using RR) than the empirical SEs of the methods and the empirical SEs based on the complete data. We comment on these observations in the discussion.

[TABLE 4 ABOUT HERE]

4 Examples

We use two example data sets from randomised, double blind, parallel-group studies comparing active drug with placebo. The first is from a trial of 172 participants with major depressive disorders, taken from the DIA page of www.missingdata.org.uk, and used in the DIA working group to demonstrate various missing data related analytical methods. The outcome variable is the 17-item Hamilton Depression Rating Scale, HAMD17. The second, kindly supplied by Devan Mehrotra, is from a pain trial with a pain score as outcome.

In the HAMD17 trial, 76% (64/84) and 74% (65/88) of randomised participants completed the final (fourth) visit in the drug and the placebo arms respectively. In the pain score trial, the completion rate at the final (sixth) visit was 70% (47/67) and 67% (36/54) in the drug and the placebo arms respectively. In both trials, participants were not followed up after treatment discontinuation. The observed trajectory means and the frequency of dropout patterns in each trial are shown in Figure 1.

[FIGURE 1 ABOUT HERE]

We used the SAS 5 macros for implementing the RBI and the causal models (Section 2.5.1). For the RBI methods, we assumed participants in the drug arm were treated similarly to the placebo arm after discontinuing the drug. To construct the joint distribution of pre- and post-discontinuation drug-arm data under the RBI methods, we first used the variance-covariance matrix from the placebo arm (RBI1 analyses) and then repeated the methods with the variance-covariance matrix from the drug arm (RBI2 analyses).

Table 5 shows the estimated treatment effect on HAMD17 and pain score at the final visit from standard MI, MMRM and RBI methods. The standard MI and MMRM methods estimate the de jure estimand. These differ slightly for HAMD17 because of a small incompatibility between the imputation and analysis models: the imputation model used all visits to estimate a common the effect of the baseline covariate PoolInv, but the analysis model used only the final visit. The RBI methods estimate the de facto estimand and showed, as expected, treatment estimates of smaller magnitude than the de jure estimand, with J2R giving the smallest magnitude of treatment effect followed by CR. Using the variance-covariance matrix from the drug arm rather than from the placebo arm gave slightly more conservative estimates.

[TABLE 5 ABOUT HERE]

We demonstrate tipping point sensitivity analyses using causal models (5) and (6). In model (5), a fraction $k_{0}$ of the treatment effect is maintained at all times after discontinuation. Figure 2 shows the de facto estimates and 95% CI over a range of $k_{0}$ from -0.5 to 2.5. As shown in the theory and the simulation results, the J2R and CIR estimates correspond to using the causal model with $k_{0}=0$ (i.e., when participants in the active arm lost their treatment effect after the discontinuation) and $k_{0}=1$ (i.e., when participants in the active arm maintained their treatment effect after the discontinuation), respectively. Values $k_{0}<0$ mean that the effect of treatment after discontinuation is harmful, while values $k_{0}>1$ mean that the effect of treatment after discontinuation is greater than before discontinuation. The tipping point analysis on HAMD17 shows that statistical significance is lost when $k_{0}<0$ (with variance-covariance from the placebo arm) or $k_{0}<0.05$ (with variance-covariance from the active arm). In both cases, this suggests that the de facto estimate of treatment effect on HAMD17 is non-significant only if any benefit of the active treatment is lost immediately following discontinuation. For the pain score, statistical significance is lost when $k_{0}<1.1$ or $k_{0}<1.3$ (depending on whether the variance-covariance matrix is taken from the placebo or the drug arm, respectively). This suggests that, in order for the de facto estimate of treatment effect to be statistically significant, there would need to be a delayed benefit such that the treatment effect was greater after discontinuation than before discontinuation. In both trials, comparing Table 5 with Figure 2 shows that the MAR analyses give estimates of the de jure estimand that are numerically similar to the causal model estimates of the de facto estimand when values around $k_{0}=2$ are assumed.

[FIGURE 2 ABOUT HERE]

In model (6), the treatment effect decays exponentially after discontinuation. Here, $k_{1}=0$ for J2R and $1$ for CIR. Figure 3 shows the de facto estimates of treatment effect at the final visit and its 95% CI from the causal model over a range of $k_{1}$. This model does not accommodate the effect of treatment after discontinuation being either harmful or greater than before discontinuation, and because of the more limited range of $k_{1}$, the tipping point is not reached: all results are statistically significant for HAMD17 and not significant for the pain score.

[FIGURE 3 ABOUT HERE]

5 Discussion

We have considered longitudinal RCTs with quantitative outcomes in which participants who discontinue an active treatment are not followed up thereafter, but are assumed to receive a standard treatment which is often the control treatment. We have focused on estimating the effect of assignment to treatment in the actual treatment circumstances of the trial (de facto estimand) rather than the treatment effect if all participants had tolerated or adhered to trial protocol (de jure estimand). We have investigated the underlying assumptions of the RBI methods and proposed a generalised causal modelling approach to account for treatment discontinuation in the estimation of the de facto estimand. The proposed causal model makes an explicit assumption about the causal effect of a given treatment history, and provides flexibility to perform sensitivity analyses to the causal assumption.

The proposed causal model specifies how much of the treatment effect is maintained after treatment discontinuation, which we represent by the parameter matrix $\boldsymbol{K}_{t}$. We illustrated this with two examples of $\boldsymbol{K}_{t}$: (5) with the maintained treatment effect independent of time since discontinuation, and (6) with the maintained treatment effect decaying exponentially with visits since discontinuation. A simple extension would allow $\boldsymbol{K}_{t}$ to depend on the reason for treatment discontinuation. Ideally sponsors should justify the choice of $\boldsymbol{K}_{t}$ in the trial protocol based on the nature of the trial and the treatments.

The choice of $\boldsymbol{\beta}_{t}(t)$ (reflecting within-subject dependence of post-discontinuation outcomes on pre-discontinuation outcomes) should similarly be pre-specified; it is hard to recommend a single choice and perhaps both $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(0)$ and $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(t_{max})$ should be implemented. If the analyst is willing to assume equal variance-covariance matrices across trial arms then the situation is simpler and $\boldsymbol{\beta}_{t}(t)=\boldsymbol{\beta}_{t}(0)=\boldsymbol{\beta}_{t}(t_{max})$ is the obvious choice.

Our model has been presented for the case of active arm treatment discontinuation, where subjects who discontinue do not then receive rescue medication (in the case of placebo controlled trials), or only receive standard of care (in the case where the control arm is standard of care). It can be extended to handle the control arm starting active treatment: an assumption like A5 still holds, but the $\boldsymbol{K}_{t}$ matrix must be replaced by assumptions about how the treatment effect builds up over time. We are currently working on extending the model to trials comparing two active treatments. An unresolved problem is how to handle initiation of rescue medications.

We have focussed on varying assumption A5, but we should also assess a number of other assumptions. The MAR assumptions A1 and A2, and the related assumption A4, could be made more plausible if the model could be extended to include further time-dependent covariates. Alternatively one could explore sensitivity to these assumptions by methods like those of Ratitch et al. (2013). It is less clear how to assess departures from the linearity assumption A3.

We have focussed on quantitative outcomes. Extension to other outcomes would be useful, and we are currently working on a causal version of the method of Keene et al. (2014) for recurrent events.

We proposed approximations to the RR variance in options 2 and 3 of Section 2.5. From the simulation and empirical results, we see that the SEs are not much affected by the value of $\boldsymbol{K}_{t}$. This means that our approximations are adequate if we want the RR variance but not if we want the repeated sampling variance. The repeated sampling variance of the estimated treatment effect tends to be smaller than the RR estimate of variance for a given $\boldsymbol{K}_{t}$. The repeated sampling variance can be approximated in practice using the delta method (Oehlert, 1992; Liu and Pang, 2016). Carpenter et al. (2014) argue that the repeated sampling variance is not appropriate, since it is typically smaller than the complete-data variance (to an extent which strongly depends on the value of $\boldsymbol{K}_{t}$). They also argue that the RR estimate of variance of the treatment effect is larger than the the complete-data variance, because of the information lost due to the missing data, and this makes it an appropriate variance (Carpenter et al., 2014; Cro et al., 2015). We point out that the type I error rate is correct for the repeated sampling variance and too small for the RR variance, meaning that the RR variance carries a loss of power: therefore the repeated-sampling variance may be appropriate for a primary analysis.

In summary, whilst multiple imputation is an attractive and powerful method for handling missing data in both experimental and observational studies, it is not always clear what estimand is being targeted or what assumptions are being made about how outcomes for subjects who discontinue randomised treatment relate to those who remain on study. The recent estimands debate (ICH steering committee, 2014) has led to a growing recognition that more complex estimation approaches that do not rely on randomisation are needed to handle post-randomisation events that lead to missing data, and there are calls for causal inference methods to become more widely adopted (e.g. Akacha et al. (2016); Little and Kang (2015)). We join this call to encourage greater understanding and application of ideas from the causal inference literature to help support the definition and estimation of estimands of interest in a randomised clinical trial. We hope that this paper illustrates how a causal inference framework can provide clarity and rigour in stating estimands, stating assumptions, and performing estimation.

Acknowledgements

IRW and RJ were supported by the Medical Research Council (Unit Programme number U105260558). RJ was also supported by GSK. Nicky Best is an employee and shareholder of GSK. We thank James Roger for very helpful comments on an earlier draft and for allowing us to modify his software, and Devan Mehrotra for providing the pain score dataset.

Appendix Appendix A Proof of equation (8)

We first write the de facto estimand at time $t_{max}$ as a weighted average over treatment discontinuation times:

[TABLE]

where we define $\boldsymbol{\mu}(s|D_{a}=s)=E[\boldsymbol{Y}(s)|D_{a}=s]$ etc.

We next write the terms on the right hand side of (10) for $t<t_{max}$ using the result in Table 2:

[TABLE]

where $\boldsymbol{e}_{t}$ is used to extract the $t_{max}$th element.

Selection effects cancel out over treatment discontinuation times $t$. To use this, note that A3 implies

[TABLE]

for each $t<t_{max}$, and the sum of the right hand side over $t$ (weighted by $\alpha_{t}$) is zero. We sum the left hand side (weighted by $\alpha_{t}$) and take the last element, giving

[TABLE]

and hence

[TABLE]

Finally, substituting equations (Appendix A) and (11) into (10) gives

[TABLE]

which is equivalent to equation (8) in the main text.

Appendix Appendix B Causal model implementation using the SAS 5 macros

The modified macro file ‘Part2A_38_causal.sas’ and the HAMD17 data are available on the DIA section of the website www.missingdata.org.uk. The HAMD17 data are stored as a SAS data set named Chapter15_example which we assume has been renamed HAMD.

The modified Part2A macro allows the user to impute under the causal model by specifying option Method=Causal. The value of $k$ in equation (9) is specified by Causalk= or defaults to 1 giving the CIR estimate. $k$ can be specified as a scalar (a constant $k$ for all participants; e.g. Causalk=0.5) or as a variable that specifies $k$ for each participant (e.g. Causalk=varname). If $k$ is specified using a variable, that variable must be identified in macro Part1A using the option id=varname. If Method=Causal is not specified then RBI methods are implemented and Causalk= is ignored.

Causal model (5) with $k_{0}=0.5$ can be fitted to the HAMD17 data by:

%Part1A(Jobname=Example, Data=HAMD, Subject=Patient, Response=Change, Time=Visit, Treat=Therapy, Catcov=PoolInv, Covbytime=Basval, Covgroup=Therapy); %Part1B(Jobname=Example, Ndraws=100, Thin=100, Seed=12345); %Part2A(Jobname=Example_causal, Inname=Example, Method=Causal, Causalk=.5, Ref=Placebo, VCMethod=Ref);

The option Covgroup=Therapy in the %Part1A call specifies that the variance-covariance matrix is estimated separately by arm, but Catcov=PoolInv and Covbytime=Basval specify that the coefficients of the baseline covariates are assumed equal across arms.

Causal model (6) with $k_{1}=0.5$ can be fitted to the HAMD17 data after creating a new variable k1power to hold the individual’s value of $k_{1}^{v_{s}-v_{t}}$ where $s=t_{max}$ and $t=D_{a}$. In these data Visit is coded 0, 4, 5, 6, 7, so we set $v_{s}=7$ and compute $v_{t}$ as the last value of Visit:

proc sql; create table HAMD2 as select , 0.5*(7-max(Visit)) as k1power from HAMD group by Patient; quit;

Although k1power is computed for all individuals, the next steps ignore values computed in the $Z=c$ arm:

%Part1A(Jobname=Example, Data=HAMD2, Subject=Patient, Response=Change, Time=Visit, Treat=Therapy, Catcov=PoolInv, Covbytime=Basval, Covgroup=Therapy, Id=k1power); %Part1B(Jobname=Example, Ndraws=100, Thin=100, Seed=12345); %Part2A(Jobname=Example_causal, Inname=Example, Method=Causal, Causalk=k1power, Ref=Placebo, VCMethod=Ref);

%Part1A and %Part1B need not be repeated for different values of $k_{0}$ or $k_{1}$. The imputed data sets in either case are analysed using

%Part2B(Jobname=Example_causal, Seed=54321); %Part3(Jobname=Example_causal, Anref=Placebo, Label=Causal with sigma from reference);