Individualized treatment rules under stochastic treatment cost constraints

4.1 Pathwise differentiability of the ATE

We first present a result regarding the pathwise differentiability of the ATE. Pathwise differentiability of the parameter of interest serves as the foundation for constructing asymptotically efficient estimators of this parameter, based on which an inferential procedure may be developed. Additional technical conditions are required and are provided in Section S1 in the Supplementary Material. For a distribution P ∈ ℳ , a function μ C : { 0 , 1 } × W → R , an ITR ρ , and a decision threshold τ ∈ R , we define pointwise the following functions:

One key condition we rely on is the following nonexceptional law assumption.

Under this condition, the true optimal ITR ρ 0 is identical to an indicator function. If all covariates are discrete, then we can plug in the empirical estimates into the identification formulae in Theorems 1 and 2 and show that the resulting estimators of the ATE are asymptotically normal by the delta method even when Condition B1 does not hold. We do not further pursue this simple case in this article, and thus need to rely on the nonexceptional law assumption, namely Condition B1, to account for continuous covariates. We list additional technical conditions in Supplement S1.

We can now provide a formal result describing the pathwise differentiability of the ATE parameter.

We note that the pathwise differentiability of P ↦ Ψ ρ P ℛ ( P ) was established in Theorem 3 in the study by Qiu et al. [9] for ℛ ∈ { FR , TP } . The other results can be proven using similar techniques. We presented the proof of these results in Supplement S4.2. In view of Theorem 3, it follows that the ATE parameter Ψ ℛ is pathwise differentiable at P 0 with nonparametric canonical gradient:

for ℛ ∈ { FR , RD , TP } .

4.2 Proposed estimator and asymptotic linearity

We next present our proposed nonparametric procedure for estimating an optimal ITR ρ 0 and the corresponding ATE ψ 0 . We will generally use subscript n to denote an estimator with sample size n and add a hat to a nuisance function estimator that is targeted toward estimating ϕ 0 .

1. Use the empirical distribution P ˆ W , n of W as an estimate of the true marginal distribution of W . Compute estimates μ n Y , μ n C , μ n T , δ n Y and δ n C of μ 0 Y , μ 0 C , μ 0 T , δ 0 Y and δ 0 C , using flexible regression methods. Recall that μ 0 Y ( t , w ) = E 0 [ Y ∣ T = t , W = w ] , μ 0 C ( t , w ) = E 0 [ C ∣ T = t , W = w ] , δ 0 Y ( v ) = E 0 [ μ 0 Y ( 1 , W ) − μ 0 Y ( 0 , W ) ∣ V = v ] , and δ 0 C ( v ) = E 0 [ μ 0 C ( 1 , W ) − μ 0 C ( 0 , W ) ∣ V = v ] . Define pointwise Δ n C ( w ) ≔ μ n C ( 1 , w ) − μ n C ( 0 , w ) .

2. Estimate an optimal ITR:

   1. Estimate ϕ 0 = E 0 [ μ 0 C ( 0 , W ) ] with a one-step correction estimator:

      ϕ n ≔ 1 n ∑ i = 1 n μ n C ( 0 , W i ) + ( 1 − T i ) ( C i − μ n C ( 0 , W i ) ) 1 − μ n T ( W i ) .

   2. Let ξ n ≔ δ n Y / δ n C , Γ n : τ ↦ 1 n ∑ i : ξ n ( V i ) > τ Δ n C ( W i ) and γ n : τ ↦ 1 n ∑ i : ξ n ( V i ) = τ Δ n C ( W i ) . For any k ∈ [ 0 , ∞ ] , define η n ( k ) ≔ inf { τ : Γ n ( τ ) ≤ k − α ϕ n } , τ n ( k ) ≔ max { η n ( k ) , 0 } and

      d n , k : v ↦ k − α ϕ n − Γ n ( η n ( k ) ) γ n ( η n ( k ) ) : if ξ n ( v ) = η n ( k ) and γ n ( η n ( k ) ) > 0 , I { ξ n ( v ) > η n ( k ) } : otherwise.

      The rule d n , k is the sample analog of an ITR that prescribes treatment to those with the highest values of ξ 0 ( V ) , regardless of whether the treatment is harmful, until treatment resources run out.

   3. Compute k n , which is used to define an estimate of ρ 0 for which the plug-in estimator is asymptotically linear under conditions, as follows:

      1. if τ n ( κ ) > 0 and there is a solution in k ∈ [ 0 , ∞ ) to

         (5) 1 n ∑ i = 1 n d n , k ( V i ) Δ n C ( W i ) + C i − μ n C ( T i , W i ) t i + μ n T ( W i ) − 1 + α ϕ n = κ ,

         then take k n to be this solution;

      2. otherwise, set k n = κ .

   4. Estimate ρ 0 using the sample analog of ρ 0 with treatment resource constraint k n , namely,

      ρ n : v ↦ k n − α ϕ n − Γ n ( τ n ( k n ) ) γ n ( τ n ( k n ) ) : if ξ n ( v ) = τ n ( k n ) and γ n ( τ n ( k n ) ) > 0 I { ξ n ( v ) > τ n ( k n ) } : otherwise.

3. Obtain an estimate ρ n ℛ of the reference ITR ρ 0 ℛ as follows:

   1. For ℛ = FR , take ρ n ℛ to be ρ FR .

   2. For ℛ = RD ,

      1. obtain a targeted estimate μ ˆ n C ( 1 , ⋅ ) of μ 0 C ( 1 , ⋅ ) : run an ordinary least-square linear regression with outcome C , covariate 1 / ( T + μ n T ( W ) − 1 ) , offset μ n C ( T , W ) , and no intercept. Take μ ˆ n C to be the fitted mean model;

      2. take ρ n ℛ to be the constant function w ↦ min 1 , ( κ − ϕ n ) / 1 n ∑ i = 1 n Δ ˆ n C ( W i ) , where we define pointwise Δ ˆ n C ( w ) ≔ μ ˆ n C ( 1 , w ) − μ ˆ n C ( 0 , w ) .

   3. For ℛ = TP , take ρ n ℛ to be μ n T .

4. Estimate ATE of ρ 0 relative to the reference ITR ρ 0 ℛ with a TMLE ψ n :

   1. obtain a targeted estimate μ ˆ n Y of μ 0 Y : run an ordinary least-square linear regression with outcome Y , covariate [ ρ n ( V ) − ρ n ℛ ( W ) ] / [ T + μ n T ( W ) − 1 ] , offset μ n Y ( T , W ) , and no intercept. Take μ ˆ n Y to be the fitted mean function.

   2. with P ˆ n being any distribution with components μ ˆ n Y and P ˆ W , n , take

      ψ n ≔ Ψ ρ n ( P ˆ n ) − Ψ ρ n ℛ ( P ˆ n ) = 1 n ∑ i = 1 n [ ρ n ( V i ) − ρ n , i ℛ ] [ μ ˆ n Y ( 1 , W i ) − μ ˆ n Y ( 1 , W i ) ] ,

      where ρ n , i ℛ is defined as ρ n ℛ ( W i ) or ρ n ℛ ( V i ) depending on the covariate used by the reference ITR.

The aforementioned procedure is similar to that proposed in the study by Qiu et al. [9]. One key difference is the use of the refined estimator k n of κ obtained via the estimating equation (5), which is a key to ensuring the asymptotic linearity of ψ n . Another difference is that the denominator of ξ n is now δ n C , which is consistent with our different definition of the unit value for solving the fractional knapsack problem (2). Similar to TMLE for other problems, when C or Y has known bounds (e.g., the closed interval [ 0 , 1 ] ), to obtain a corresponding targeted estimate that respect the known bounds, we may use logistic regression rather than ordinary least squares [36].

The aforementioned procedure has both similarities and substantial differences compared to the estimation procedure proposed by Sun et al. [11]. The main difference is that our procedure is targeted towards efficient estimation of and inference about the ATE of ψ 0 of the optimal ITR under a nonparametric model, while Sun et al. [11] focused on estimating the optimal ITR ρ 0 and does not evaluate this optimal ITR. This leads to a key difference between the two procedures when estimating the optimal ITR: we need to solve an estimating equation (5), which is crucial to ensuring that the estimator ψ n is asymptotically linear, while Sun et al. [11] do not. The requirement of solving (5) is related to the nature of the fractional knapsack problem discussed in Remark 3, and we conjecture that such a calibration on the resource used is necessary for general problems of the same nature. Our procedure is also related to the method in Sun [12]. Sun [12] relied on the availability of asymptotically normal estimators of both the average benefit and average resource used (Assumption 2.4), a nontrivial requirement when the propensity score μ 0 T is unknown in observational studies. Our procedure essentially produces such estimators: in Step 4, an asymptotically normal estimator of the ATE is constructed, whereas an asymptotically normal estimator of the expected resource is produced in Step 2 and used to calibrate the resource expenditure of the estimated optimal ITR ρ n in Step 2(c).

We now present results on the asymptotic linearity and efficiency of our proposed estimator. We state and discuss the technical conditions required by the theorem below in Supplement S1.

To conduct inference about ψ 0 , we can directly plug the estimators of nuisance functions into D ℛ ( P 0 ) to obtain a consistent estimator of D ℛ ( P 0 ) , and then take the sample variance to obtain a consistent estimator of the asymptotic variance σ 0 2 . The proof of Theorem 4 can be found in Supplements S4.3 and S4.4.

5.1 Simulation setting

In this simulation study, we investigate the performance of our proposed estimator of the ATE of an optimal ITR relative to specified reference ITRs. We focus here on the setting α = 1 . This scenario is more difficult than the case α = 0 because it requires the estimation of ϕ 0 .

We generate data from a model in which the treatment T is an IV and both the treatment cost C and outcome Y are binary. This data-generating mechanism satisfies all causal conditions and has an unobserved confounder between treatment cost and outcome. We first generate a trivariate covariate W = ( W 1 , W 2 , W 3 ) , where W 1 ∼ Unif ( − 1 , 1 ) , W 2 ∼ Bernoulli ( 0.8 ) , and W 3 ∼ N ( 0 , 1 ) are mutually independent. We also simulate an unobserved treatment-outcome confounder U ∼ Bernoulli ( 0.5 ) independently of W , and then simulate T , C , and Y as follows:

We introduce U in the data-generating mechanism to emphasize that we do not require assumptions on the joint distribution of treatment cost and outcome conditional on covariates. We consider all three reference ITRs ℛ ∈ { FR , RD , TP } , where we set ρ FR : v ↦ 0 . We set κ = 0.68 , which is an active constraint with τ 0 > 0 and ρ 0 RD < 1 .

The ITRs we consider are based on all covariates – that is, we take V ( W ) = W . We estimate the nuisance functions using the Super Learner [43] with library including a logistic regression, generalized additive model with logit link [44], gradient boosting machine [45,46, 47], support vector machine [48,49], and neural network [50,51]. Because none of the nuisance functions follows a logistic regression model, the resulting ensemble learner is not expected to achieve the parametric convergence rate. Since both C and Y are binary, we use logistic regression rather than ordinary least squares to obtain their corresponding targeted estimates in Section 4.2. We consider sample size n ∈ { 500 , 1,000 , 4,000 , 16,000 } and run 1000 Monte Carlo repetitions for each sample size. We implement the algorithm that incorporates cross-fitting discussed in Remark 6 and described in Section S3 in the Supplementary Material.

To evaluate the performance of our proposed estimator, we investigate the bias and root-mean-squared error (RMSE) of the estimator. We also investigate the coverage probability and the width of nominal 95% Wald CIs constructed using influence function-based standard error estimates. We further investigate the probability that our confidence lower limit falls below the true ATE, that is, the coverage probability of the 97.5% Wald confidence lower bound.

5.2 Simulation results

Table 1 presents the performance of our proposed estimator in this simulation. For sample sizes 500, 1,000 and 4,000, the CI coverage of our proposed method is lower than the nominal coverage 95%. When sample size is larger (16,000), the CI coverage of our proposed method increases to 90–93%. The coverage of the confidence lower bounds is much closer to nominal (97.5%) for all sample sizes considered, though, and is always approximately nominal when the sample size is large. For all reference ITRs, the bias and RMSE of our proposed estimator appear to converge to zero faster than and at the same rate as the square root of sample size, respectively. All biases are negative, which is expected in view of Remark 6. All standard errors underestimate the variation of the estimator with the extent decreasing as sample size increases.

Figure 1 presents the width of the Wald CIs scaled by the square root of sample size n . Our theory indicates that the CI width should shrink at a root- n rate, and our simulation results are consistent with this. There are some outlying cases of extremely wide or narrow CIs. This is expected for small sample sizes because the estimator of σ 0 2 in Theorem 4 resembles a sample mean and might not be close to σ 0 2 with high probability when the sample size is small. In practice, this issue might be slightly mitigated by fine-tuning the involved machine learning algorithms.

Figure 1

Boxplot of n × CI width for ATE relative to each reference IER.

As indicated in Theorem 4, theoretical guarantees for the validity of the Wald CIs rely on the nuisance function estimators converging to the truth sufficiently quickly. It appears that the undercoverage of our Wald CI in small samples may owe, in part, to poor estimation of these nuisance functions in small sample sizes. To illustrate how our procedure may perform with improved small-sample nuisance function estimators, we conducted another two simulations: one is identical to those reported earlier in all ways except that the nuisance function estimators μ n Y , μ n C , and μ n T are taken to be equal to the truth; the other is a simpler scenario under a lower dimension and a parametric model. The results are presented in Section S5 in the Supplementary material and suggest that our proposed estimator may achieve significantly better performance with improved machine learning estimators of the nuisance functions. This motivates seeking ways to optimize the finite-sample performance of the nuisance function estimators employed in future applications of the proposed method, possibly based on prior subject-matter expertise. The underestimation of standard errors in this simulation also motivates future work exploring whether there are standard error estimators with better finite-sample performance, for example, estimators based on the bootstrap.