Estimating and testing sequential causal effects based on alternative G-formula: an observational study of the influence of early diagnosis on survival of cardia cancer

Abstract

Cancer diagnosis is part of a complex stochastic process, in which patients' personal and social characteristics influence the choice of diagnosing methods, diagnosing methods in turn influence the initial assessment of cancer stage, cancer stage in turn influences the choice of treating methods, and treating methods in turn influence cancer outcomes such as cancer survival. To evaluate the performance of diagnoses, one needs to estimate and test the sequential causal effect (SCE) under a specified regime of diagnoses and treatments in such a complex observational study, where the data-generating mechanism is unknown and modeling is needed for statistical inference. In this article, we introduce a method of statistical modeling to estimate and test SCEs under regimes of treatments (diagnoses and treatments in cancer diagnosis) in complex observational studies. By applying the alternative G-formula, we express the SCE in terms of the point effects of treatments in the sequence, so that the modeling can be conducted via the point effects in the framework of single-point causal inference. We illustrate our method by a medical example of cancer diagnosis with data from a Swedish prognosis study of cardia cancer.

KEYWORDS:

• Cancer diagnosis
• G-formula
• Point effect
• Sequential causal effect
• Statistical modeling

1. Introduction

Cancer diagnosis is part of a complex stochastic process, in which patients' personal and social characteristics influence the choice of diagnosing methods, diagnosing methods in turn influence the initial assessment of cancer stage, cancer stage in turn influences the choice of treating methods, and treating methods in turn influence cancer outcomes such as cancer survival. To evaluate the performance of cancer diagnoses, one needs to estimate and test the sequential causal effect (SCE), which is the causal effect under a specified regime of diagnoses and treatments on a certain outcome of interest. An important example of SCEs is the one under optimal regime.

In sequential causal inference, Robins (1986, 1997, 2009) derived the well-known G-formula, which identifies the SCE under a specified regime of treatments (diagnoses and treatments in cancer diagnosis) via standard parameters. Based on the G-formula, a parametric method has been developed for estimating the SCE via standard parameters (Hernan and Robins 2020; Robins 1997, 2009; Taubman et al. 2009). However, this method may suffer from the curse of dimensionality and the null paradox. See also the improved versions of the parametric method. Still, these versions have the difficulty of specifying the influence of the time-dependent covariates due to the curse of dimensionality and the null paradox (Almirall, Have, and Murphy 2010; Henderson, Ansell, and Alshibani 2010; Chaffe and van der Laan 2012).

To avoid these problems, the literature focuses on semi-parametric methods. One class of methods is the marginal structural model based on the inverse probability of treatment weighting (Hernan and Robins 2020; Hernan, Brumback, and Robins 2000; Robins, Brumback, and Hernan 2000; Robins 2009) or the augmented the inverse probability of treatment weighting (Zhang et al. 2013). Another class of methods is the G-estimation based on the structural nested mean model (SNMM) or optimal SNMM (Hernan and Robins 2020; Robins 1997, 2004, 2009); the G-estimation incorporates the A-learning (for instance, Murphy 2003) and Q-learning (for instance, Sutton and Barto 1998). See also a summary of the development of these methods in the area of optimal regimes (Kosorok et al. 2021). These methods achieve consistent estimates of SCEs. In a recent review, Kosorok and Laber (2019) highlighted the need for methods of evaluating the uncertainty in estimating SCEs in the context of optimal regimes (Sec. 6, statistical inference).

Recently, Wang and Yin (2015, 2020) derived an alternative G-formula, which expresses the SCE in terms of the point effects of treatments in the sequence. The point effect is simply the effect of single-point treatment in single-point causal inference (Rosenbaum and Rubin 1983). In designed experiments, where data-generating mechanisms are known, they estimated the SCE by estimating point effects of treatments without the need for modeling. In observational studies such as the influence of early cancer diagnosis, where data-generating mechanism is unknown, one needs to conduct modeling to estimate and test the SCE.

In this article, we will introduce a method of statistical modeling to estimate and test the SCE in observational studies and illustrate our method by a medical example using data from a Swedish prognosis study of cardia cancer. In Sec. 2, we describe our motivating example. In Sec. 3, we introduce the general framework of our method. We apply the alternative G-formula to express the SCE in terms of point effects of treatments in the sequence and describe how to conduct modeling to estimate and test the point effects and then the SCE in the framework of single-point causal inference. In Sec. 4, we apply our method to estimate and test the SCE under any regime of diagnosing and treating hospitals on one-year survival in our medical example. In Sec. 5, we compare our method with the three available methods mentioned earlier by applying all these four methods to the data. Finally, in Sec. 6, we conclude the article with discussion.

2. Data and the medical background

In Sweden, patients usually seek medical help at hospitals near their residential areas, namely, catchment areas. When cancer is diagnosed, they may stay at the diagnosing hospital or transfer to another hospital for treatment. The hospital diagnosing cancer is called the diagnosing hospital, while the one treating cancer is called the treating hospital. To evaluate the performance of diagnosing and treating hospitals, we may study cancer outcomes such as survival under various regimes of diagnosing and treating hospitals among cancer patients after adjusting for patients' differences.

The data used in this study is from a prognosis study conducted during a period between 1988 and 1995 in hospitals in central and northern Sweden (Hansson et al. 2000). It contained information of 157 patients of cardia cancer. Cardia cancer is highly malignant with bad prognosis, and one-year survival is a good measure of the performance of both diagnosing and treating hospitals. A question of relevance to the public health policy is which types of the diagnosing and treating hospitals, large versus small, perform better on cancer outcomes, where the large type refers to the regional or county hospitals and the small type to local hospitals.

In Sec. 4 below, we will estimate and test the SCE on one-year survival under any regime of diagnosing and treating hospitals and then use the SCEs to evaluate the performance of diagnosing and treating hospitals.

3. General framework of estimating and testing SCEs based on the alternative G-formula

In this section, we will describe the alternative G-formula, which identifies the SCE by the point effects of treatments, and the general framework of our method for estimating and testing SCEs in observational studies.

3.1. Potential variables and the causal effect from a sequence of treatments

Let $D_t$ at time $t=1,\mathrm{ }\dots ,\mathrm{ }T$ be the treatment variable, which potentially and deterministically assigned treatments $z_t$ to each unit of the population. A regime is a sequence of such treatment variables, ${\mathit{D}}_1^T=(D_1,\dots ,\mathrm{ }{D}_T).$ Prior to $D_1,$ there exists a set of stationary covariates ${\mathit{X}}_1.$ Under ${\mathit{D}}_1^T,$ each unit could have a set of potential time-dependent covariates ${\mathit{X}}_t({\mathit{D}}_1^{t-1})$ between $D_{t-1}$ and $D_t$ ($t=2,\mathrm{ }\dots ,\mathrm{ }T$) and a potential outcome $Y({\mathit{D}}_1^T)$ of interest after the last treatment $D_T.$ Noticeably, $D_t$ can be dynamic, namely, the assignment of the treatment depends on the earlier treatments and covariates of each unit.

Let ${\mathit{X}}_1^T\left({\mathit{D}}_1^{T-1}\right)=\{{\mathit{X}}_1,{\mathit{X}}_2(D_1),\dots ,\mathrm{ }{\mathit{X}}_T({\mathit{D}}_1^{T-1})\}$ be the set of all stationary covariates and potential covariates, and {${\mathit{X}}_1^T\left({\mathit{D}}_1^{T-1}\right),{\mathit{D}}_1^T,Y({\mathit{D}}_1^T)\}=\{{\mathit{X}}_1,D_1,{\mathit{X}}_2\left(D_1\right),D_2,\dots ,{\mathit{X}}_T\left({\mathit{D}}_1^{T-1}\right),D_T,Y({\mathit{D}}_1^T)\}$ be the set of all stationary and potential covariates, treatment variables, and the potential outcome. The realization of these variables is $\{{\mathit{x}}_1^T,{\mathit{z}}_1^T,y\}=\{{\mathit{x}}_1,z_1,{\mathit{x}}_2,z_2,\dots ,{\mathit{x}}_T,z_T,y\}.$ In stratum $({\mathit{x}}_1^t,{\mathit{z}}_1^{t-1}),$ we document the regime as ${\mathit{D}}_t^T=(D_t,\dots ,\mathrm{ }{D}_T),$ its potential covariates as ${\mathit{X}}_{t+1}^T\left({\mathit{D}}_t^{T-1}\right)$ and potential outcome as $Y({\mathit{D}}_t^T).$

Let ${\mathit{D}}_t^T={\mathbf{0}}_t^T=(0,\mathrm{ }0,\mathrm{ }\dots ,\mathrm{ }0)\mathrm{ }$ be the control regime given $({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}).$ The blip effect $\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)$ of treatment $z_t$ in stratum $({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1})$ is an increase of the mean of potential outcome $Y({\mathit{D}}_t^T)$ under regimes ${\mathit{D}}_t^T=(z_t,0_{t+1}^T)$ versus ${\mathit{D}}_t^T=0_t^T$ in the stratum,$\mathrm{ }$ that is $$\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)=E\left\{Y\left(z_t,0_{t+1}^T\right)\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}\right\}-E\left\{Y\left(0_t^T\right)\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}\right\}.$$The blip effect is called the net effect (Wang and Yin 2015). Clearly, $\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t=0\right)=0.$ The SCE under regime ${\mathit{D}}_t^T$ in stratum $({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}),$ ${\text{SCE}\mathrm{(}{\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }\mathit{D}}_t^T),\mathrm{ }$ is an increase of the mean of potential outcome $Y\left({\mathit{D}}_t^T\right)$ under ${\mathit{D}}_t^T$ versus $Y\left({\mathbf{0}}_t^{\mathit{T}}\right)$ under ${\mathbf{0}}_t^T$ in the stratum, namely, $$\text{SCE}\left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{\mathit{D}}_t^T\right)=E\left\{Y\left({\mathit{D}}_t^{\mathit{T}}\right)\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}\right\}-E\left\{Y\left({\mathbf{0}}_t^{\mathit{T}}\right)\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}\right\}.$$Here we see that if ${\mathit{D}}_t^T=(z_t,{\mathbf{0}}_{t+1}^T),$ then $\text{SCE}\left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{\mathit{D}}_t^T\right)=\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right).\mathrm{ }$ Therefore the blip effect is a special type of SCEs. As it will become clear in Sec. 3.4, the blip effects determine the SCE under any regime of treatments.

The SCE under regime ${\mathit{D}}_1^T$ in the population, ${\text{SCE}\mathrm{(}\mathit{D}}_1^T),$ is an increase of the mean of potential outcome $Y\left({\mathit{D}}_1^T\right)$ under ${\mathit{D}}_1^T$ versus $Y\left({\mathbf{0}}_1^{\mathit{T}}\right)$ under ${\mathbf{0}}_1^T$ in the population, namely, $${\text{SCE}\mathrm{(}\mathit{D}}_1^T)=E\{Y\left({\mathit{D}}_1^{\mathit{T}}\right)\}-E\{Y\left({\mathbf{0}}_1^{\mathit{T}}\right)\}$$

3.2. Observable variables and the point effects of treatments in the sequence

Now we consider observable variables. Let $Z_t$ be the treatment variable at $t=1,\mathrm{ }\dots ,\mathrm{ }T;$ ${\mathit{X}}_1\mathrm{ }$ be the set of stationary covariates, ${\mathit{X}}_t$ be the set of time-dependent covariates between $Z_{t-1}$ and $Z_t;$ and $Y$ be the outcome of interest. Let ${\mathit{Z}}_1^T=\{Z_1,\dots ,\mathrm{ }{Z}_T\}$ be the sequence of treatment variables; ${\mathit{X}}_1^T=\{{\mathit{X}}_1,\dots ,\mathrm{ }{\mathit{X}}_T\}$ be the set of all stationary and time-dependent covariates; $\left\{{\mathit{X}}_1^T,{\mathit{Z}}_1^T,\mathrm{ }Y\right\}=\{{\mathit{X}}_1,Z_1,{\mathit{X}}_2,Z_2,\dots ,\mathrm{ }{\mathit{X}}_T,Z_T,Y\}$ be the set of all treatments, covariates and outcome.

Suppose that the observable variables $\left\{{\mathit{X}}_1^T,{\mathit{Z}}_1^T,\mathrm{ }Y\right\}$ have the same support as the potential ones ${\{\mathit{X}}_1^T\left({\mathit{D}}_1^{T-1}\right),{\mathit{D}}_1^T,Y({\mathit{D}}_1^T)\}.$ The observed values of these variables are $\left\{{\mathit{x}}_1^T,{\mathit{z}}_1^T,y\right\}=\{{\mathit{x}}_1,z_1,{\mathit{x}}_2,z_2,\dots ,\mathrm{ }{\mathit{x}}_T,z_T,y\}.$ In the following, we will use $P(.)$ to denote the probability distribution of a discrete variable or the density distribution of a continuous variable. The joint distribution of $\left\{{\mathit{X}}_1^T,{\mathit{Z}}_1^T,\mathrm{ }Y\right\}$ is given by (1) $$P\left({\mathit{x}}_1^T,{\mathit{z}}_1^T,y\right)=P\left({\mathit{x}}_1\right)P{z}_1\mathrm{ }\left|\mathrm{ }{\mathit{x}}_1\right)\cdots P{\mathit{x}}_T\mathrm{ }|{\mathrm{ }\mathit{x}}_1^{T-1},{\mathit{z}}_1^{T-1})P\left(z_T\mathrm{ }\right|{\mathrm{ }\mathit{x}}_1^T,{\mathit{z}}_1^{T-1}P\left(y\mathrm{ }{\mathrm{ }\mathit{x}}_1^T,{\mathit{z}}_1^T\right)\mathrm{ }$$(1) The standard parameters for these conditional distributions are simply their conditional means, for instance, the standard parameter for $Py\mathrm{ }|\mathrm{ }{\mathit{x}}_1^T,{\mathit{z}}_1^T)$ is the conditional mean $\mu \left({\mathit{x}}_1^T,{\mathit{z}}_1^T\right)=E(Y\mathrm{ }|{\mathrm{ }\mathit{x}}_1^T,{\mathit{z}}_1^T).$

Now we describe the point effects of treatments in the sequence, which are also parameters for the outcome (Wang and Yin 2015). Let $\mu \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},\mathrm{ }{z}_t\right)=E\left(Y\mathrm{ }\right|\mathrm{ }{\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},\mathrm{ }{z}_t)$ be the conditional mean of observable outcome $Y$ given $({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},\mathrm{ }{z}_t).$ Without loss of generality, we take $z_t=0$ as the control treatment. Then the point effect of treatment $z_t$ in the stratum ${(\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1})\mathrm{ }$ is the following difference between the means (2) $$\vartheta \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)=\mu \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},\mathrm{ }{z}_t\right)-\mu \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},\mathrm{ }0\right)$$(2) Clearly, $\vartheta \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t=0\right)=0.$ The point effect is simply the point effect of treatment in single-point causal inference and its estimation and hypothesis testing is well studied (for instance, Rosenbaum and Rubin 1983).

To identify the SCEs by the observable variables $\left\{{\mathit{X}}_1^T,{\mathit{Z}}_1^T,\mathrm{ }Y\right\},$ Robins (1986, 1997, 2009) introduced the identifying condition, which consists of the consistency assumption, the positivity assumption for treatments and the assumption of no unmeasured confounders. The identifying condition is satisfied in sequential randomized experiments where treatment $z_t$ is randomly assigned according to a history ${(\mathit{x}}_1^{\mathit{t}},\mathrm{ }{\mathit{z}}_1^{t-1}$) of the earlier covariates and treatments. It is approximately satisfied in observational studies with a sufficient number of covariates. Throughout, we assume the identifying condition.

Under the identifying condition, Robins (1986, 1997, 2009) derived the well-known G-formula, which expresses the SCE in terms of the standard parameters $\mu ({\mathit{x}}_1^T,{\mathit{z}}_1^T).$ If ${\mathit{x}}_t$ ($t>1$) has influence from earlier treatments ${\mathit{Z}}_1^{t-1}$ and influence on subsequent treatments ${\mathit{Z}}_t^T,$ then $\mu ({\mathit{x}}_1^T,{\mathit{z}}_1^T)$ are essentially all different for different $({\mathit{x}}_1^{\mathit{T}},{\mathit{z}}_1^{T-1}),$ even if each $Z_t$ has null effect on $Y.$ So any model imposing equalities between $\mu ({\mathit{x}}_1^T,{\mathit{z}}_1^T)$ for different $({\mathit{x}}_1^T,{\mathit{z}}_1^{T-1})$ is misspecified (null paradox). Obviously, the number of $\mu ({\mathit{z}}_1^T,{\mathit{x}}_1^T)$ is huge in case of large $T$ or even uncountable (curse of dimensionality).

Under the same identifying condition, Wang and Yin (2020) derived an alternative G-formula, which expresses the blip effect and SCE in terms of point effects instead of the standard parameters. In the rest of this section, we will first describe the G-formula, then describe our method of conducting modeling to estimate and test blip effects and SCEs based on the alternative G-formula.

3.3. Alternative G-formula to estimate and test the blip effect in observational studies

Intuitively, the point effect $\vartheta \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)$ is the result of treatment $z_t$ and the subsequent treatments ${\mathit{Z}}_{\mathit{t}+1}^{\mathit{T}},$ and it may decompose into a sum of the blip effects of these treatments. Therefore, we should have the following decomposing formula (3) $$\begin{array}{c}\vartheta \left({\mathit{x}}_1^t, {\mathit{z}}_1^{t-1}; z_t\right)=\varphi \left({\mathit{x}}_1^t, {\mathit{z}}_1^{t-1}; z_t\right)+\\ \sum _{s=t+1}^{T}E\{\varphi \left({\mathit{x}}_1^s, {\mathit{z}}_1^{s-1}; z_s\right) | {\mathit{x}}_1^t, {\mathit{z}}_1^{t-1},z_t\}-\sum _{s=t+1}^{T}E\left\{\varphi \left({\mathit{x}}_1^s, {\mathit{z}}_1^{s-1}; z_s\right) \right| {\mathit{x}}_1^t, {\mathit{z}}_1^{t-1},z_t=0\}\end{array}$$(3) where the expectations are with respect to the conditional distributions of observable variables $P({\mathit{x}}_{t+1}^s,\mathrm{ }{\mathit{z}}_{\mathit{t}+1}^{s-1},\mathrm{ }{z}_s\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},z_t)$ and $P({\mathit{x}}_{t+1}^s,\mathrm{ }{\mathit{z}}_{\mathit{t}+1}^{s-1},\mathrm{ }{z}_s\mathrm{ }|\mathrm{ }{\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},z_t=0).$ Wang and Yin (2020) theoretically proved (3) under the identifying condition. Formula (3) is the alternative G-formula and expresses the point effects in terms of blip effects.

Now we describe our method of conducting modeling to estimate and test blip effects. First, the point effect is simply the point effect of treatment in single-point causal inference, so we can estimate it by modeling the mean $\mu {(\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^t)=EY\mathrm{ }|{\mathrm{ }\mathit{x}}_1^t,{\mathit{z}}_1^t)$ in the framework of single-point causal inference. As an illustration, if ${(\mathit{x}}_1^{t-1},\mathrm{ }{\mathit{z}}_1^{t-1})$ influences the outcome Y via ${\mathit{x}}_t,$ then we have $\mu {(\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^t)=EY\mathrm{ }|{\mathrm{ }\mathit{x}}_t,z_t),$ denoted by $\mu \left({\mathit{x}}_t,z_t\right),$ i.e., $\mu {(\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^t)=\mu ({\mathit{x}}_t,z_t).$ Thus, we may rewrite (2) as $$\vartheta \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};{\mathrm{ }z}_t\right)=\vartheta \left({\mathit{x}}_t;{\mathrm{ }z}_t\right)=\mu \left({\mathit{x}}_t,z_t\right)-\mu \left({\mathit{x}}_t,z_t=0\right).$$Second, Robins (1997, 2009) and Hernan and Robins (2020) pointed out that the blip effects follow a certain pattern described by structural nested mean model (SNMM), as an illustration, $$\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)=f\left({\mathit{x}}_t,z_t;\mathrm{ }\mathit{\gamma }\right),\mathrm{ }t=1,\dots ,T$$where $f(.)$ is a deterministic function of $({\mathit{x}}_t,z_t)$ and is indexed by a parameter vector $\mathit{\gamma }$ of small dimension. With $\vartheta \left({\mathit{x}}_t;z_t\right)$ and SNMM, the alternative G-formula becomes the following decomposing formula $$\vartheta \left({\mathit{x}}_t;\mathrm{ }{z}_t\right)=f({\mathit{x}}_t,z_t;\mathrm{ }\mathit{\gamma })+\sum _{s=t+1}^{T}E\{f\left({\mathit{x}}_s,z_s;\mathrm{ }\mathit{\gamma }\right)\mathrm{ }|\mathrm{ }{\mathit{x}}_t,z_t\}-\sum _{s=t+1}^{T}E\left\{f\left({\mathit{x}}_s,z_s;\mathrm{ }\mathit{\gamma }\right)\mathrm{ }\right|\mathrm{ }{\mathit{x}}_t,z_t=0\}$$where the expectations are with respect to the conditional distributions of observable variables $P({\mathit{x}}_s,z_s\mathrm{ }|{\mathrm{ }\mathit{x}}_t,z_t)$ and $P({\mathit{x}}_s,z_s\mathrm{ }|{\mathrm{ }\mathit{x}}_t,z_t=0).$

To estimate and test $\mathit{\gamma },$ we treat the above decomposing formula as a regression model, where the response variables are the estimates $\widehat{\vartheta }\left({\mathit{x}}_t;{\mathrm{ }z}_t\right)$ and the explanatory variables are the observed proportions $\widehat{P}({\mathit{x}}_s,z_s\mathrm{ }|{\mathrm{ }\mathit{x}}_t,z_t),$ which are obtained from the data without modeling. In the regression, we need the conditional variance of the estimated point effect given all covariates and treatments $({\mathit{x}}_1^T,\mathrm{ }{\mathit{z}}_1^T),$ which is obtained by specifying a constant dispersion parameter for the distribution $P\left(y\mathrm{ }\right|{\mathrm{ }\mathit{x}}_1^T,{\mathit{z}}_1^T)$. Noticeably, the conditional covariance between the estimated point effects at different times is negligible, for instance, it is equal to zero for the normal outcome (Wang and Yin 2015, 2020). The bootstrap method is used to obtain the covariance matrix $\text{cov}\mathrm{(}\widehat{\mathit{\gamma }}\mathrm{)}$ incorporating the variability of $({\mathit{x}}_1^T,\mathrm{ }{\mathit{z}}_1^T).$

If the estimates for the point effects are consistent and asymptotically normal as typically required in single-point causal inference and if $\mathit{\gamma }$ is of finite dimension, then the estimate $\widehat{\mathit{\gamma }}\mathrm{ }$ is consistent and asymptotically normal. In this case, with $\widehat{\mathit{\gamma }}$ and $\text{cov}\mathrm{(}\widehat{\mathit{\gamma }}\mathrm{)},$ we conduct the Wald test on $\mathit{\gamma }.$

3.4. Alternative G-formula to estimate and test SCEs in observational studies

For the sake of explication, we first describe the alternative G-formula for SCEs in the population. Intuitively, it may be decomposed into a sum of the blip effects of treatments in the regime, as follows. Prior to $D_1,$ the covariate distribution in the population is $P({\mathit{x}}_1),$ so when $D_1$ potentially and deterministically assigned $z_1$ to the population, the contribution of $D_1$ to SCE would be the expectation $E\{\varphi ({\mathit{x}}_1;\mathrm{ }{z}_1)\}$ with respect to $P({\mathit{x}}_1),$ as well-known in single-point causal inference (for instance, Rosenbaum and Rubin 1983). Between $D_1$ and $D_2,$ the covariate distribution in the population under $D_1$ would be $P({\mathit{x}}_1)P{\mathit{x}}_2\mathrm{ }|\mathrm{ }{\mathit{x}}_1,z_1),$ so the contribution of $D_2$ to SCE would be $E\{\varphi ({\mathit{x}}_{1,}{z}_1,{\mathit{x}}_2;\mathrm{ }{z}_2)\}$ with respect to $P({\mathit{x}}_1)P{\mathit{x}}_2\mathrm{ }|{\mathrm{ }\mathit{x}}_1,z_1),$ and so forth. Therefore, we have the following decomposing formula, (4) $$\text{SCE}\left({\mathit{D}}_1^T\right)=\sum _{s=1}^{T}E\left\{\mathrm{\varphi }\left({\mathit{x}}_1^s,\mathrm{ }{\mathit{z}}_1^{s-1};{\mathrm{ }z}_s\right)\right\},\mathrm{ }$$(4) where the expectation is with respect to the distribution of observable covariates ${\mathit{x}}_1^s:$ ${\prod }_{k=1}^{s}P({\mathit{x}}_k\mathrm{ }|{\mathrm{ }\mathit{x}}_1^{k-1},{\mathit{z}}_1^{k-1}).$ Similarly, for the SCE in stratum $({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}),$ we have the decomposing formula (5) $$\text{SCE}\left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{\mathit{D}}_t^T\right)=\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)+\sum _{s=t+1}^{T}E\mathrm{ }\left\{\varphi \left({\mathit{x}}_1^s,\mathrm{ }{\mathit{z}}_1^{s-1};{\mathrm{ }z}_s\right)\mathrm{ }\right|\mathrm{ }{\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1},z_t\},\mathrm{ }$$(5) where the expectation is with respect to the distribution of observable covariates ${\mathit{x}}_{t+1}^s:$ ${\prod }_{k=t+1}^{s}P({\mathit{x}}_k\mathrm{ }|{\mathrm{ }\mathit{x}}_1^{k-1},{\mathit{z}}_1^{k-1}).$ Theoretically, Wang and Yin (2020) proved (4) and (5). The formula (4) or (5) is the alternative G-formula for the SCE and expresses the SCE in terms of the blip effects.

Now we describe our method of conducting modeling to estimate and test SCEs in observational studies. The blip effects $\varphi \left({\mathit{x}}_1^s,\mathrm{ }{\mathit{z}}_1^{s-1};{\mathrm{ }z}_s\right)\mathrm{ }(s=t,\mathrm{ }t+1,\mathrm{ }\cdots ,\mathrm{ }T)$ are already estimated from the previous subsection. Then we only need to model the probabilities $P({\mathit{x}}_k\mathrm{ }|\mathrm{ }{\mathit{x}}_1^{k-1},{\mathit{z}}_1^{k-1})$ by regression in the framework of single-point causal inference. As an illustration of our method, suppose that we obtain SNMM $\varphi \left({\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1};\mathrm{ }{z}_t\right)=f\left({\mathit{x}}_t,z_t;\mathrm{ }\mathit{\gamma }\right)$ and $P({\mathit{x}}_k\mathrm{ }|{\mathrm{ }\mathit{x}}_1^{k-1},{\mathit{z}}_1^{k-1})=P({\mathit{x}}_k\mathrm{ }|{\mathrm{ }\mathit{x}}_{k-1},z_{k-1}).$ By decomposing the SCE into the blip effects or applying (4) and (5), we obtain $$\text{SCE}\left({\mathit{D}}_1^T\right)=\sum _{s=1}^{T}E\left\{f\right({\mathit{x}}_t,z_t;\mathrm{ }\mathit{\gamma }\left)\right\},$$with respect to ${\prod }_{k=1}^{s}P({\mathit{x}}_k\mathrm{ }|\mathrm{ }{\mathit{x}}_{k-1},z_{k-1})$ and $$\text{SCE}\left({\mathit{x}}_t;\mathrm{ }{\mathit{D}}_t^T\right)=f\left({\mathit{x}}_t,z_t;\mathrm{ }\mathit{\gamma }\right)+\sum _{s=t+1}^{T}E\left\{f\right({\mathit{x}}_s,z_s;\mathrm{ }\mathit{\gamma }\left)\mathrm{ }\right|{\mathrm{ }\mathit{x}}_t,z_t\},$$with respect to ${\prod }_{k=t+1}^{s}P({\mathit{x}}_k\mathrm{ }|{\mathrm{ }\mathit{x}}_{k-1},z_{k-1}).\mathrm{ }$

With the estimates $\widehat{P}({\mathit{x}}_k\mathrm{ }|\mathrm{ }{\mathit{x}}_{k-1},z_{k-1})\mathrm{ }$ and $\widehat{\mathit{\gamma }},$ we obtain the estimate $\widehat{\mathrm{SCE}}\left({\mathit{x}}_t;\mathrm{ }{\mathit{D}}_t^T\right)$ under regime ${\mathit{D}}_t^T$ in stratum ${\mathit{x}}_t$ and $\widehat{\mathrm{SCE}}\left({\mathit{D}}_1^T\right)$ under regime ${\mathit{D}}_1^T$ in the population. The covariance matrix between the estimated SCEs under different regimes is estimated by the bootstrap method. If $\widehat{P}({\mathit{x}}_k\mathrm{ }|\mathrm{ }{\mathit{x}}_{k-1},z_{k-1}),$ in addition to $\widehat{\mathit{\gamma }},$ are consistent and asymptotically normal as typically required in single-point causal inference, so are $\widehat{\mathrm{SCE}}\left({\mathit{x}}_t;\mathrm{ }{\mathit{D}}_t^T\right)$ and $\widehat{\mathrm{SCE}}\left({\mathit{D}}_1^T\right).$ With the estimates and the covariance matrix, we conduct the Wald test on the SCE under any regime. We may also compare different regimes by testing the differences between their SCEs.

Several advantages of our method of modeling to estimate and test blip effects and SCEs are summarized as follows.

• When estimating the blip effect and SNMM, there is little risk for the model misspecification due to the following reasons. First, the regression model for the blip effect is the alternative G-formula (the decomposing formula) for blip effects, which is subject to no model misspecification. Second, the point effects as the response variables are estimated by regression in the framework of single-point causal inference. Third, the probability of treatments and covariates as the explanatory variables are estimated by the observed proportion in the data without modeling. Fourth, there is no restriction on the functional form of SNMM, and SNMM is testable in the framework of regression. SNMM improves the efficiency of estimating and testing blip effects. Notably, the modeling does not suffer from the curse of dimensionality and the null paradox as compared to the modeling via standard parameters (Hernan and Robins 2020; Robins 1997, 2009; Taubman et al. 2009).

• The covariate probabilities in the alternative G-formula (the decomposing formula) for SCEs are estimated by regression in the framework of single-point causal inference. With the estimated blip effects/SNMM and the estimated covariates probabilities, we use the alternative G-formula to calculate the estimate of SCE under any regime of treatments. We do not need to specify models for the baselines such as $E\left\{Y\left(0_t^{\mathit{T}}\right)|{\mathit{x}}_1^t,\mathrm{ }{\mathit{z}}_1^{t-1}\right\},$ which are subject to high risks for model misspecification (Hernan and Robins 2020; Robins 1997, 2004, 2009).

4. Estimating and testing SCEs of diagnosing and treating hospitals on cardia cancer survival

4.1. Setting

In the medical example of Sec. 2, the diagnosing hospital is the treatment variable $Z_1,$ which takes the value $z_1=0$ for small type and $z_1=1$ for large type. The treating hospital is the treatment variable $Z_2,$ which takes $z_2=0$ for small type and $z_2=1$ for large type. The following stationary covariates before $Z_1$ are measured: gender $X_{11},$ geographic area $X_{12}$ and age $X_{13}.$ Gender is $x_{11}=0$ for female and $x_{11}=1$ for male. Geographic area is categorized into $x_{12}=0$ for rural and $x_{12}=1$ for urban. Age takes continuous values $x_{13}.$ The time-dependent covariate between $Z_1$ and $Z_2$ is cancer stage $X_2$ taking the values $x_2=1,\mathrm{ }2,\mathrm{ }3,\mathrm{ }4.$ The outcome of interest is $Y,$ which takes $y=0$ for death and $y=1$ for survival within one year after diagnosis. In Supporting Material, the data is available together with the code for analyses in this article. The descriptive statistics are given in Table 1.

Table 1. Frequencies or means (standard deviations) of covariates and outcome across the diagnosing and treating hospitals for 157 cardia cancer patients.

	Cardia cancer: frequency or mean (standard deviation) on levels of $z_1$ and $z_2$
	$z_1=0$	$z_1=1$	$z_2=0$	$z_2=1$
$N$	69	88	39	118
$x_{11}=0$$=1$	14 55	26 62	10 29	30 88
$x_{12}=0$$=1$	43 26	39 49	27 12	55 63
$x_{13}$	68.04 (11.37)	65.38 (13.22)	71.90(9.20)	64.78(12.93)
$x_2=1$$=2$$=3$$=4$	7 7 21 34	12 13 17 46	3 3 9 24	16 17 29 56
One year survival $y=0$ $=1$	43 26	58 30	31 8	70 48

Hospital types:$z_1=0, 1$ for small and large diagnosing hospitals; $z_2=0, 1$ for small and large treating hospitals.

Stationary covariates: $x_{11}=0, 1$ for female and male; $x_{12}=0, 1$ for rural and urban residential areas; $x_{13}$ for age.

Time-dependent covariate: $x_2=1, 2, 3, 4$ for cancer stages.

Outcome: $y=0, 1$ death and survival within one year after diagnosis.

Due to the long-term social welfare system and relatively uniform culture in Sweden, most of the stationary covariates, such as education and socioeconomic status, have similar distributions across different hospitals and thus do not confound the causal effect. As a common practice in many epidemiologic researches in Sweden, we assume no unmeasured confounders for diagnosing hospitals $Z_1$ at least after conditional on gender $x_{11},$ residential area $x_{12},$ and age $x_{13}.$ Similarly, we assume no unmeasured confounders for treating hospitals $Z_2$ conditional on $z_1$ and cancer stage $x_2$ besides $x_{11},$ $x_{12},$ and $x_{13}.$ Together with the positivity and consistency assumptions, the assumption of no hidden confounders is the identifying condition (Robins 1986, 1997, 2009). Under the identifying condition, we identify the SCE under any regime of diagnosing and treating hospitals including the optimal regime by our data.

In the rest of this section, we will apply our method described in Sec. 3 to estimate and test the SCE under any regime of diagnosing and treating hospitals and evaluate the diagnosing hospitals.

4.2. Estimating and testing point effects of diagnosing and treating hospitals on one-year survival

In the area of cancer diagnosis, one usually studies one-year, three-year and five-year survivals. One-year survival is often the focus of highly malignant cancers such as cardia cancer.

To estimate the point effect of diagnosing hospital $z_1,$ we model the mean $\mu (x_{11},x_{12},x_{13},z_1),$ which is the probability of one-year survival $P(Y=1\mathrm{ }|\mathrm{ }{x}_{11},x_{12},x_{13},z_1)$ in stratum $\left(x_{11},x_{12},x_{13},z_1\right),$ where $z_1$ is the variable of interest while all others are possible confounders of $z_1.$ Notably, the model does not contain any posttreatment variable after $z_1,$ so the null paradox does not necessarily occur, and the model can be unsaturated without biasing the point effect. At the significance level of 0.10, we exclude the geographic area $x_{12}$ and age $x_{13}$ and obtain (6) $$\text{Logit}\left\{\mu \left(x_{11},z_1\right)\right\}=\log \left\{\frac{\mu \left(x_{11},z_1\right)}{1-\mu \left(x_{11},z_1\right)}\right\}={\beta }_0+{\beta }_1{x}_{11}+{\vartheta }_1{z}_1.\mathrm{ }$$(6)

From this model, we obtain the estimate for $$\mu \left(x_{11},z_1\right)=\frac{\hspace{0.17em}\text{exp}({\beta }_0+{\beta }_1{x}_{11}+{\vartheta }_1{z}_1)}{1+\hspace{0.17em}\text{exp}({\beta }_0+{\beta }_1{x}_{11}+{\vartheta }_1{z}_1)}$$

and then the point effect of large diagnosing hospital $z_1=1$ $$\vartheta \left(x_{11};\mathrm{ }{z}_1=1\right)=\mu \left(x_{11},\mathrm{ }{z}_1=1\right)-\mu \left(x_{11},z_1=0\right)\approx {\vartheta }_1.$$Here, the point effects for both $x_{11}=0$ and $x_{11}=1$ are approximately equal, i.e., $\vartheta \left(x_{11};\mathrm{ }{z}_1=1\right)\approx \mathrm{ }{\vartheta }_1$ with the p-value $=$ 0.6.

The estimate ${\widehat{\vartheta }}_1$ is presented in Table 2. Please note that the point effect is a mixed effect of both $z_1=1$ and $z_2=1$ and thus of little medical interest.

Table 2. Point effects, blip effects and optimal SCEs of diagnosing and treating hospitals on one-year survival of cardia cancer: estimate, p-value and 95 % CI.

	Estimate, p-value and 95 % CI in column for causal effects
			${\vartheta }_1$	${\vartheta }_2$
Point effect			$-$0.020 0.799 $(-0.171, 0.132)$	0.194 0.009 $(0.048, 0.340)$
Blip effect	${\gamma }_{10}$	${\gamma }_{11}$	${\gamma }_1$	${\gamma }_2$
	$-0.158$ 0.090 $(-0.341, 0.024)$	$-0.109$ 0.179 $(-0.267, 0.050)$	$-$0.121 0.142 $(-0.283, 0.040)$	0.194 0.009 $(0.048, 0.340)$
Optimal SCE	${\Gamma }_{10}$	${\Gamma }_{11}$	${\Gamma }_1$	${\Gamma }_2$
	0.194 0.009 $(0.048, 0.340)$	0.194 0.009 $(0.048, 0.340)$	0.194 0.009 $(0.048, 0.340)$	0.194 0.009 $(0.048, 0.340)$

Two point effects: ${\vartheta }_1$ of large diagnosing hospital; ${\vartheta }_2$ of large treating hospital.

Four blip effects: ${\gamma }_{10}$ of large diagnosing hospital among women; ${\gamma }_{11}$ of large diagnosing hospital among men; ${\gamma }_1$ of large diagnosing hospital in the population; ${\gamma }_2$ of large treating hospital.

Four optimal SCEs: ${\Gamma }_{10}$ under optimal regime of diagnosing and treating hospital among women; ${\Gamma }_{11}$ under optimal regime of diagnosing and treating hospital among men; ${\Gamma }_1$ under optimal regime of diagnosing and treating hospital in the population; ${\Gamma }_2$ under optimal regime of treating hospital in population. The estimated optimal treating hospital is ${\widehat{O}}_2=1$ (the large treating hospital). The estimated optimal diagnosing hospital is ${\widehat{O}}_1=0$ (small diagnosing hospital).

To estimate the point effect of treating hospital $z_2,$ we model the mean $\mu (x_{11},x_{12},x_{13},z_1,x_2,z_2),$ which is the probability of one-year survival $P(Y=1\mathrm{ }|\mathrm{ }{x}_{11},x_{12},x_{13},z_1,x_2,z_2)$ in stratum $\left(x_{11},x_{12},x_{13},z_1,x_2,z_2\right),$ where $z_2$ is the variable of interest while all others are possible confounders of $z_2.$ At the significance level of 0.10, we exclude $x_{11},x_{12}$ and $x_{13}$ and obtain (7) $$\text{Logit}\mathrm{ }\left\{\mu \left(z_1,x_2,z_2\right)\right\}=\log \left\{\frac{\mu \left(z_1,x_2,z_2\right)}{1-\mu {(z}_1,x_2,z_2)}\right\}={\alpha }_0+{\alpha }_1{z}_1+{\alpha }_2{x}_2+{\vartheta }_2{z}_2,$$(7)

From this model, we obtain the estimate for the point effect of large treating hospital $z_2=1:$ $$\vartheta \left(z_1,x_2;z_2=1\right)=\mu \left(z_1,x_2,z_2=1\right)-\mu \left(z_1,x_2,z_2=0\right)\approx {\vartheta }_2$$Here, the point effects in all strata $\left(z_1,x_2\right)$ are approximately equal, i.e., ${\vartheta \left(z_1,x_2;z_2=1\right)\approx \vartheta }_2$ with the p-value $0.2.$ The estimate ${\widehat{\vartheta }}_2$ is presented in Table 2 together with ${\widehat{\vartheta }}_1.$

As it will be clear in the next subsection, we will need the conditional variances of the estimated point effects given all the covariates, the diagnosing and treating hospitals. The conditional variances are obtained by specifying a constant dispersion parameter for the binomial distribution and applying the obtained models (6) and (7). Noticeably, the conditional covariance between the estimated point effects of diagnosing and treating hospitals is approximately zero (Wang and Yin 2015, 2020).

4.3. Estimating and testing blip effect of diagnosing and treating hospitals on one-year survival

Because neither (6) nor (7) contains geographic area $x_{12}$ and age $x_{13},$ we do not include them in the blip effects and SCEs below. Because the blip effects are a special type of SCEs and determine all other SCEs (see the next subsection), we will estimate the blip effects here and SCEs in the next subsection. The blip effect $\varphi \left(x_{11};\mathrm{ }{z}_1\right)$ of diagnosing hospital $z_1$ in the stratum ${(x}_{11})$ is an increase of the mean of potential survival $Y\left(D_1,D_2\right)$ under regimes ${(D}_1,D_2)=(z_1,0)$ versus to ${(D}_1,D_2)=(0,\mathrm{ }0)$ in the stratum, namely $$\varphi \left(x_{11};\mathrm{ }{z}_1\right)=E\left\{Y\left(z_1,\mathrm{ }0\right)\mathrm{ }|\mathrm{ }{x}_{11}\right\}-E\left\{Y\left(0,\mathrm{ }0\right)\mathrm{ }|\mathrm{ }{x}_{11}\right\}.$$Clearly, we have $\varphi \left(x_{11};\mathrm{ }{z}_1=0\right)=0.$ The blip effect $\varphi \left(x_{11},z_1,x_2;\mathrm{ }{z}_2\right)$ of treating hospital $z_2$ in the stratum ${(x}_{11},z_1,x_2)$ is an increase of the mean of potential survival $Y\left(D_2\right)$ under regimes ${D_2=z}_2$ versus $D_2=0\mathrm{ }$ in the stratum, namely, $$\varphi \left(x_{11},z_1,x_2;\mathrm{ }{z}_2\right)=E\left\{Y\left(z_2\right)\mathrm{ }|{\mathrm{ }x}_{11},z_1,x_2\right\}-E\left\{Y\left(0\right)\mathrm{ }|\mathrm{ }{x}_{11},z_1,x_2\right\}.$$clearly, $\varphi \left(x_{11},z_1,x_2;\mathrm{ }{z}_2=0\right)=0.$ Because the treating hospital $Z_2$ is the last treatment variable, the blip effect $\varphi \left(x_{11},z_1,x_2;\mathrm{ }{z}_2=1\right)\mathrm{ }$ is equal to the point effect $\vartheta \left(z_1,x_2;\mathrm{ }{z}_2=1\right)={\vartheta }_2$ obtained in the previous subsection. Furthermore, a question of medical interest is whether the blip effect of $z_1$ depends on gender $x_{11}.$ Motivating by all these together, we suppose that the blip effects satisfy SNMM of the form (8) $$\begin{array}{c}\varphi \left(x_{11}=0;{\mathrm{ }z}_1\right)=z_1{\gamma }_{10},\\ \varphi \left(x_{11}=1;{\mathrm{ }z}_1\right)=z_1{\gamma }_{11},\\ \mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2\right)={\mathrm{z}}_2{\gamma }_2.\end{array}$$(8)

Let $\mathit{\gamma }=({\gamma }_{10},{\gamma }_{11},{\gamma }_2)$ be the vector of all parameters in (9). The parameter ${\gamma }_{10}$ is the blip effect of large diagnosing hospital among women ($x_{11}=0)\mathrm{ }$ and ${\gamma }_{11}$ is blip effect of large diagnosing hospital among men ($x_{11}=1)\mathrm{ }.$ They describe modification of the blip effect of diagnosing hospital by gender. The average of the blip effect $\varphi \left(x_{11};z_1=1\right)$ with respect to $P(x_{11})$ is blip effect of large diagnosing hospital in the population and equal to ${\gamma }_1={\gamma }_{10}P{(x}_{11}=0)\mathrm{ }+{\gamma }_{11}P{(x}_{11}=1).$

Intuitively, the point effect $\vartheta \left(x_{11};{\mathrm{ }z}_1=1\right)$ results from the blip effects of $z_1=1$ and $z_2=1.$ For $z_1=1,$ the contribution is ${\gamma }_{10}$ if $x_{11}=0$ and ${\gamma }_{11}$ if $x_{11}=1.$ For $z_2=1,$ the contribution comes from the two strata ${(x}_{11},z_1=1,x_2)$ and ${(x}_{11},\mathrm{ }{z}_1=0,x_2),$ and it is equal to ${\gamma }_{2}C\left(x_{11}\right),$ where $$C\left(x_{11}\right)={\sum }_{x_2}\left[P\right\{{(\mathrm{x}}_2,\mathrm{ }{\mathrm{z}}_2=1)\mathrm{ }\left|\mathrm{ }{x_{11},\mathrm{ }z}_1=1\right\}-P\left\{\right({\mathrm{x}}_2,{\mathrm{z}}_2=1\left)\mathrm{ }\right|{{\mathrm{ }x}_{11},\mathrm{ }z}_1=0\}.$$Thus, we have

$\vartheta \left(x_{11}=0;z_1=1\right)={\gamma }_{10}+{\mathrm{\gamma }}_{2}C(x_{11}=0)\mathrm{ }$ and $\vartheta \left(x_{11}=1;z_1=1\right)={\gamma }_{11}+{\mathrm{\gamma }}_{2}C\left(x_{11}=1\right).$ The $C\left(x_{11}\right)$ is estimated by the observed proportion in the data without modeling. For the last treatment variable $Z_2,$ the point effect of $z_2=1$ is equal to the blip effect of $z_2=1,$ namely, $${\gamma }_2=\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;{\mathrm{z}}_2=1\right)=\vartheta \left(z_1,x_2;z_2=1\right)={\vartheta }_2.$$Taking these together, we have the following decomposing formula under SNMM (8) (9) $$\begin{array}{c}\vartheta \left(x_{11}=0;\mathrm{ }{z}_1=1\right)={\gamma }_{10}+{\mathrm{\gamma }}_{2}C\left(x_{11}=0\right),\mathrm{ }\\ \vartheta \left(x_{11}=1;{\mathrm{ }z}_1=1\right)={\gamma }_{11}+{\mathrm{\gamma }}_{2}C\left(x_{11}=1\right),\mathrm{ }\\ {\mathrm{\vartheta }}_2={\gamma }_2.\end{array}$$(9)

The formula (9) is the alternative G-formula for the blip effect under SNMM (8) and it can also be derived by applying the formula (3) and SNMM (8) as follows. By applying (3) to $\mathrm{t}=1,\mathrm{ }\mathrm{T}=2\mathrm{ }$ we obtain $$\begin{array}{c}\vartheta \left(x_{11};{ z}_1=1\right)=\\ \varphi \left(x_{11}; z_1=1\right)+E\left\{\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2\right)\mathrm{ }\right|{ x}_{11}, z_1=1\}-E\{\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2\right) |x_{11}, z_1=0\}\end{array}$$Applying (8), we obtain for$\mathrm{ }{x}_{11}=0$ $$\begin{array}{c}\begin{array}{c}\vartheta \left(x_{11}=0; z_1=1\right)\\ ={\gamma }_{10}+{\gamma }_2 {\sum }_{x_2}\left[P\right\{{(\mathrm{x}}_2,\mathrm{ }{\mathrm{z}}_2=1)\mathrm{ }\left|{{ x}_{11}=0, z}_1=1\right\}-P\left\{\right({\mathrm{x}}_2,{\mathrm{z}}_2=1\left)\right|{{ x}_{11}=0, z}_1=0\}]\end{array}\\ ={\gamma }_{10}+{\gamma }_{2}C(x_{11}=0)\end{array}$$Thus, we have $\vartheta \left(x_{11}=0;\mathrm{ }{z}_1=1\right)={\gamma }_{10}+{\mathrm{\gamma }}_{2}C\left(x_{11}=0\right).$ Similarly we have $\vartheta \left(x_{11}=1;z_1=1\right)={\gamma }_{11}+{\mathrm{\gamma }}_{2}C(x_{11}=1).$ By applying (3) to $\mathrm{t}=\mathrm{ }\mathrm{T}=2\mathrm{ }$ we obtain ${\vartheta }_2=\vartheta \left(z_1,x_2;\mathrm{ }{z}_2=1\right)=\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2=1\right)={\gamma }_2.$ That is ${\vartheta }_2={\gamma }_2.$

Because $\vartheta \left(x_{11};\mathrm{ }{z}_1=1\right)={\vartheta }_1$ for both $x_{11}=0\text{ and }1$ as seen from model (6), we have $\widehat{\vartheta }\left(x_{11};\mathrm{ }{z}_1)=1\right)={\widehat{\vartheta }}_1,$ but the variance $\text{var}\left\{\widehat{\vartheta }\right(x_{11};{\mathrm{ }z}_1)=1\left)\right\}$ is obtained by adjusting $\text{var}({\widehat{\vartheta }}_1)$ to the size of stratum $x_{11}.$ All probabilities in $C(x_{11})$ are estimated by the corresponding proportions in the data without modeling. The ${\widehat{\vartheta }}_2$ is estimated from model (7).

Now, conditional on all covariates and diagnosing and treating hospitals, we use (9) as a regression model to estimate $\mathit{\gamma }={(\gamma }_{10},$ ${\gamma }_{11},$ ${\gamma }_2$), where the response variables are $\widehat{\vartheta }\left(x_{11}=0;\mathrm{ }{z}_1=1\right),$ $\widehat{\vartheta }\left(x_{11}=1;\mathrm{ }{z}_1=1\right),$ and ${\widehat{\vartheta }}_2$ while the explanatory variables are the estimates $\widehat{C}(x_{11}=0),$ $\widehat{C}\left(x_{11}=1\right)$ and the ones. The bootstrap method is used to obtain the covariance matrix $\text{cov}\mathrm{(}\widehat{\mathit{\gamma }}\mathrm{)}$ incorporating the variability of all covariates and diagnosing and treating hospitals. With $\widehat{\mathit{\gamma }}$ and $\text{cov}\mathrm{(}\widehat{\mathit{\gamma }}\mathrm{)},$ we conduct the Wald test on $\mathit{\gamma }.$ The result is presented in Table 2.

4.4. Estimate and testing SCE under any regime of diagnosing and treating hospitals on one-year survival

First, we consider $\text{SCE}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right).$ It is an increase of the mean of potential survival under regimes $D_2$ versus $0$ in stratum ${(x}_{11},z_1,x_2),$ namely $$\text{SCE}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right)=E\left\{Y\left(D_2\right)\mathrm{ }|\mathrm{ }{x}_{11},z_1,x_2\right\}-E\left\{Y\left(0\right)\mathrm{ }|{\mathrm{ }x}_{11},z_1,x_2\right\}$$where $D_2$ potentially and deterministically assigned treating hospital $z_2$ to each patient. Noticeably, $D_2$ can be dynamic, for instance, $D_2=1$ for patients younger than 55 years and $D_2=0$ otherwise. By applying (5) to $\mathrm{t}=\mathrm{ }\mathrm{T}=2,\mathrm{ }$ we see that $\text{SCE}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right)=\mathrm{ }\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2\right).$ Thus, by applying $\varphi \left(x_{11},z_1,x_2;\mathrm{ }{z}_2\right)$ in SNMM (8), we obtain following decomposing formula for the SCE (10) $$\text{SCE}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right)={\gamma }_2{z}_2$$(10)

Second, we consider $\text{SCE}\left(x_{11};{\mathrm{ }D}_1,D_2\right),$ which is an increase of the mean of potential survival under regime ${(D}_1,D_2)$ versus $(0,\mathrm{ }0)$ in stratum $(x_{11}),$ namely $$\mathrm{SCE}\left(x_{11};{\mathrm{ }D}_1,D_2\right)=E\left\{Y\left({D_1,\mathrm{ }D}_2\right)\mathrm{ }|\mathrm{ }{x}_{11}\right\}-E\left\{Y\left(0,\mathrm{ }0\right)\mathrm{ }|\mathrm{ }{x}_{11}\right\}$$where ${(D}_1,D_2)$ potentially and deterministically assigned ${(z}_1,\mathrm{ }{z}_2)$ of diagnosing and treating hospitals and it can be dynamic. If regime ${(D}_1,D_2)=(z_1,\mathrm{ }0),$ then $\mathrm{SCE}\left(x_{11};\mathrm{ }{D}_1,D_2\right)=\varphi \left(x_{11};\mathrm{ }{z}_1\right),$ which implies that the blip effect is a special type of SCEs. Intuitively the SCE should be a result of the blip effects of $z_1=1$ and $z_2=1$ in the regime. The contribution of $z_1$ is ${\gamma }_{10}{z}_1$ if $x_{11}=0$ and ${\gamma }_{11}{z}_1$ if $x_{11}=1.$ The contribution of $z_2$ is ${\gamma }_2{z}_2.$ Therefore, we have the following decomposing formula (11) $$\begin{array}{c}\text{SCE}\left(x_{11}=0;{\mathrm{ }D}_1,D_2\right)={\gamma }_{10}{z}_1+{\mathrm{\gamma }}_2{z}_2,\mathrm{ }\\ \mathrm{ }\text{SCE}\left(x_{11}=1;\mathrm{ }{D}_1,D_2\right)={\gamma }_{11}{z}_1+{\mathrm{\gamma }}_2{z}_2.\mathrm{ }\end{array}$$(11)

The formulas (10) and (11) implies that all SCEs are determined by blip effects. They are the alternative G-formula for the SCE under SNMM (8) and can also be derived by applying the formula (5) and the SNMM (8) as follows. By applying (5) to $\mathrm{t}=1,\mathrm{ }\mathrm{T}=2,\mathrm{ }$ we obtain $$\text{SCE}\left(x_{11};{\mathrm{ }D}_1,D_2\right)=\varphi \left(x_{11};\mathrm{ }{z}_1\right)+\mathrm{ }E\left\{\mathrm{\varphi }\left({\mathrm{x}}_{11},{\mathrm{z}}_1,{\mathrm{x}}_2;\mathrm{ }{\mathrm{z}}_2\right)x_{11},z_1\right\}$$Now by applying (8), we obtain(11).

Taking the average of (11) with respect to ${P(x}_{11}),$ we obtain the alternative G-formula for $\text{SCE}\left(D_1,D_2\right),$ which is an increase of the mean of the potential survival under regimes $(D_1,D_2)$ versus $(0,\mathrm{ }0)$ in the population.

In (10) replacing $\mathit{\gamma }={(\gamma }_{10},$ ${\gamma }_{11},$ ${\gamma }_2$) by estimates$\mathrm{ }\widehat{\mathit{\gamma }}$ obtained from (9), we obtain the estimate for $\text{SCE}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right)$ under any treating hospital $D_2.$ In (11) replacing $\mathit{\gamma }={(\gamma }_{10},$ ${\gamma }_{11},$ ${\gamma }_2$) by estimates$\mathrm{ }\widehat{\mathit{\gamma }}$ obtained from (9), we obtain the estimate for $\text{SCE}\left(x_{11};{\mathrm{ }D}_1,D_2\right)$ under any regime ${(D}_1,D_2).$ Averaging $\widehat{\mathit{SCE}}\left(x_{11};\mathrm{ }{D}_1,D_2\right)$ with respect to proportion ${\widehat{P}(x}_{11}),$ we obtain the estimate for $\text{SCE}\left(D_1,D_2\right).$ The variance of the estimated SCE is obtained by the bootstrap method. With the estimate and its variance, we obtain the confidence interval of the SCE. The result is presented in Table 2.

4.5. Estimating and testing optimal regime of diagnosing and treating hospitals on one-year survival

By the dynamic programming procedure (Bellman, 1957), we estimate and test optimal regime ${(O}_1,O_2)$ of diagnosing and treating hospitals by estimating and testing SCEs, as follows. We first estimate the optimal treating hospital $O_2$ by $${\widehat{O}}_2=\text{arg}\mathrm{ }\text{ma}{x}_{D_2\in \{0,1\}}\widehat{\mathrm{SCE}}\left(x_{11},z_1,x_2;\mathrm{ }{D}_2\right),$$which is the treating hospital such that $\widehat{\mathrm{SCE}}\left(x_{11},z_1,x_2;{\mathrm{ }D}_2\right)={\widehat{\gamma }}_2{z}_2$ achieves the maximum value among all possible values (0 or 1) that $D_2$ assigns for given $(x_{11},z_1,x_2).$ Because ${\widehat{\gamma }}_2=0.194>0$ (see Table 2), we have ${\widehat{O}}_2=1$ (the large treating hospital) and $\widehat{\mathrm{SCE}}\left(x_{11},z_1,x_2;{\widehat{O}}_2\right)={\widehat{\gamma }}_2$ for all given $(x_{11},z_1,x_2).$ We test ${\widehat{O}}_2$ by testing ${\gamma }_2$ and obtain the 95% CI and the p-value. Now we estimate the optimal diagnosing hospital $O_1$ by $${\widehat{O}}_1=\text{arg}\mathrm{ }\text{ma}{x}_{D_1\in \{0,1\}}\widehat{\mathrm{SCE}}\left(x_{11};{\mathrm{ }D}_1,{\widehat{O}}_2\right).$$

Because ${\widehat{\gamma }}_{10}=-0.158$ and ${\widehat{\gamma }}_{11}=-0.109$ we have $$\widehat{\mathrm{SCE}}\left(x_{11}=0;{\mathrm{ }D}_1,{\widehat{O}}_2\right)={\widehat{\gamma }}_{10}{z}_1+{\widehat{\gamma }}_2=-0.158z_1+0.194,$$ $$\widehat{\mathrm{SCE}}\left(x_{11}=1;{\mathrm{ }D}_1,{\widehat{O}}_2\right)={\widehat{\gamma }}_{11}{z}_1+{\widehat{\gamma }}_2=-0.109z_1+0.194.$$

We see that $\widehat{\mathrm{SCE}}\left(x_{11};\mathrm{ }{D}_1,{\widehat{O}}_2\right)$ take maximum value when $z_1=0$ for both $x_{11}=0$ and $x_{11}=1.$ Therefore ${\widehat{O}}_1=0\mathrm{ }$(small diagnosing hospital) and $\widehat{\mathrm{SCE}}\left(x_{11};\mathrm{ }{\widehat{O}}_1,{\widehat{O}}_2\right)=0.194\mathrm{ }$ for any given gender $x_{11}.$ Thus, the estimates for $${\mathrm{\Gamma }}_{10}=\mathrm{SCE}\left(x_{11}=0;{\mathrm{ }\widehat{O}}_1,{\widehat{O}}_2\right), {\mathrm{\Gamma }}_{11}=\mathrm{SCE}\left(x_{11}=1;{\mathrm{ }\widehat{O}}_1,{\widehat{O}}_2\right),\mathrm{\hspace{0.5em}and\hspace{0.5em}}{\mathrm{\Gamma }}_1=\mathrm{SCE}\left({\widehat{O}}_1,{\widehat{O}}_2\right)$$are all equal to $0.194.$ We test ${(\widehat{O}}_1,{\widehat{O}}_2)$ by testing $\mathrm{SCE}\left(x_{11};\mathrm{ }{\widehat{O}}_1,{\widehat{O}}_2\right)$ and obtain the 95% CI and the p-value for $\mathrm{SCE}\left(x_{11};\mathrm{ }{\widehat{O}}_1,{\widehat{O}}_2\right).$ The result is presented in Table 2.

4.6. Causal analysis of diagnosing and treating hospitals based on Table 2

For the treating hospital, the point effect ${\vartheta }_2$ or the blip effect ${\gamma }_2$ of large treating hospital is estimated at 0.194 with the p-value equal to 0.009, implying that the large treating hospital is superior to one-year survival. There is no effect modification by the cancer stage, gender and age. This observation is consistent with the medical knowledge. Cardia cancer is highly malicious regardless of age and gender. It progressed quickly regardless of cancer stage. Patients typically have comorbidities and large treating hospitals are probably better in dealing with comorbidities.

For the diagnosing hospital, the point effect ${\vartheta }_1$ is estimated at $-$0.020 with the p-value 0.799, suggesting large and small diagnosing hospitals perform equally well. But it is misleading, because the point effect is a mixed effect of diagnosing and treating hospitals. The blip effect ${\gamma }_1$ is estimated at $-$0.121 with the p-value 0.142, implying that the small diagnosing hospital performs better for one-year survival. This observation indicates that the early diagnosis of such a malicious cancer as cardia cancer depends on the short waiting queue and attention typically at small diagnosing hospitals. The blip effect ${\gamma }_{10}$ for women is estimated at $-0.158$ with the p-value equal to 0.090, and the blip effect ${\gamma }_{11}$ for men at $-0.109$ with the p-value equal to 0.179, suggesting some modification of the blip effect of diagnosing hospital by gender.

The optimal regime is estimated as the small diagnosing hospital and large treating hospital. The optimal SCE is estimated at 0.194 with the p-value 0.009. This reveals the potential improvement for one-year survival by the optimal regime.

Here, all effects are measured by the difference in mean survival, but they can also be measured by the difference in a function of mean survival. We may apply the alternative G-formula (namely, the decomposing formula) to estimate the causal effect measured by odds ratio, rate ratio and hazard ratio. Though the estimates are consistent, the problem of non-collapsibility of a non-linear measure becomes far worse in the context of sequential treatments, leading to a biased estimate for finite sample. Therefore, it is recommended that one uses the linear measure for the SCE.

5. Comparison of our method with available methods

Method (i) is our method described in Sec. 4. Method (ii) is based on the well-known G-formula (Hernan and Robins 2020; Robins 1997, 2009; Taubman et al. 2009). Method (iii) is the marginal structural model based on inverse probability of treatment weighting (Hernan and Robins 2020; Hernan, Brumback, and Robins 2000; Robins, Brumback, and Hernan 2000; Robins 2009). Method (iv) is the G-estimation based on SNMM or optimal SNMM (Hernan and Robins 2020; Robins 1997, 2004, 2009). Method (i) and (ii) are parametric whereas methods (iii) and (iv) are semi-parametric. Methods (ii)-(iv) are also described in the context of this article in Supplemental Material.

Because methods (i)-(iv) may lead to the same inference of causal effect of treating hospital, we only present the causal effect of diagnosing hospital. In Sec. 4, we have applied method (i) to estimate the blip effect ${\gamma }_{10}$ of large diagnosing hospital among women; ${\gamma }_{11}$ of large diagnosing hospital among men; ${\gamma }_1\mathrm{ }$ of large diagnosing hospital in the population; optimal SCE $({\mathrm{\Gamma }}_1)$ under optimal regime of diagnosing and treating hospitals versus small ones in the population. Here, we additionally estimate the optimal blip effect of large diagnosing hospital: ${\mathrm{{\rm Y}}}_1=\text{SCE}\left(D_1=1,{\widehat{O}}_2\right)-\text{SCE}\left(D_1=0,{\widehat{O}}_2\right),$ because method (iv) does not estimate ${\mathrm{\Gamma }}_1.$ The result from these methods is presented in Table 3. From this table, we have the following observations.

Table 3. Comparison of our method with available methods in Sec. 5: estimate, p-value and 95 % confidence interval (95% CI) for causal effects of diagnosing and treating hospitals on one-year survival of cardia cancer.

	Estimate, p-value and 95 % CI in column of the causal effect
	${\gamma }_{10}$	${\gamma }_{11}$	${\gamma }_1$	${\mathrm{{\rm Y}}}_1$	${\Gamma }_1$
Method (i)	$-0.158$ 0.090 $(-0.341, 0.024)$	$-0.109$ 0.179 $(-0.267, 0.050)$	$-$0.121 0.142 $(-0.283, 0.040)$	$-$0.121 0.142 $(-0.283, 0.040)$	0.194 0.009 $(0.048, 0.340)$
Method (ii)			$-$0.066 0.605 $(-0.316, 0.184)$	$-$0.101 0.392 $(-0.333, 0.131)$	0.200 0.048 $(0.002, 0.399)$
Method (iii)			0.463 0.259 $(-0.342, 1.267)$	$-$0.167 0.087 $(-0.359, 0.024)$	0.249 0.012 $(-0.056, 0.443)$
Method (iv)	$-0.242$ 0.078 $(-0.509, 0.027)$	$-0.092$ 0.369 $(-0.291, 0.108)$	$-$0.127 0.162 $(-0.305, 0.051)$	$-$0.127 0.162 $(-0.305, 0.051)$

Four estimation methods: method (i) our method; method (ii) the parametric method based on the well-known G-formula; (iii) the marginal structural model based on inverse probability of treatment weighting; method (iv) the G-estimation based on SNMM or optimal SNMM; Empty cells imply that they are not easily estimable by the method.

Causal effects: ${\gamma }_{10}$ blip effect of large diagnosing hospital among women, ${\gamma }_{11}$ blip effect of large diagnosing hospital among men; ${\gamma }_1$ the blip effect of large diagnosing hospital in the population; ${\mathrm{{\rm Y}}}_1$ the optimal blip effect of large diagnosing hospital in the population; ${\Gamma }_1$ the SCE under the optimal regime of diagnosing and treating hospitals in the population. The estimated optimal treating hospital is ${\widehat{O}}_2=1$ (the large treating hospital). The estimated optimal diagnosing hospital is ${\widehat{O}}_1=0$ (small diagnosing hospital).

Method (i) yields the estimates for all five causal effects: $${\widehat{\gamma }}_{10}=-0.158,{\widehat{\gamma }}_{11}=-0.109,{\widehat{\gamma }}_1=-0.121,{\widehat{\mathrm{{\rm Y}}}}_1=-0.121,\mathrm{ }\text{and\hspace{0.5em}}{\widehat{\mathrm{\Gamma }}}_1=0.194.$$Here we only need to estimate a total of two parameters: ${\vartheta }_1$ for the point effect of $z_1$ and$\mathrm{ }{\vartheta }_2$ for the point effect of $z_2,$ to estimate the blip effect and the SCE under any regime, as described in Sec. 4. The estimates are in consistency with the medical knowledge. Because SCEs are all estimated from the same estimated SNMM, we can compare different regimes by testing their SCEs, for instance, we compare ${(\widehat{O}}_1,{\widehat{O}}_2)\mathrm{ }$ with $(0,\mathrm{ }0)$ by testing ${\mathrm{\Gamma }}_1.$

Method (ii) only yields the estimates for three causal effects: $${\widehat{\gamma }}_1=-0.066,{\widehat{\mathrm{{\rm Y}}}}_1=-0.101,\mathrm{ }\text{and\hspace{0.5em}}{\widehat{\Gamma }}_1=-0.200,$$because of imbalanced and sparse data between exposure sequence ${(z}_1,z_2)$ and gender $x_{11}.$ Here one needs to estimate a total of $2\times 2\times 4\times 2=32$ standard parameters$\mathrm{ }\mu \left(x_{11},{z_1,x}_2,z_2\right)$ to estimates these causal effects (curse of dimensionality), which leads to wide confidence intervals. The estimate ${\widehat{\Gamma }}_1=-0.200$ is biased, because the optimal SCE should be larger than or equal to zero. The bias may be due to the null paradox (Hernan and Robins 2020; Robins 1997, 2009; Taubman et al. 2009). For a long sequence of cancer diagnoses and treatments, the problems with the null paradox and the curse of dimensionality in method (ii) become far worse. In this case, method (i) may provide a useful tool for estimating and testing the causal effect of cancer diagnosis.

Method (iii) also yields the estimates for only three causal effects: $${\widehat{\gamma }}_1=0.463,{\widehat{\mathrm{{\rm Y}}}}_1=-0.167,\mathrm{ }\text{and\hspace{0.5em}}{\widehat{\Gamma }}_1=0.249.$$

The estimate ${\widehat{\gamma }}_1=0.463$ is not medically sensible. Also notice the wide confidence intervals. The reason is due to imbalanced and sparse data between exposure sequence ${(z}_1,z_2)$ and gender $x_{11}$ (Hernan and Robins 2020; Hernan, Brumback, and Robins 2000; Robins, Brumback, and Hernan 2000; Robins 2009).

Method (iv) yields estimates for four causal estimates: $${\widehat{\gamma }}_{10}=-0.241,{\widehat{\gamma }}_{11}=-0.092,{\widehat{\gamma }}_1=-0.127,\mathrm{ }\text{and\hspace{0.5em}}{\widehat{\mathrm{{\rm Y}}}}_1=-0.127.$$The estimates ${\widehat{\gamma }}_1=-0.127$ and ${\widehat{\mathrm{{\rm Y}}}}_1=-0.127\mathrm{ }$ are rather close to those obtained from method (i), albeit with lower significance. However, it does not estimate ${\mathrm{\Gamma }}_1.$ Generally, method (iv) does not facilitate comparison between different regimes.

6. Conclusion

In many practices, a single-point treatment can hardly achieve the desirable result. More often, a sequence of treatments are assigned or implemented, where a new situation often arises from the early treatments and also influences the assignment of subsequent treatments. A typical example is the stochastic process of cancer diagnosis. In such circumstances, it is highly difficult to design and conduct sequential randomized trials in comparison to randomized trials for single-point treatments. Consequently, data arising from a sequence of treatments are often observational, and this requires statistical modeling when estimating and testing SCEs under specified regimes of treatments.

In this article, we introduce a parametric method of estimating and testing SCEs in the observational study settings by modeling based on the alternative G-formula. Our method is implemented in three steps: first to estimate and test the point effect, second to use the estimated point effect to estimate and test the blip effect, and finally to use the estimated blip effects to estimate and test SCEs. Each of these steps is intuitive and can be carried out using regressions familiar to epidemiologists in single-point causal inference.

With its simple setting, we can readily extend the example in this article to other observational studies. As to cancer diagnosis, we can similarly study various cancer types, different cancer outcomes such as cancer progression, quality of life and others, different diagnosing techniques such as the biomarkers or a sequence of screening steps, and different modification factors such as social-economic status, comorbidity, and family history.

Ethics approval

The proposed research is covered by the ethical committee approval (DNR880113/13, §121) from the ethical review board of Uppsala University.

This work was supported by Vetenskapsrådet, Sweden (2019-02913).

Supplemental material

1. Data and code available in [Zenodo], at https://doi.org/10.5281/zenodo.6367502

2. A description of available methods (ii), (iii) and (iv) in the context of the medical example.

Additional information

Funding

This work is partially supported by the Swedish Research Council with the grant number 2019-02913.