Let D_t (t = 1, …, T) be a treatment plan, which would deterministically assign treatments z_t to each unit of the population, although possibly contrary to fact. A treatment regime is a sequence of such treatment plans, D1T={D1,…,DT}. Prior to D₁, there exists a stationary covariate vector X₁. Under D1T, each unit could have a potential time-dependent covariate vector Xt(D1t-1) between D_t−1 and D_t (t = 2, …, T) and a potential outcome Y(D1T) of interest, which is assumed to occur after the last treatment plan D_T without loss of generality. The stochastic process is illustrated by

These potential variables are denoted by {X1T(D1T-1),D1T,Y(D1T)}.. Their realizations are {x1T,z1T,y}={x1,z1,x2,…,xT,zT,y}. Let Ht(D1t-1)={X1t(D1t-1),D1t-1} be the history of potential covariates and treatment plans before treatment plan D_t. Its realization is ht={x1t,z1t-1}. Given h_t, we document the (sub) regime in stratum h_t as DtT={Dt,…,DT}, the potential covariate vector between D_s−1 and D_s as Xs(Dts-1) (s = t + 1, …, T), and the potential outcome as Y(DtT).

Without loss of generality, we take z_t = 0 as the control treatment and DtT=0={0,0,…,0} the control regime. Given h_t, consider two regimes DtT={zt,0} and {0, 0} in stratum h_t. We have the potential outcomes Y(z_t, 0) and Y(0, 0). According to Robins [1] and Hernan and Robins [5], the blip effect of treatment z_t in stratum h_t is defined as

where the (conditional) expectations are with respect to the conditional distributions P{Y(z_t, 0)∣h_t} and P{Y(0, 0)∣h_t}. Clearly, ϕ(h_t; 0) = 0. The blip effect describes the net effect of treatment in sequence [2] and its effect modification by the previous covariates and treatments [3, 4]. Furthermore, the blip effects of treatments are the determining factors for the causal effect of any regime of treatments according to Robins [1] and Hernan and Robins [5].

Corresponding to the potential variables, we have a sequence of observable treatment variables Z1T={Z1,Z2,…,ZT}, a sequence of observable covariate vectors X1T={X1,X2,…,XT}, and the observable outcome Y. The stochastic process is illustrated by

Suppose that the observable variables {X1T,Z1T,Y} have the same support as the potential variables {X1T(D1T-1),D1T,Y(D1T)}, that is, the observable variables take the same values as the potential variables. Like the realizations of the potential variables, the observed values of the observable variables are denoted by {z1T,x1T,y}. Let Ht={X1t,Z1t-1} be the history of covariates and treatment variables before treatment variable Z_t. Its realization is ht={x1t,z1t-1}. In the following, we will use P(.) to denote the probability distribution of discrete variables or the density distribution of continuous variables. The joint distribution of {X1T,Z1T,Y} factorizes into

To identify the blip effect by the observable variables {X1T,Z1T,Y}, Robins introduced the identifying condition [1, 6]: (a) The consistency assumption: if the observed treatments z1T are equal to the realizations of D1T, then the observed covariates x1T are equal to the realizations of the potential covariates X1T(D1T-1), and the observed outcome y is equal to the realization of the potential outcome Y(D1T). (b) The assumption of no unmeasured confounders: given h_t, treatment Z_t is conditionally independent of Xt+1T(DtT-1) and Y(DtT). (c) The positivity assumption: if P(H_t = h_t)>0, then 0 < P(Z_t∣h_t) < 1.

The standard parameter for the conditional distribution P(Y∣h_T, z_T) is the conditional mean μ(h_T, z_T) = E(Y∣h_T, z_T). Under the identifying condition, Robins expressed the blip effect in terms of the standard parameters [1, 6]. If the time-dependent covariate X_t is a posttreatment variable from the earlier treatments Z1t-1 and a confounder for the subsequent treatments ZtT, however, it is highly difficult to specify a model for the standard parameters. Alternatively, Almirall et al. expressed the standard parameter in terms of the blip effects of treatments and the effects of covariates [3]. However, it is highly difficult to specify models for the covariate effects.

Now consider the mean μ(h_t, z_t) = E(Y∣h_t, z_t) (t = 1, 2, ⋯ , T); for t = T, it is the standard parameter. Then the point effect of treatment z_t in stratum h_t is defined as

Notably, the treatment variable Z_t does not have posttreatment variables. Thus, this point effect refers to the effect of treatment in single-point causal inference and can be estimated by available methods [5, 13–15].

Under the identifying condition as Robins [1, 6], Wang and Yin expressed the point effect in terms of the blip effects [2, 16]. In designed experiments, where the treatment assignment condition is known and may reduce the dimensionality of the point effects, they demonstrated that the blip effect can be estimated via the point effects. In observational studies, where the treatment assignment is unknown, we will express the blip effect in terms of the standardized point effects and estimate and test the blip effect via the standardized point effects, as will be shown below.

Standardization is a procedure in epidemiology for estimating the point effect in strata of interest in the population. Here we extend it to a sequence of treatments. We divide the population at time t into a small number of disjoint strata. Let S_t indicate strata consisting of h_t at time t. The standardized point effect of treatment z_t in stratum S_t is defined as

where the expectation is with respect to P(H_t = h_t∣S_t). The set of the standardized point effects Θ(S_t; z_t) is far smaller than that of the point effects θ(h_t; z_t). Clearly, Θ(h_t; z_t) = θ(h_t; z_t) in the case of S_t = h_t.

In the Appendix, we prove

Theorem 1. Under the identifying condition, the standardized point effect is expressed in terms of the blip effects by

where the first expectation is with respect to P(H_t = h_t∣S_t), the second one to P(Xt+1s=xt+1s,Zt+1s-1=zt+1s-1,Zs=zs∣ht,zt) P(H_t = h_t∣S_t) , and the third one to P(Xt+1s=xt+1s,Zt+1s-1=zt+1s-1,Zs=zs∣ht,0) P(h_t∣S_t). Notably, the equation is Θ(S_T; z_T) = E{ϕ(h_T; z_T)∣S_T} at t = T.

Equation 2 is true without any treatment assignment conditions and thus applicable to observational studies. It is interesting to examine Equation 2 in the presence of a certain treatment assignment condition. Generally, if the treatment assignment satisfies P(Z_t∣h_t) = P(Z_t∣S_t), then we have the point effect of z_t in S_t given by θ(S_t; z_t) = μ(S_t, z_t)−μ(S_t, 0), where μ(S_t, z_t) = E(Y∣S_t, z_t).

In the Appendix, we prove

Theorem 2. Suppose that treatment assignment satisfies the condition P(Z_t∣h_t) = P(Z_t∣S_t) besides the identifying condition. Then the point effect is expressed in terms of the blip effects by

where the first expectation is with respect to P(H_t = h_t∣S_t, z_t), the second one to P(H_s = h_s, Z_s = z_s∣S_t, z_t), and the third one to P(H_s = h_s, Z_s = z_s∣S_t, 0). Notably, the equation is θ(S_T; z_T) = E{ϕ(h_T; z_T)∣S_T, z_T} at t = T.

Equation 3 can be used with a certain treatment assignment condition, applicable only for designed experiments. Sometimes in observational studies, the treatment assignment condition can be approximated by subclassification [13], and then, Equation 3 is also applicable.

First, we specify strata S_t (t = 1, …, T) in accordance with the blip effects of interest, aiming to improve the estimation and hypothesis testing of the blip effects of interest (i.e., the targeted analysis of blip effects). For instance, if we aim to analyze the modification of the blip effect by age and sex, then we specify strata by disjoint ranges of age and sex at times t = 1, …, T. Second, we estimate the point effect θ(h_t; z_t) in stratum S_t by available methods, for instance, modeling the mean μ(h_t, z_t). Finally, we estimate the standardized point effect in stratum S_t according to Equation 1,

where the expectation is with respect to the distribution P(H_t = h_t∣S_t).

For a finite sample, the expectation becomes a sum of the point effect with respect to the probability distribution P(H_t = h_t∣S_t) if h_t is discrete, or with respect to the probability distribution P(H_t = h_t∣S_t) dν(h_t) if h_t is continuous, where ν(.) is a proper measure of H_t. In both cases, the probability is estimated by the corresponding proportion denoted by P^(Ht∣St). Together with θ^(ht;zt), we obtain Θ^(St;zt).

There is considerable flexibility in estimating the standardized point effect Θ(S_t; z_t). Besides the usual regression, one can also estimate the point effect and the standardized counterpart by applying the propensity score-based subclassification method [13], the doubly robust method [14, 15], and others.

While it is straightforward to estimate Θ(S_t; z_t), it may not be easy to estimate the correlations between Θ^(St;zt) at different times. In contrast, the correlations between point effects θ^(ht;zt) at different times are negligible conditional on all treatments and covariates {h_T, z_T}, for instance, they are equal to zero for normally distributed outcomes [2, 16], so are the correlations between Θ^(St;zt) at different times conditional on all treatments and covariates {h_T, z_T}.

In most practices, the blip effects follow a certain pattern described by SNMM. Because the blip effect is a linear effect, that is, the difference in mean between potential outcomes, SNMM is often of the form

which is indexed by the blip effect vector γ = (γ₁, …, γ_k) of small dimension, where it is required that f_j(h_t, z_t = 0) = 0. For instance, ϕ(h_t; z_t) = γ_tz_t (t = 1, ..., T), in which f_t(h_t, z_t) = z_t and f_s(h_t, z_t) = 0 if s ≠ t. In this case, the blip effects at different times are different and can be estimated recursively at t = T, ..., 1 [5, 12]. However, the blip effects may be the same across different times, for instance, ϕ(h_t; z_t) = γ₁z_t, in which f₁(h_t, z_t) = z_t. It is far more difficult to estimate γ₁ than γ_t [3, 12].

With SNMM, Equation 2 becomes

where b_j(S_T, z_T) = E{f_j(h_T, z_T)∣S_T}, and for t = 1, …, T − 1,

Here, the first expectation is with respect to P(H_t = h_t∣S_t), the second to P(Xt+1s=xt+1s,Zt+1s-1=zt+1s-1,Zs=zs∣ht,zt) P(H_t = h_t∣S_t), and the third to P(Xt+1s=xt+1s,Zt+1s-1=zt+1s-1,Zs=zs∣ht,0) P(H_t = h_t∣S_t). For a finite sample, these expectations become summations, where the probabilities are P(H_t∣S_t) and P(Xt+1s,Zt+1s-1,Zs∣ht,zt) if covariates and treatments are discrete, or P(H_t∣S_t) dν(h_t) and P(Xt+1s,Zt+1s-1,Zs∣ht,zt) dν(xt+1s,zt+1s-1,zs) if covarites and treatments are continuous. In both cases, the probabilities are estimated by the corresponding proportions denoted by P^(Ht∣St) and P^(Xt+1s,Zt+1s-1,Zs∣ht,zt). This equation does not need any treatment assignment conditions and thus can be applied to observational studies.

Conditional on all treatments and covariates {h_T, z_T} prior to outcome Y, we apply Equation 4 as a regression model to estimate the blip effect vector γ. The Θ(S_t; z_t) has been estimated in Section 3.1. The conditional correlation between Θ^(St;zt) at diffrent times is approximately zero, as described in Section 3.1. The probabilities are estimated by the corresponding proportions, which are subject to no variability. The regression yields γ^. However, the covariates and treatments {h_T, z_T} are not ancillary to γ. Therefore, we need to incorporate the variability of {h_T, z_T} when estimating the covariance matrix cov(γ^). This can be achieved by, say, the bootstrap method.

Notably, the estimands Θ(S_t; z_t), E{f_j(h_t, z_t)∣S_t}, E{f_j(h_s, z_s)∣S_t, z_t}, and E{f_j(h_s, z_s)∣S_t, 0} are those in the framework of causal inference for single-point treatments. The conditions for the consistency and asymptotic normality of their estimates are well studied and often assumed to be satisfied in practice; see, for instance, Rosenbaum and Rubin [13] and Hernan and Robins [5]. Therefore, if these estimates are consistent and asymptotically normal, it follows from Equation 4 that γ^ is consistent and asymptotically normal. With γ^ and cov(γ^), we perform the Wald test on γ.

Now, suppose the treatment assignment condition P(Z_t∣h_t) = P(Z_t∣S_t). With SNMM, Equation 3 becomes

where b_j(S_T, z_T) = E{f_j(h_T, z_T)∣S_T, z_t}, and for t = 1, …, T − 1,

Here, E{f_j(h_t, z_t)∣S_t, z_t} is with respect to P(H_t = h_t∣S_t, z_t), E{f_j(h_s, z_s)∣S_t, z_t} to P(H_s = h_s, Z_s = z_s∣S_t, z_t), and E{f_j(h_s, z_s)∣S_t, 0} to P(H_s = h_s, Z_s = z_s∣S_t, 0). This equation needs treatment assignment conditions and thus applies only to designed studies.

In analogous to Equation 4, we apply Equation 5 to estimate and test γ. Conditional on all treatments and covariates {h_T, z_T} prior to outcome Y, we use this equation as a regression model to estimate the blip effect vector γ. The θ(S_t; z_t) is estimated by available methods, for instance, calculating the difference between the averages of the outcome in strata {S_t, z_t} and {S_t, 0}. The conditional correlation between θ^(St;zt) at different times is approximately zero, as described in Section 3.1. The probabilities are estimated by the corresponding proportions, which are subject to no variability. The regression yields γ^. We apply the bootstrap method to estimate the covariance matrix cov(γ^) incorporating the variability of {h_T, z_T}. If the estimates E^{fj(ht,zt)∣St,zt}, E^{fj(hs,zs)∣St,zt}, and E^{fj(hs,zs)∣St,0} are consistent and asymptotically normal as often assumed in causal inference for single-point treatments, so is γ^. With γ^ and cov(γ^), we may conduct the Wald test on γ.

Several statements can be made comparing the two methods based on Equations 4 and 5. First, the former method is applicable to observational studies, whereas the latter is applicable only to designed studies. Second, both methods reduce the dimension of point effects, simplifying the regression. Third, the former allows for standardizing the point effects in accordance with the blip effects of interest, leading to a targeted analysis of blip effects, in comparison to the latter. For instance, a complete randomized trial of treatments in sequence does not allow for an analysis of the modification of blip effects by any covariates. Finally, the designed studies are still superior due to the fact that the identifying conditions are satisfied by design in comparison to observational studies.

Suppose a treatment sequence of length T = 3. The treatment variables Z_t are dichotomous with z_t = 0, 1 (t = 1, 2, 3). The time-dependent covariates X_t are polytomous with x_t = 0, 1, 2, 3 (t = 2, 3). Conditional on all treatment and time-dependent covariates {z₁, x₂, z₂, x₃, z₃}, the outcome Y follows the normal, Bernoulli, or Poisson distribution. A summary of the variables is {Z₁, X₂, Z₂, X₃, Z₃, Y} in the temporal order, with observations {z₁, x₂, z₂, x₃, z₃, y}. In the Supplementary material, we describe the simulation study in detail. The relevant codes developed for the simulation are given by Yin [17].

From the treatment sequence, we have 73 = 1 + 8 + 64 point effects of treatments,

We divide the population at time t = 2 into strata {x₂ = j} (j = 0, 1, 2, 3) and the population at time t = 3 into strata {x₃ = j} (j = 0, 1, 2, 3). Then we may obtain the following nine standardized point effects,

Our method, based on Equation 4, is applicable when estimating the blip effects without treatment assignment conditions. To evaluate the performance of this method, however, we will compare it with a method based on Equation 5, which is only applicable with a certain treatment assignment condition. Thus, we choose a strong treatment assignment such that the assignment of z₂ depends only on x₂ and that of z₃ only on x₃; however, the simulation can readily be extended to weaker treatment assignment conditions. Thus, the treatment assignment satisfies the condition P(Z_t∣h_t) = P(Z_t∣x_t) (t = 2, 3), so we have the following nine standardized point effects,

Suppose two specific SNMMs, SNMM1 and SNMM2, are of the following form. With SNMM1, the blip effects are

Thus, for SNMM1, we have γ = (γ₁, γ₂₀, γ₂₁, γ₂₂, γ₂₃, γ₃₀, γ₃₁, γ₃₂, γ₃₃). With SNMM2, it is further required that γ_2j = γ_3j, so there are only four different blip effects in addition to γ₁. Thus for SNMM2, we have γ = (γ₁, γ₂₀, γ₂₁, γ₂₂, γ₂₃).

Under SNMM1 or SNMM2, we will estimate and test the blip effect by applying Equation 4 to Θ(z₁ = 1), Θ(x₂; z₂ = 1), and Θ(x₃; z₃ = 1) or by applying Equation 5 to θ(z₁ = 1), θ(x₂; z₂ = 1), and θ(x₃; z₃ = 1). Specifically, method (i) uses Equation 5 and SNMM1, with treatment assignment condition. Method (ii) uses Equation 4 and SNMM1, without a treatment assignment condition. Method (iii) uses Equation 5 and SNMM2, with treatment assignment condition. Method (iv) uses Equation 4 and SNMM2, without a treatment assignment condition. With methods (i)–(iv), we obtain the estimate and variance for the estimated blip effect as well as the coverage probability and the power for the hypothesis testing of the blip effect. The result is presented in Table 1a.

From Table 1b, the following observations can be made for estimating and testing the blip effects. First, all four methods achieve unbiased estimates and the nominal level of the coverage probability. Second, as compared to methods (i) and (ii) with SNMM1, methods (iii) and (iv) may impose SNMM2 across times t = 2, 3, reducing the number of blip effects and thus resulting in a smaller variance and a greater power for estimating and testing the blip effects. Third, methods (ii) and (iv) are based on Equation 4 and achieve nearly the same results as methods (i) and (iii) based on Equation 5, demonstrating that our method performs equally well with or without treatment assignment condition.

As described above, our method (iv) uses SNMM2 and Equation 4 when estimating and testing the blip effect without any treatment assignment conditions. Here, we compare this method with the following three available methods in the literature. Method (v) is the parametric method based on the well-known G-formula expressing the blip effect in terms of the standard parameters [5, 7]. Method (vi) is the marginal structural model based on the inverse probability of treatment weighting [5, 9]. Method (vii) is the G-estimation based on SNMM1 [1, 5, 11, 12]. These methods are reviewed in the introduction and also described in the context of the simulation in the Supplementary material.

With methods (iv)–(vii), we obtain the estimate and variance for the blip effect as well as the coverage probability and power for the hypothesis testing of the blip effect. The result is presented in Table 2. As seen from Table 2, all four methods achieve unbiased estimates and the nominal level of the coverage probability. Our method (iv) achieves the smallest variance and largest power due to SNMM2. Methods (v) and (vi) yield sizable variances and low powers due to the difficulty of imposing any SNMMs. Method (vii) yields a smaller variance and a larger power than methods (v) and (vi) due to SNMM1, but a larger variance and smaller power than method (iv) due to the difficulty of imposing SNMM2. In general, it is difficult to introduce SNMMs across different times with the G-estimation.

It is also interesting to compare our method (iv) with method (vii) in terms of other modeling conditions than SNMM. With method (iv), we only need models for the point effects, which in this simulation are the three models for μ(z₁), μ(z₁, x₂, z₂), and μ(z₁, x₂, z₂, x₃, z₃) plus SNMM2. Standardization does not need additional models. With method (vii), we need the following models. First, a model for μ(x₃, z₃), which is smaller than μ(z₁, x₂, z₂, x₃, z₃) due to the treatment assignment condition. Second, a model for the baseline E{Y(D₂ = 0, D₃ = 0)∣x₂} together with SNMM1 at t = 2. Third, a model for baseline E{Y(D₁ = 0, D₂ = 0, D₃ = 0) together with SNMM1 at t = 1. In general, the baseline E{Y(DtT=0)∣ht} with the G-estimation is typically subject to model misspecification in comparison to the baseline E(Y∣h_t, z_t = 0) with our method; for this reason, the doubly robust version of the G-estimation may be used [5, 12]. Although it is not needed in the simple setting of this article, we may have a doubly robust version of our method by obtaining the doubly robust estimates for the point effects.

In Sweden, patients usually seek medical help at hospitals near their residential areas. When cancer is diagnosed, they may stay at the diagnosing hospital or transfer to another hospital for treatment. The hospital that diagnoses cancer is called the diagnosing hospital, while the one that treats cancer is called the treating hospital. To evaluate the performance of diagnosing and treating hospitals, one may study the blip effects of diagnosing and treating hospitals among cancer patients after adjusting for patients' differences.

The data used in this study contain the information on 1, 070 stomach cancer patients from a clinical study during the period between 1988 and 1995 in hospitals located in central and northern Sweden [18]. Stomach cancer is highly malignant with a poor prognosis, so the 1-year survival is a good measure of the performance of both diagnosing and treating hospitals. A question of medical relevance is which types of diagnosing and treating hospitals, large vs. small, perform better on cancer outcomes, where the large type refers to the regional or county hospitals and the small type to local hospitals. One concern is that young patients diagnosed at local hospitals tend to have poor prognoses. This phenomenon is known as doctors' delay in the area of cancer diagnosis, but little studied statistically [19].

The diagnosing hospital is the treatment variable Z₁ at time t = 1: z₁ = 0 for small type and z₁ = 1 for large type. The treating hospital is the treatment variable Z₂ at t = 2: z₂ = 0 for small type and z₂ = 1 for large type. The outcome of interest is Y: y = 1 for a successful 1-year survival and y = 0 otherwise. The stationary covariates X₁ = (X₁₁, X₁₂, X₁₃) before Z₁ were measured with gender X₁₁, geographic area X₁₂ and age X₁₃. Gender was x₁₁ = 0 for female and x₁₁ = 1 for male. Geographic area was categorized into rural x₁₂ = 0 vs. urban x₁₂ = 1. Age took continuous values x₁₃. The time-dependent covariate between Z₁ and Z₂ was cancer stage X₂, taking the values x₂ = 1, 2, 3, 4 for cancer stages. The descriptive statistics are given in Table 3. The data and code are given by Yin [17].

Notably, the data are unbalanced between diagnosing and treating hospitals: 981 patients did not transfer, 80 of them from small diagnosing hospitals to large ones, and only 6 from large to small diagnosing hospitals. According to Swedish medical experts, it was a typical transferal pattern. The small sample size and unbalanced data contribute to the small significance of our results. The confounding situation is described below.

Due to the long-term social welfare system and the 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 blip effect. As a common practice in many epidemiological research studies in Sweden, the assumption of no unmeasured confounders is approximately satisfied for diagnosing hospital Z₁, at least after conditioning on gender x₁₁, residential area x₁₂, and age x₁₃. Similarly, the assumption of no unmeasured confounders is also approximately satisfied for treating hospital Z₂ conditional on the cancer stage x₂ and the diagnosing hospital z₁ in addition to (x₁₁, x₁₂, x₁₃). The causal DAG for diagnosing and treating hospital types is given in Figure 1.

Without further examining the validity of the identifying condition, we will focus on the inference part of the medical example; interested readers may find a large body of literature on sensitivity analysis of the causal effect to the identifying condition (e.g., [20, 21]).

Because we aim at the effect modification of the blip effect by age, we divide the population at t = 1 into two strata U and L, where S₁ = U is the one with age x₁₃ smaller than the median and S₁ = L with age x₁₃ larger than or equal to the median. So, we have two standardized point effects Θ(S₁; z₁ = 1).

A large variety of methods are available for estimating the point effect in the framework of causal inference for single-point treatments (e.g., [5, 13–15]). As an illustration, we use the usual regression to estimate the point effects of diagnosing hospital z₁. Because the sample is small and only three stationary covariates—that is, x₁₁, x₁₂, x₁₃—are involved in the estimation, we model the conditional mean μ(x₁₁, x₁₂, x₁₃, z₁) = E(Y∣x₁₁, x₁₂, x₁₃, z₁) in the whole data set. For the sake of presentation, we use the linear model to estimate the point effect. Notably, we may use the logistic model, which only improves the estimation slightly.

We exclude the residential area x₁₂ at a significance level of 0.05, consistent with medical observations indicating that the residential area is less influential than gender and age. Finally, we obtain the regression model,

From this model, we obtain the point effects of z₁ = 1 in stratum (x₁₁, x₁₃),

which is the same for all (x₁₁, x₁₃). The standardized point effects are obtained by Θ(S₁; z₁ = 1) = E{θ(x₁₁, x₁₃; z₁ = 1)∣S₁} with respect to P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁), where the point effect is estimated above while the probability is estimated by the corresponding proportion. Clearly, Θ(S₁; z₁ = 1) = θ₁ from Model 8. Therefore, we have the estimate Θ^(S1;z1=1)=θ^1, but the variance of Θ^(S1;z1=1) is obtained by adjusting the variance of θ^1 to the size of sub sample S₁.

To estimate the point effects of treating hospital z₂, we model the mean μ(x₁₁, x₁₂, x₁₃, z₁, x₂, z₂). We exclude x₁₂ at the significance level of 0.05. Furthermore, age x₁₃ has rather different influences on cancer survival for different cancer stages x₂, so we model the conditional mean separately for different x₂. Finally, we obtain the regression model,

From this model, we obtain the point effects of z₂ = 1 in stratum (x₁₁, x₁₃, x₂) for early cancer stage x₂ = 1,

Averaging it with respect to P(X₁₁ = x₁₁, X₁₃ = x₁₃∣x₂ = 1),we obtain the average point effect of z₂ = 1 in stratum x₂ = 1,

In cancer stage x₂ = j = 2, 3, 4, we have,

Because the point effects are equal to the blip effects for the last treatment z₂ = 1, we use the point effects of z₂ = 1 as a special case of the standardized point effects in estimating and testing the blip effects of z₁ = 1 and z₂ = 1.

As seen from Models 8 and 9, residential area x₁₂ does not appear in μ(x₁₁, x₁₃, z₁) and μ(x₁₁, x₁₃, x₂, z₂) and thus are irrelevant to the standardized point effects and the blip effects of treatments. Hence, we remove x₁₂ from the following development.

The point effect of diagnosing hospital results from both diagnosing and treating hospitals and thus cannot be used to evaluate the diagnosing hospital. In comparison, the blip effect of the diagnosing hospital represents the causal effect of the diagnosing hospital, while setting the treating hospitals as small ones, and thus can be used for evaluation.

To study the phenomenon of doctors' delay in cancer diagnosis [19], we suppose that the blip effect of large diagnosing hospital ϕ(x₁₁, x₁₃; z₁ = 1) is a linear function of age x₁₃, that is, age modifies the blip effect. Furthermore, because Z₂ is the last treatment variable in the treatment sequence, the blip effect of z₂ is equal to the point effect of z₂ as obtained from Model 9. Summarizing these observations, we specify an SNMM of the following form,

Thus, SNMM 8 is indexed by γ = (γ_{1, 0}, γ_{1, 3}, γ_{21, 0}, γ_{21, 3}, γ₂₂, γ₂₃, γ₂₄). Let γ₁ = E{ϕ(x₁₁, x₁₃; z₁ = 1)} be the blip effect of a large diagnosing hospital z₁ = 1 in the whole population. Then we have γ₁ = γ_{1, 0} + E(X₁₃)γ_{1, 3}, so γ_{1, 3} is the modification of the blip effect γ₁ by age. Let γ₂₁ = E{ϕ(x₁₁, x₁₃, z₁, x₂; z₂ = 1)∣x₂ = 1} be the blip effect of large treating hospital z₂ = 1 in low cancer stage x₂ = 1. Then we have γ₂₁ = γ_{21, 0} + E(X₁₃∣x₂ = 1)γ_{21, 3}, so γ_{21, 3} is the modification of the blip effect γ₂₁ by age in low cancer stage x₂ = 1. The blip effect of z₂ = 1 in cancer stage x₂ = j = 2, 3, 4 is γ_2j.

Because the blip effect ϕ(x₁₁, x₁₃; z₁ = 1) is indexed by two parameters γ_{1, 0} and γ_{1, 3} under SNMM, as expressed in Equation 8, we need two standardized point effects Θ(S₁; z₁ = 1) = E{θ(x₁₁, x₁₃; z₁ = 1)∣S₁} with respect to P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁) for S₁ = U, L, as described and estimated in Section 5.2. Now by applying Equation 2 or Equation 4 to Θ(S₁; z₁ = 1), we obtain

Here the first expectation is with respect to the P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁), so we have

The second expectation is with respect to P(X₂ = x₂, Z₂ = z₂∣x₁₁, x₁₃, z₁ = 1) P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁) and the third one to P(X₂ = x₂, Z₂ = z₂∣x₁₁, x₁₃, z₁ = 0) P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁). By using ϕ(x₁₁, x₁₃, z₁, x₂; z₂ = 0) = 0, the second and third expectations are, for z₁ = 1, 0,

with respect to P(X₁₁ = x₁₁, X₁₃ = x₁₃∣S₁). Now by inserting SNMM, as provided in Equation 8, we obtain