Instrumental variable estimation of early treatment effect in randomized screening trials

We introduce the necessary notation in the context of the US National Lung Screening Trial, where $53,452$ heavy smokers, currently smoking or having quit within the last $15$ years, with a smoking history of $30+$ pack years, aged between $55-74$, were assigned at random to undergo three annual screening examinations or control. The participants were asymptomatic of lung cancer at the time of randomization. Participants in the screening arm received three annual low-dose helical CT scans, whereas participants in the control arm received a standard chest X-ray, and both arms were followed up for $7$ years for the main analysis (NLST Research Team, 2011). There were $649$ early diagnoses in the screening arm and $469$ and $552$ cancer deaths in the screening arm and control arm, respectively. In the following, we assume the X-ray screening to have negligible mortality benefits, as has been demonstrated by Oken et al. (2011), and equate it to a placebo. We characterize the causal effect of screening-induced early treatments through a simplified structural multi-state model where the mortality benefits manifest after an early detection through a CT scan. From the healthy state (numbered as state 1), participants in the screening arm have three possible transitions: early diagnosed through CT screening (state 2), cancer-specific death (state 3) and other-cause death (state 4). From the early diagnosed state, the participants can move to cancer-specific death or other-cause death states, with potentially different transition rates. We assume that the control arm does not have early diagnosis through CT screening, and thus only direct transitions from state 1 to states 3 and 4 are possible. For concise notation, we use counting processes. We note that since each particular type of transition in a multi-state model can be characterized through a competing risks model, the transition intensities are analogous to the counterfactual cause-specific hazards, as outlined by Young et al. (2020), but conditional on the event histories taken place so far.

Let $N^{z}_{kl}(t)\in\{0,1\}$ be the potential (underlying, in the absence of censoring) counting process that counts the transitions of a given individual from state $k$ to $l$ by time $t$ in arm $z$, where $k\neq l;\;k,l\in\{1,2,3,4\}$ and $z=1$, or $0$ indicate the screening or control arms, respectively. In the current context each transition can take place at most once, and thus each counting process only counts to one, with states 3 and 4 being absorbing states. The corresponding observed process is formulated as $N_{kl}(t)=ZN_{kl}^{1}(t)+(1-Z)N_{kl}^{0}(t)$, where $Z$ is the arm that the patient actually was assigned to in the randomized screening trial (counterfactual consistency). The counting processes are characterized by the transition intensities defined as

$\displaystyle\lambda_{kl}^{z}(t)=\lim_{\Delta t\to 0}\frac{P(N_{kl}^{z}(t+\Delta t)-N_{kl}^{z}(t^{-})=1\mid H^{z}(t^{-}))}{\Delta t}$

where $H^{z}(t^{-})$ indicates the history of states until time $t^{-}$. We note that the model could be alternatively specified in discrete time through the conditional probabilities $P(N_{kl}^{z}(t+\Delta t)-N_{kl}^{z}(t^{-})=1\mid H^{z}(t^{-}))\approx\lambda_{kl}^{z}(t)\Delta t$, but we will use continuous time notation for simplicity. The corresponding cumulative transition intensity function is defined as $\Lambda_{kl}^{z}(t)=\int_{s}^{t}\lambda_{kl}^{z}(t)\,\textrm{d}t$. The transition intensities can be collected into a matrix

$$\lambda^{z}(t)=\blockarray{c@{\hskip 1.0pt}rrrr@{\hskip 4.0pt}}&1234\\ \block{r@{\hskip 1.0pt}|@{\hskip 1.0pt}|@{\hskip 1.0pt}rrrr@{\hskip 1.0pt}|@{\hskip 1.0pt}|}1.\lambda_{12}^{z}(t)\lambda_{13}^{z}(t)\lambda_{14}^{z}(t)\\ 20.\lambda_{23}^{z}(t)\lambda_{24}^{z}(t)\\ 300.0\\ 4000.\\$$

In the conventional screening trials, the ITS estimand would be defined by comparing the outcomes between the two arms at time $t$ after sufficiently long follow-up, in terms of the cumulative incidences for cancer-specific mortality. For example, the absolute reduction in cancer-specific mortality risk would be measured by

$$E[N_{13}^{0}(t)]-E[N_{13}^{1}(t)+N_{23}^{1}(t)]$$

(1)

and the proportional reduction by

$$1-\frac{E[N_{13}^{1}(t)+N_{23}^{1}(t)]}{E[N_{13}^{0}(t)]}$$

(2)

The estimation of these quantities does not require a multi-state model, as the components could be estimated through empirical cumulative incidences in the two arms. However, the multi-state model is needed to specify and estimate the causal effect of the early treatment among the early detectable subgroup. This subgroup is latent at baseline and accumulates in the screening arm during the follow-up of the trial through the repeated screening examinations. Since screening itself is not an intervention, the mortality benefits can manifest only after an early detection. To specify the corresponding causal effect, we need to introduce a well-defined intervention.

Following and extending Miettinen (2015), we can envision a trial that has similar screening regimen as the screening arm of the conventional trial, but instead of randomizing into screening and non-screening, screens everyone, and at the time of an early detection, randomizes into referral to early treatment vs delayed treatment through withholding the screening result. Delayed treatment refers treatments following a later symptomatic diagnosis. Since the early detection can take place only once for a given individual, the randomization for the assignment in the event of early detection can already take place at study baseline, which is relevant for the indexing of the corresponding potential outcomes. A schematic illustration of the hypothetical trial is presented in Figure 1, and the corresponding multi-state model in Figure 2. While this kind of trial could not be implemented in practice, it is important as a thought experiment to define the causal quantities that we are interested in.

Among $N_{12}^{1}(t)=1$ we define $N^{1r}_{2l}(t)\in\{0,1\}$, $l\in\{3,4\}$ to be potential (underlying) counting processes for subsequent death corresponding to early ($r=1$) versus delayed ($r=0$) treatment, characterized by transition intensities

$\displaystyle\lambda_{2l}^{1r}(t)=\lim_{\Delta t\to 0}\frac{P(N_{2l}^{1r}(t+\Delta t)-N_{2l}^{1r}(t^{-})=1\mid H^{1}(t^{-}))}{\Delta t}.$

The modified transition matrix for screening detectable subgroup is now

$$\lambda^{1r}(t)=\blockarray{c@{\hskip 1.0pt}rrrr@{\hskip 4.0pt}}&1234\\ \block{r@{\hskip 1.0pt}|@{\hskip 1.0pt}|@{\hskip 1.0pt}rrrr@{\hskip 1.0pt}|@{\hskip 1.0pt}|}1.\lambda_{12}^{1}(t)\lambda_{13}^{1}(t)\lambda_{14}^{1}(t)\\ 20.\lambda_{23}^{1r}(t)\lambda_{24}^{1r}(t)\\ 300.0\\ 4000.\\$$

Under this model, the absolute cancer mortality risk reduction in the subgroup could be measured by

$$E[N_{23}^{10}(t)-N_{23}^{11}(t)\mid N_{12}^{1}(t)=1]$$

(3)

and the proportional reduction by

$$1-\frac{E[N_{23}^{11}(t)=1\mid N_{12}^{1}(t)=1]}{E[N_{23}^{10}(t)=1\mid N_{12}^{1}(t)=1]}.$$

(4)

We note that these quantities are distinct from (1) and (2) as they are defined in the early diagnosed subpopulation. Conditioning on this subgroup membership also conditions on survival until the early diagnosis, so (3) and (4) measure mortality reductions after the early diagnosis. In a trial where the screening is discontinued at some point (as it was in NLST after three annual rounds), and the follow-up is long enough so that all the mortality benefits have been realized, these quantities approximate the absolute and proportional case-fatality reduction, contrasting the screening-induced early versus delayed treatment among the screening detectable cases. The estimation of these quantities was discussed by Saha et al. (2018).

Herein we aim to measure the effect by contrasting the transition intensities $\lambda_{23}^{11}(t)$ and $\lambda_{23}^{10}(t)$. We note that the formulation in terms of transition intensities does not yet assume anything about the early treatment effect, for example generally this effect could be a function of time and/or function of the time of the early diagnosis. However, we can simplify the model by characterizing the effect as a constant ratio $\lambda_{23}^{10}(t)=\theta\lambda_{23}^{11}(t)$, characterizing the early treatment effect on cancer-specific mortality through a single parameter, corresponding to the proportionality of the two cause-specific mortality hazards. In addition, although this is not generally required, we take $\lambda_{24}^{10}(t)=\lambda_{24}^{11}(t)$, meaning that the early treatments do not have effect on other-cause mortality. In what follows, we focus on the estimation of the quantity $\theta$ as our effect measure of interest. We note that we parametrized this as the ratio of the mortality under delayed vs. early treatment to be consistent with the earlier defined risk reduction measures. Alternatively, the inverse of this can be reported. Herein we take the time scale for all of the transition intensities to be the time since the baseline/study entry. However, it would also be equally possible to use the time since early detection as the time scale for the transitions after the early detection, and to define proportionality of the effect on this time scale.

Since the transition intensities $\lambda_{2l}^{10}(t)$ are not directly estimable from data observed under a conventional randomized screening trial, the estimation requires instrumental variable type methods, using the randomized assignment in the conventional trial as an instrument. We will propose estimation methods in Section 3, but will first discuss the use of the model for simulation, and outline the assumptions required for the instrumental variable estimation.

In the methods proposed in Section 3, the random assignment in the conventional screening trial is used as an instrumental variable. We outline the identifying assumptions here as they are also used to generate simulated data from the model. In particular, we assume that (i) $N_{12}^{1}(t)=1\Rightarrow N_{1l}^{0}(t)=N_{2l}^{10}(t)$, $l\in\{3,4\}$, meaning that for someone who would get early detected by time $t$, but not receive early treatment, the potential mortality outcome would be the same as in the control arm of the conventional trial. Further, we assume that (ii) $N_{12}^{1}(t)=0\Rightarrow N_{1l}^{0}(t)=N_{1l}^{1}(t)$, $l\in\{3,4\}$, meaning that for someone who would not get early detected, the potential mortality outcome would be the same between the two arms of the conventional trial. Essentially, assumptions (i) and (ii) together state that screening examinations themselves do not have impact on mortality in the absence of the intervention. They broadly correspond to the exclusion restriction assumption of conventional instrumental variable analysis (Angrist et al., 1996), stating that the instrumental variable does not have direct causal effect on the outcome, outside the effect mediated by the intervention. In the present context this assumption could be violated for example if the screening examinations cause the participants to modify their health behaviour in a way that affects their cancer risk. In the NLST also the control arm received screening, so any such effects should be more similar between the two arms, making (i) and (ii) more plausible.

In addition, we assume that (iii) $N_{12}^{1}(t)=1\Rightarrow N_{2l}^{1}(t)=N_{2l}^{11}(t)$, $l\in\{3,4\}$, meaning that in practice in the screening arm early treatment follows early diagnosis, which is reasonable at least in the lung cancer context. From this it follows that $\lambda_{2l}^{1}(t)=\lambda_{2l}^{11}(t)$, meaning that the mortality rates under early treatment following early diagnosis can be estimated from the screening arm of the conventional trial. For the purpose of identification based on data from the two arms of the conventional trial, we also need counterfactual consistency (iv), formulated in Section 2.1. Finally, we take that $\lambda_{12}^{0}(t)=\lambda_{23}^{0}(t)=\lambda_{24}^{0}(t)=0$, meaning that early detection through the screening technology offered in the screening arm is not available in the control arm.

To simulate data from the proposed multi-state model, one needs to fix the transition intensities in the matrix $\lambda^{1r}(t)$, to simulate event histories in the two arms in the hypothetical intervention trial, corresponding to $r=1$ and $r=0$. These event histories are also used to represent event histories in the conventional screening trial. In particular, under condition (iii) the intervention arm ($r=1$) event histories can be directly taken to be event histories in the screening arm of the conventional trial. Also, the control arm $r=0$ event histories can be directly taken to be event histories in the control arm by taking $N_{12}^{0}(t)=0$ and $N_{1l}^{0}(t)=N_{1l}^{1}(t)+N_{2l}^{10}(t)$, $l\in\{3,4\}$, which follows from (i) and (ii). If the outcomes are simulated without covariates, the data generating mechanism and the structural model are the same. However, for the data generation, all the intensities can be made conditional on observed and/or unobserved simulated covariates, in which case the structural model will parametrize the marginal effect of interest. In this case, the marginal effect is not directly specified by the data generating mechanism, but can be approximated by fitting a marginal proportional hazards model to the two arms of the hypothetical trial.

For a simulation study, for example exponential or Weibull functional forms can be assumed for the intensities, or the observable ones ($\lambda_{1l}^{1}(t)$, $l\in\{2,3,4\}$ and $\lambda_{2l}^{11}(t)$, $l\in\{3,4\}$) can be estimated from the screening arm of an existing trial such as the NLST, with the control arm ($r=0$) counterparts obtained by fixing the effect parameter $\theta$ to a value. With the transition intensity matrix $\lambda^{1r}$, simulation can then proceed through usual algorithms for simulating event histories from a multi-state model as described for example in Beyersmann et al. (2011, p.  45). Briefly, this proceeds by simulating the time and type of the next event from a competing risks model, moving to a new state, and continuing the same. Starting from the healthy state, this would mean simulating the time of the first event, $T$, from the total hazard

$\displaystyle\Lambda^{1}(t)=\Lambda_{12}^{1}(t)+\Lambda_{13}^{1}(t)+\Lambda_{14}^{1}(t),$

which specifies the event time distribution $F(t)=P(T\leq t)=1-\exp(-\Lambda^{1}(t))$. Using the inverse cumulative distribution function method, we can take $u\sim U[0,1]$ and solve $\Lambda^{1}(t)+\log(1-u)=0$ with respect to $t$ to get the time. At time $T=t$ the type of the event is then drawn with multinomial probabilities

$\displaystyle\frac{\lambda_{1l}^{1}(t)}{\lambda_{12}^{1}(t)+\lambda_{13}^{1}(t)+\lambda_{14}^{1}(t)}$

for $l\in\{2,3,4\}$. If the new state is terminal ($l\in\{3,4\}$), the algorithm stops, otherwise it continues similarly from state 2 after possible modification of the subsequent transition intensities based on the timing of the first event, where the new total hazard is given by $\Lambda^{1r}(t)=\Lambda_{23}^{1r}(t)+\Lambda_{24}^{1r}(t)$, and the event type probabilities by

$\displaystyle\frac{\lambda_{2l}^{1r}(t)}{\lambda_{23}^{1r}(t)+\lambda_{24}^{1r}(t)}$

for $l\in\{3,4\}$. We used the NLST as the basis of the simulation study presented in Section 4.1 by estimating the quantities $\lambda^{11}(t)$ from the screening arm of the trial. For the observable transition intensities, we calculated smooth hazard estimates using the muhaz (Hess and Gentleman, 2019) package in R and integrated them numerically to get the corresponding cumulative intensities.

To derive a connection between the observable and unobservable quantities under the identifying assumptions listed in Section 2.3, we can re-express the cumulative incidence of cancer-specific mortality in the control arm of a conventional trial $E[N_{13}^{0}(t)]$ as the sum of two different alternative event histories in terms of the latent subgroup membership of the individual, that is, the event history under early detection and no early treatment and the event history under no early detection by time $t$. This gives

	$\displaystyle E[N_{13}^{0}(t)]$	$\displaystyle=\int_{0}^{t}\exp\left\{-(\Lambda_{13}^{0}(u)+\Lambda_{14}^{0}(u))\right\}\lambda_{13}^{0}(u)\,\textrm{d}u$
		$\displaystyle=E[N_{13}^{0}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
		$\displaystyle\quad+E[N_{13}^{0}(t)\mid N_{12}^{1}(t)=0]P(N_{12}^{1}(t)=0)$
		$\displaystyle\stackrel{{\scriptstyle\textrm{(i),(ii)}}}{{=}}E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
		$\displaystyle\quad+E[N_{13}^{1}(t)\mid N_{12}^{1}(t)=0]P(N_{12}^{1}(t)=0)$
		$\displaystyle=\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{12}^{1}(u)$
		$\displaystyle\quad\times\int_{u}^{t}\exp\left\{-\left(\int_{u}^{s}[\lambda_{23}^{10}(s)+\lambda_{24}^{10}(s)]\,\textrm{d}s\right)\right\}\lambda_{23}^{10}(v)\,\textrm{d}v\,\textrm{d}u$
		$\displaystyle\quad+\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{13}^{1}(u)\,\textrm{d}u$

Here the second equality used assumptions (i) and (ii). Finally, applying assumption (iii), as well as the modeling assumption of proportional hazards, we get

	$\displaystyle E[N_{13}^{0}(t)]$	$\displaystyle\stackrel{{\scriptstyle\textrm{(iii)}}}{{=}}\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{12}^{1}(u)$
		$\displaystyle\quad\times\int_{u}^{t}\exp\left\{-\left(\int_{u}^{v}[\theta\lambda_{23}^{1}(s)+\lambda_{24}^{1}(s)]\,\textrm{d}s\right)\right\}\theta\lambda_{23}^{1}(v)\,\textrm{d}v\,\textrm{d}u$
		$\displaystyle\quad+\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{13}^{1}(u)\,\textrm{d}u$		(5)

A corresponding equation using time since the early diagnosis as the time scale after transition to state 2 would be obtained by replacing $\lambda_{2l}^{1}(v)$ with $\lambda_{2l}^{1}(v-u)$, $l\in\{3,4\}$. A related equation was considered by McIntosh (1999), but without considering the hazard ratio as the parameter of interest. We note that in the case of only single baseline screening examination at time $t=0$ and the absence of competing risks, (3.1) would simplify to

	$\displaystyle E[N_{13}^{0}(t)]$	$\displaystyle=E[N_{23}^{10}(t)\mid N_{12}^{1}(0)=1]P(N_{12}^{1}(0)=1)$
		$\displaystyle\quad+E[N_{13}^{1}(t)\mid N_{12}^{1}(0)=0]P(N_{12}^{1}(0)=0)$
		$\displaystyle=P(N_{12}^{1}(0)=1)\int_{0}^{t}\exp\left\{-\theta\Lambda_{23}^{1}(u)\right\}\theta\lambda_{23}^{1}(u)\,\textrm{d}u$
		$\displaystyle\quad+P(N_{12}^{1}(0)=0)\int_{0}^{t}\exp\left\{-\Lambda_{13}^{1}(u)\right\}\lambda_{13}^{1}(u)\,\textrm{d}u,$

which corresponds to the equality given by Martinussen et al. (2019, p. 69).

In addition to the parameter of interest $\theta$, (3.1) involves other unknown quantities, needing to be estimated. We note that under counterfactual consistency (iv) all the quantities (except $\theta$) on the right hand side of (3.1) can be estimated from the screening arm of the conventional trial, while the left hand size can be estimated from the control arm of the conventional trial. In principle, any associational model can be used to obtain estimators to substitute in for the quantities on the right hand side. For example, proportional cause-specific hazard models where the mortality after the early diagnosis is allowed to depend on time $u$ of the early diagnosis can be fitted for $\lambda_{2l}^{1}(v)$, $l\in\{3,4\}$. However, here we focus on the special case of Markov multi-state model, where $\lambda_{2l}^{1}(v)$ does not depend on $u$, in which case we can use non-parametric estimators to substitute in for the unknown quantities.

State occupation probabilities in a Markov multi-state model can be generally estimated using the non-parametric Aalen-Johansen approach (Borgan, 2014), which in the case of a competing risks setting reduces to the non-parametric cumulative incidence estimator for the left hand side of the equation. The required inputs for the right hand side are Nelson-Aalen/Breslow estimates for the cumulative hazards/hazard increments; using these as inputs (with one of them modified by the multiplicative factor $\theta$), the Aalen-Johansen state occupation probabilities can be calculated for example using the msfit and probtrans functions of the mstate package (Putter et al., 2007; de Wreede et al., 2011). Equation (3.1) can then be solved numerically with respect to $\log\theta$ to estimate the early treatment effect.

Equation (3.1) was written for the underlying counting processes in the absence of censoring. However, by substituting in estimators that can accommodate independent censoring (unconditionally on covariates), the proposed approach can accommodate independent censoring. We discuss in Section 6 how covariate-dependent censoring can be accommodated. The fixed time $t$ at which the estimating equation is evaluated can be chosen as the maximum follow-up length, or, as we demonstrate in Section 5, the effect can be estimated as a function of $t$ to see if the estimates stabilize over the follow-up period. Alternatively, if reporting a single estimate is desirable, the timepoint for the estimation can be chosen to minimize the empirical variance as $\hat{\theta}=\arg\min_{t}\hat{V}(\hat{\theta}_{t})$ as suggested by Martinussen et al. (2019), or as the inverse variance weighted average $\hat{\theta}=\left(\sum_{t}\hat{\theta}_{t}/\hat{V}(\hat{\theta}_{t})\right)/\left(\sum_{t}1/\hat{V}(\hat{\theta}_{t})\right)$ over an equally spaced grid of times $t$. The variance estimates can be obtained through the bootstrap.

We note that we can obtain an analogous equation

	$\displaystyle E[N_{14}^{0}(t)]$	$\displaystyle=\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{12}^{1}(u)$
		$\displaystyle\quad\times\int_{u}^{t}\exp\left\{-\left(\int_{u}^{v}[\theta\lambda_{23}^{1}(s)+\lambda_{24}^{1}(s)]\,\textrm{d}s\right)\right\}\lambda_{24}^{1}(v)\,\textrm{d}v\,\textrm{d}u$
		$\displaystyle\quad+\int_{0}^{t}\exp\left\{-(\Lambda_{12}^{1}(u)+\Lambda_{13}^{1}(u)+\Lambda_{14}^{1}(u))\right\}\lambda_{14}^{1}(u)\,\textrm{d}u.$		(6)

for other-cause mortality. However, this is not needed for estimation as we have only one unknown quantity to solve under the assumption that early treatments do not have effect on other-cause mortality. We note that the proposed methodology can be extended to relax this assumption, that is $\lambda^{10}_{24}(t)=\lambda^{11}_{24}(t)$. In that case, we can solve the two equations for two unknowns, characterizing the effect on cancer-specific and other cause mortality, respectively.

Alternatively to the estimating equation approach, we can consider a multinomial likelihood expression for the parameter $\theta$ using data from the mortality outcomes in the conventional trial. This can be written as

	$\displaystyle L(\theta)$	$\displaystyle=\prod_{\{i:Z_{i}=1\}}\Big{\{}{E[N_{13}^{1}(t)+N_{23}^{1}(t)]}^{N_{13}^{1}(t)+N_{23}^{1}(t)}{E[N_{14}^{1}(t)+N_{24}^{1}(t)]}^{N_{14}^{1}(t)+N_{24}^{1}(t)}$
		$\displaystyle\qquad\times\left(1-E[N_{13}^{1}(t)+N_{23}^{1}(t)]-E[N_{14}^{1}(t)+N_{24}^{1}(t)]\right)^{1-N_{13}^{1}(t)-N_{23}^{1}(t)-N_{14}^{1}(t)-N_{24}^{1}(t)}\Big{\}}$
		$\displaystyle\quad\times\prod_{\{i:Z_{i}=0\}}E[N_{13}^{0}(t)]^{N_{13}^{0}(t)}E[N_{14}^{0}(t)]^{N_{14}^{0}(t)}$
		$\displaystyle\qquad\times\left(1-E[N_{13}^{0}(t)]-E[N_{14}^{0}(t)]\right)^{(1-N_{13}^{0}(t)-N_{14}^{0}(t))}$
		$\displaystyle\stackrel{{\scriptstyle\theta}}{{\propto}}\prod_{\{i:Z_{i}=0\}}E[N_{13}^{0}(t)]^{N_{13}^{0}(t)}E[N_{14}^{0}(t)]^{N_{14}^{0}(t)}$
		$\displaystyle\qquad\times\left(1-E[N_{13}^{0}(t)]-E[N_{14}^{0}(t)]\right)^{1-N_{13}^{0}(t)-N_{14}^{0}(t)},$		(7)

where $E[N_{13}^{0}(t)]$ and $E[N_{14}^{0}(t)]$ are as in (3.1) and (3.1). The proportionality followed because, considering the transition intensities in the screening arm as fixed quantities, the likelihood contributions of the individuals in the screening arm do not depend on $\theta$. Similarly to the estimating equation in Section 3.1, (3.2) involves unknown quantities that under the stated assumptions can be replaced with estimators from the screening arm of the trial, and the expression $L(\theta)$ then maximized numerically with respect to $\log\theta$. The limitation of the multinomial likelihood expression (3.2) is that it can only accommodate type I censoring. However, instead of evaluating the likelihood expression at a fixed time $t$, independent censoring can be accommodated by evaluating the likelihood contributions at the minimum of the time of death or censoring time. Such a likelihood expression could also be informative of time-dependent treatment effects. We will briefly outline this in Section 6.

To compare the estimators proposed in Section 3 in a setting resembling a real randomized screening trial, we generated data using the algorithm described in Section 2.4, using the NLST as the basis of the simulation study. In particular, because the NLST was powered for the ITS effect, we were interested in determining if such a study is also powered for detecting the subgroup effect using tests based on the proposed estimators. For the purpose of comparison, we also included tests based on these to the two estimators proposed earlier for absolute and proportional subgroup mortality risk reduction (Saha et al., 2018). In the present notational framework, the absolute mortality risk reduction (3) can be expressed as

	$\displaystyle E[N_{13}^{0}(t)-\{N_{13}^{1}(t)+N_{23}^{1}(t)\}]$
	$\displaystyle=E[N_{13}^{0}(t)-\{N_{13}^{1}(t)+N_{23}^{1}(t)\}\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
	$\displaystyle\quad+E[N_{13}^{0}(t)-\{N_{13}^{1}(t)+N_{23}^{1}(t)\}\mid N_{12}^{1}(t)=0]P(N_{12}^{1}(t)=0)$
	$\displaystyle\stackrel{{\scriptstyle\textrm{(i),(ii),(iii)}}}{{=}}E[N_{23}^{10}(t)-\{N_{13}^{1}(t)+N_{23}^{11}(t)\}\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
	$\displaystyle=E[N_{23}^{10}(t)-N_{23}^{11}(t)\mid N_{12}^{1}(t)=1]E[N_{12}^{1}(t)]$

Thus, under (iv) we get,

$\displaystyle\eqref{equation:acfr}=\frac{E[N_{13}(t)\mid Z=0]-E[N_{13}(t)+N_{23}(t)\mid Z=1]}{E[N_{12}(t)\mid Z=1]},$

(8)

where the expectations can be estimated non-parametrically as cumulative incidences. Further, we can write

	$\displaystyle E[N_{13}^{0}(t)]$
	$\displaystyle=E[N_{13}^{0}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
	$\displaystyle\quad+E[N_{13}^{0}(t)\mid N_{12}^{1}(t)=0]P(N_{12}^{1}(t)=0)$
	$\displaystyle\stackrel{{\scriptstyle\textrm{(i),(ii)}}}{{=}}E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)$
	$\displaystyle\quad\quad+E[N_{13}^{1}(t)\mid N_{12}^{1}(t)=0]P(N_{12}^{1}(t)=0)$
	$\displaystyle=E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)+P(N_{13}^{1}(t)=1,N_{12}^{1}(t)=0)$
	$\displaystyle=E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)+E[N^{1}_{13}(t)]$

Using the above two results, we can now re-express (4) as,

	$\displaystyle 1-\frac{E[N_{23}^{11}(t)\mid N_{12}^{1}(t)=1]}{E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]}$
	$\displaystyle=\frac{(E[N_{23}^{10}(t)\mid N_{12}^{1}(t)=1]-E[N_{23}^{11}(t)\mid N_{12}^{1}(t)=1])P(N_{12}^{1}(t)=1)}{E[N_{23}^{10}(t)=1\mid N_{12}^{1}(t)=1]P(N_{12}^{1}(t)=1)}$
	$\displaystyle=\frac{E[N^{0}_{13}(t)]-E[N^{1}_{13}(t)+N^{1}_{23}(t)]}{E[N_{13}^{0}(t)]-E[N^{1}_{13}(t)]}$

Finally, under (iv) we can write that

$\displaystyle\eqref{equation:pcfr}=\frac{E[N_{13}(t)\mid Z=0]-E[N_{13}(t)+N_{23}(t)\mid Z=1]}{E[N_{13}(t)\mid Z=0]-E[N_{13}(t)\mid Z=1]},$

(9)

where again the expectations can be replaced with appropriate non-parametric cumulative incidence estimators. While (8) and (9) are estimating mortality risk reductions rather than the hazard ratio, the four estimators are comparable in terms of the power of a Wald-type or confidence interval based test on them. This will help answer the question of whether conventional screening trial such as NLST is powered to detect the subgroup early treatment effect. For this reason, we also calculated estimates for the ITS mortality reductions as

$$\eqref{equation:itt}=E[N_{13}(t)\mid Z=0]-E[N_{13}(t)+N_{23}(t)\mid Z=1]$$

(10)

and the proportional reduction by

$$\eqref{equation:ittp}=1-\frac{E[N_{13}(t)+N_{23}(t)\mid Z=1]}{E[N_{13}(t)\mid Z=0]},$$

(11)

as the trial is powered for these.

To compare the four estimators, we generated $1000$ NLST-like datasets, by fixing the transition intensities to those estimated from the NLST, and fixing the $\log$ hazard ratio to $0.4804$, estimated from the real data using the estimating equation (3.1). While we could have fixed the hazard ratio to any value in the simulation, we used the empirical value to obtain similar death counts to the real data. For each dataset we simulated $n=53,452$ individuals assigned randomly with equal probability to screening or control. Since in the NLST the censoring is mostly administrative at the end of the follow-up at 7 years, we did not simulate random censoring, so the simulated datasets had only type I censoring at 7 years. The true values of the subgroup mortality reductions were evaluated by calculating the quantities (3) and (4) under the specified transition intensities, that is, evaluating the cumulative incidences through numerical integration using the true transition intensity functions used in the simulation. In the estimators, the cumulative incidences were calculated using the cuminc function of the cmprsk R package (Gray, 2019). For more complex state occupation probabilities needed for evaluating the estimating equation and the likelihood expression, we first calculated Nelson-Aalen/Breslow estimators, which were used as inputs to Aalen-Johansen estimators for state occupation probabilities using the functions available in the mstate R package. The estimating equation was solved numerically using the uniroot function and the likelihood maximized numerically using the optim function of R (R Core Team, 2017).

The mortality reductions, as well as the estimating equation (3.1) and the likelihood expression (3.2) were evaluated at 7 years. To estimate standard errors for each estimator, we bootstrapped each simulated dataset 50 times, and took the standard deviation of the point estimates as the standard error. This was compared to the Monte Carlo standard deviation of the point estimates over the simulation rounds. 95% confidence intervals were constructed using the normal approximation and compared to the nominal coverage probability. The power of the confidence interval based test was calculated as the proportion of simulation rounds where the interval did not cover the null value.

The simulation results are presented in Table 1. These indicate that all the estimators could capture the corresponding true value relatively well, based on low bias and near nominal coverage probability. The power of the subgroup treatment effect estimators was also comparable to the ITS analysis. This is consistent with the similar property of g-estimation methods (White, 2005). This can be explained by the fact that all the mortality reduction due to screening is realized through the effect of the early treatments, and thus focusing on the subgroup does not lose any of the effect. The test based on the estimator (9) had the highest power. The two hazard ratio estimators gave very similar results, so we only use the estimating equation for the data analysis in the next section.

To check the performance of the proposed methods under settings not resembling a typical randomized screening trial with large size and rare events, we also simulated scenarios with smaller sample size and relatively more common events. In particular, we simulated scenarios with more than $50\%$ in the screening arm early diagnosed through (compared to $2.4\%$ in NLST) and among them majority experiencing cancer-specific death. The screening arm state transitions were simulated from the following constant transition rates, with follow-up type I censored at $t=7$:

$$\lambda^{11}(t)=\blockarray{c@{\hskip 1.0pt}rrrr@{\hskip 4.0pt}}&1234\\ \block{r@{\hskip 1.0pt}|@{\hskip 1.0pt}|@{\hskip 1.0pt}rrrr@{\hskip 1.0pt}|@{\hskip 1.0pt}|}1.0.22800.11480.0168\\ 20.0.19800.0111\\ 300.0\\ 4000.\\$$

For the control arm, all the intensities in the matrix were the same, except for $\lambda^{10}_{23}(t)=0.1980\,\times\theta$, where $\log\theta$ was fixed to $0.4700$ (i.e. $\theta=1.6$). From the intensity matrices, we generated screening trial data using the algorithm of Section 2.4 with sample sizes $n=500,800,$ and $1000$ and estimated the hazard ratios in the screening detectable subgroup using the two methods described in Section 3. We did not consider the risk reduction type measures here. The results for a scenario without unmeasured confounding are presented in Table 2 (rows with $\beta=0$).

Results in Table 2 indicate that the two estimators again perform very similarly. With $n=500$, they show some small sample bias, which disappears with increasing sample size. The bootstrap standard errors give a reasonable approximation of the sampling distribution standard deviation. These results indicate that when the modeling and identifying assumptions are satisfied, the methods perform reasonably also with smaller sample size and more common events.

To study how the methods perform in the presence of unmeasured confounding, we modified the above transition intensities for transitions $1\rightarrow 2$, $1\rightarrow 3$ and $2\rightarrow 3$ as $\lambda_{12}^{1}(t)\exp(\beta U)$, $\lambda_{13}^{1}(t)\exp(\beta U)$ and $\lambda_{23}^{11}(t)\exp(\beta U)$, where $U\sim\textrm{Bernoulli}(0.5)$ is an unmeasured baseline covariate. This represents a scenario where some individuals are at higher risk for cancer diagnosis and cancer mortality, but under which the instrumental variable assumptions are still satisfied. However, because the data are simulated from a conditional generating mechanism, some degree of model misspecification is expected in the estimators. In addition to the proportionality of the marginal effect, the estimators use a Markov multi-state model for calculation of the state occupancy probabilities. While the Markov assumption can be relaxed, here we study its effect on the estimates. The true marginal effects were approximated as explained in Section 2.4.

The results with $\beta=0.34$ and $\beta=0.47$ are also summarized in Table 2. They show as expected that the true marginal effects are now smaller compared to the conditional effect used in the simulation. The points estimates were close to the true marginal effects under most of the scenarios and sample sizes, with coverage probabilities close to nominal. This suggests that the estimates are free of confounding bias, as is expected of instrumental variable estimators.