Optimal Individualized Treatments in Resource-Limited Settings

1 Introduction

Suppose one wishes to maximize the population mean of some outcome using some binary point treatment, where for each individual clinicians have access to (some subset of) measured baseline covariates. Such a treatment strategy is termed an individualized treatment regime (ITR), and the (counterfactual) population mean outcome under an ITR is referred to as the value of an ITR. The ITR which maximizes the value is referred to as the optimal ITR or the optimal rule. There has been much recent work on this problem in the case where treatment is an unlimited resource (see Murphy [1] and Robins [2] for early works on the topic, and Chakraborty and Moodie [3] for a recent overview). It has been shown that the optimal treatment in this context is given by checking the sign of the average treatment effect conditional on (some subset of) the baseline covariates, also known as the blip function [2].

The optimal ITR assigns treatment to people from a given strata of covariates for which treatment is on average beneficial, and does not assign treatment to this strata otherwise. If treatment is even slightly beneficial to all subsets of the population, then such a treatment strategy would suggest treating the entire population. There are many realistic situations in which such a treatment strategy, or any strategy that treats a large proportion of the population, is not feasible due to limitations on the total amount of the treatment resource. In a discussion of Murphy [1], Arjas observed that resource constraints may render optimal ITRs of little practical use when the treatment of interest is a social or educational program, though no solution to the constrained problem was given [4].

The mathematical modeling literature has considered the resource allocation problem to a greater extent. Lasry et al. [5] developed a model to allocate the annual CDC budget for HIV prevention programs to subpopulations which would benefit most from such an intervention. Tao et al. [6] consider a mathematical model to optimally allocate screening procedures for sexually tranmitted diseases subject to a cost constraint. Though Tao et al. do not frame the problem as a statistical estimation problem, they end up confronting similar optimization challenges to those that we will face. In particular, they confront the (weakly) NP-hard knapsack problem from the combinatorial optimization literature [7, 8]. We will end up avoiding most of the challenges associated with this problem by primarily focusing on stochastic treatment rules, which will reduce to the easier fractional knapsack problem [9, 8]. Stochastic ITRs allow the treatment to rely on some external stochastic mechanism for individuals in a particular strata of covariates.

We consider a resource constraint under which there is a maximum proportion of the population which can be treated. We primarily focus on evaluating the public health impact of an optimal resource-constrained (R-C) ITR via its value. The value function has been shown to be of interest in several previous works (see, e.g., Zhang et al. [10], van der Laan and Luedtke [11], Goldberg et al. [12]. Despite the general interest of this quantity, estimating this quantity is challenging for unconstrained deterministic regimes at so-called exceptional laws, i.e. probability distributions at which the blip function is zero in some positive probability strata of covariates [2]; a slightly more general assumption is given in Luedtke and van der Laan [13]. Chakraborty et al. [14] showed that one can develop confidence intervals for this parameter using m-out-of-n bootstrap, though these confidence intervals shrink at a slower than root-n rate. Luedtke and van der Laan [13] showed that root-n rate confidence intervals can be developed for this quantity under reasonable conditions in the large semiparametric model which at most places restrictions on the treatment mechanism.

We develop a root-n rate estimator for the optimal R-C value and corresponding confidence intervals in this same large semiparametric model. We show that our estimator is efficient among all regular and asymptotically linear estimators under conditions. When the baseline covariates are continuous and the resource constraint is active, i.e. when the optimal R-C value is less than the optimal unconstrained value, these conditions are far more reasonable than the non-exceptional law assumption needed for regular estimation of the optimal unconstrained value.

We now give a brief outline of the paper. Section 2 defines the statistical estimation problem of interest, gives an expression for the optimal deterministic rule under a condition, and gives a general expression for the optimal stochastic rule. Section 3 presents our estimator of the optimal R-C value. Section 4 presents conditions under which the optimal R-C value is pathwise differentiable, and gives an explicit expression for the canonical gradient under these conditions. Section 5 describes the properties of our estimator, including how to develop confidence intervals for the optimal R-C value. Section 6 presents our simulation methods. Section 7 presents our simulation results. Section 8 closes with a discussion and areas of future research. All proofs are given in the Appendix.

2 Optimal R-C rule and value

Suppose we observe n independent and identically distributed (i.i.d.) draws from a single time point data structure (W,A,Y)∼P0, where the vector of covariates W has support W, the treatment A has support {0,1}, and the outcome Y has support in the closed unit interval. Our statistical model is nonparametric, beyond possible knowledge of the treatment mechanism, i.e. the probability of treatment given covariates. Little generality is lost with the bound on Y, given that any continuous outcome bounded in [b, c] can be rescaled to the unit interval with the linear transformation (y−b)/(c−b). Suppose that treatment are resources are limited so that at most a κ∈(0,1) proportion of the population can receive the treatment A=1. Let V be some function of W, and denote the support of V with V. A deterministic treatment rule d˜ takes as input a function of the covariates v∈V and outputs a binary treatment decision d˜(v). The stochastic treatment rules considered in this work are maps from U∈V to {0,1}, where U is the support of some random variable U∼PU. If d is a stochastic rule and u∈U is fixed, then d(u,⋅) represents a deterministic treatment rule. Throughout this work we will let U be drawn independently of all draws from P0.

For a distribution P, let QˉP(a,w)=ΔEP[Y|A=a,W=w]. For notational convenience, we let Qˉ0≜QˉP0. Let d˜ be a deterministic treatment regime. For a distribution P, let Ψ˜d˜=ΔEP0[QˉP(d˜(V),W)] represent the value of d˜. Under causal assumptions, this quantity is equal to the counterfactual mean outcome if, possibly contrary to fact, the rule d˜ were implemented in the population [15, 16]. The optimal R-C deterministic regime at P is defined as the deterministic regime d˜ which solves the optimization problem

For a stochastic regime d, let ΨdP=ΔEPUΨ˜dU,.P represent the value of d. Under causal assumptions, this quantity is equal to the counterfactual mean outcome if, possibly contrary to fact, the stochastic rule d were implemented in the population (see Ref. [17] for a similar identification result). The optimal R-C stochastic regime at P is defined as the stochastic treatment regime d which solves the optimization problem

We call the optimal value under a R-C stochastic regime Ψ(P). Because any deterministic regime can be written as a stochastic regime which does not rely on the stochastic mechanism U, we have that Ψ(P)≥Ψ˜(P). The constraint EPU×P[d(U,V)]≤κ above is primarily meant to represent a clinical setting where each patient arrives at the clinic with covariate summary measure V, a value of U is drawn from PU for this patient, and treatment is then assigned according to d(U,V). By Fubini’s theorem, this is like rewriting the above constraint as EPEPU[d(U,V)]≤κ. Nonetheless, this constraint also represents the case where a single value of U=u is drawn for the entire population, and each individual is treated according to the deterministic regime d(u,⋅), i.e. EPUEP[d(U,V)]≤κ. This case appears less interesting because, for a fixed u, there is no guarantee that EP[d(u,V)]≤κ.

For a distribution P, define the blip function as

Let SP represent the survival function of Qˉb,P, i.e. τ↦PrP(Qˉb,P>τ). Let

For notational convenience we let Qˉb,0=ΔQˉb,P0, S0=ΔSP0, η0=ΔηP0, and τ0=ΔτP0.

Define the deterministic treatment rule d˜P as v↦I(Qˉb,P(v)>τP), and for notational convenience let d˜0=Δd˜P0. We have the following result.

One can in fact show that d˜P is the P almost surely unique optimal deterministic regime under the stated condition. We do not treat the case where PrP(Qˉb,P(V)=τP)>0 for deterministic regimes, since in this case (1) is a more challenging problem: for discrete V with positive treatment effect in all strata, eq. (1) is a special case of the 0–1 knapsack problem, which is NP-hard, though is considered one of the easier problems in this class [7, 8]. In the knapsack problem, one has a collection of items, each with a value and a weight. Given a knapsack which can only carry a limited weight, the objective is to choose which items to bring so as to maximize the value of the items in the knapsack while respecting the weight restriction. Considering the optimization problem over stochastic rather than deterministic regimes yields a fractional knapsack problem, which is known to be solvable in polynomial time [9, 8]. The fractional knapsack problem differs from the 0–1 knapsack problem in that one can pack partial items, with the value of the partial items proportional to the fraction of the item packed.

Define the stochastic treatment rule dP by its distribution with respect to a random variable drawn from PU:

We will let d0=ΔdP0. Note that d˜P(V) and dP(U,V) are PU×P almost surely equal if PrP(Qˉb,P(V)=τP)=0 or if τP≤0, and thus have the same value in these settings. It is easy to show that

The following theorem establishes the optimality of the stochastic rule dP in a resource-limited setting.

Note that the above theorem does not claim that dP is the unique optimal stochastic regime. For discrete V, the above theorem is an immediate consequence of the discussion of the knapsack problem in Dantzig [9].

In this paper we focus on the value of the optimal stochastic rule. Nonetheless, the techniques that we present in this paper will only yield valid inference in the case where the data are generated according to a distribution P0 for which Pr0(Qˉb,0(V)=τ0)=0. This is analogous to assuming a non-exceptional law in settings where resources are not limited [13, 2], though we note that for continuous covariates V this assumption is much more likely if τ0>0. It seems unlikely that the treatment effect in some positive probability strata of covariates will concentrate on some arbitrary (determined by the constraint κ) value τ0. Nonetheless, one could deal with situations where Pr0(Qˉb,0(V)=τ0)>0 using similar martingale-based online estimation techniques to those presented in Luedtke and van der Laan [13].

3 Estimating the optimal optimal R-C value

We now present an estimation strategy for the optimal R-C rule. The upcoming sections justify this strategy and suggest that it will perform well for a wide variety of data generating distributions. The estimation strategy proceeds as follows:

1. Obtain estimates Qˉn, Qˉb,n, and gn of Qˉ0, Qˉb,0, and g0 using any desired estimation strategy which respects the fact that Y is bounded in the unit interval.

2. Estimate the marginal distributions of W and V with the corresponding empirical distributions.

3. Estimate S0 with the plug-in estimator Sn given by τ↦1n∑i=1nIQˉb,n(vi)>τ.

4. Estimate η0 with the plug-in estimator ηn=Δinfτ:Sn(τ)≤κ.

5. Estimate τ0 with the plug-in estimator given by τn=Δmax{ηn,0}.

6. Estimate d0 with the plug-in estimator dn with distribution

   PrPU(dn(U,v)=1)={κ−Sn(τn),ifQˉb,n(v)=τnandτn>0I(Qˉb,n(v)>τn),otherwise.

7. Run a TMLE for the parameter Ψdn(P0):

   1. For a˜∈{0,1}, define H(a,w)=ΔPrPU(dn(U,v)=a)gn(a|w). Run a univariate logistic regression using:

      Outcome:(yi:i=1,...,n)

      Offset:logitQˉn(ai,wi):i=1,...,n

      Covariate:H(ai,wi):i=1,...,n.

      Let εn represent the estimate of the coefficient for the covariate, i.e.

      εn=Δargmaxε∈R1n∑i=1nQˉnε(ai,wi)logyi+1−Qˉnε(ai,wi)log(1−yi),

      where Qˉnε(a,w)=Δlogit−1logitQˉn(a,w)+εH(a,w).

   2. Define Qˉn∗=ΔQˉnεn.

   3. Estimate Ψdn(P0) using the plug-in estimator given by

We use Ψdn(Pn∗) as our estimate of Ψ(P0). We will denote this estimator Ψˆ, where we have defined Ψˆ so that Ψˆ(Pn)=Ψdn(Pn∗). Note that we have used a TMLE for the data dependent parameter Ψdn(P0), which represents the value under a stochastic intervention dn. Nonetheless, we assume that PrP0(Qˉb,0(V)=τ0)=0 for many of the results pertaining to our estimator Ψˆ, i.e. we assume that the optimal R-C rule is deterministic. We view estimating the value under a stochastic rather than deterministic intervention as worthwhile because one can give conditions under which the above estimator is (root-n) consistent for Ψ(P0) at all laws P0, even if non-negligible bias invalidates standard Wald-type confidence intervals for the parameter of interest at laws P0 for which PrP0(Qˉb,0(V)=τ0)>0.

We will use Pn∗ to denote any distribution for which QˉPn∗=Qˉn∗, gPn∗=gn, and Pn∗ has the marginal empirical distribution of W for the marginal distribution of W. We note that such a distribution Pn∗ exists provided that Qˉn∗ and gn fall in the parameter spaces of P↦QˉP(W) and P↦gP, respectively.

In practice we recommend estimating Qˉ0 and Qˉb,0 using an ensemble method such as super-learning to make an optimal bias-variance trade-off (or, more generally, minimize cross-validated risk) between a mix of parametric models and data adaptive regression algorithms [18, 19]. If the treatment mechanism g0 is unknown then we recommend using similar data adaptive approaches to obtain the estimate gn. If g0 is known (as in a randomized controlled trial without missingness), then one can either take gn=g0 or estimate g0 using a correctly specified parametric model, which we expect to increase the efficiency of estimators when the Qˉ0 part of the likelihood is misspecified [20, 21].

There is typically little downside to using data adaptive approaches to estimate the needed portions of the likelihood, though we do give a formal empirical process condition in Section 5.1 which describes exactly how data adaptive these estimators can be. If one is concerned about the data adaptivity of the estimators of the needed portions of the likelihood, then one can consider a cross-validated TMLE approach such as that presented in van der Laan and Luedtke [20] or an online one-step estimator as that presented in Luedtke and van der Laan [13]. These two approaches make no restrictions on the data adaptivity of the estimators of Qˉ0, Qˉb,0, or g0.

We now outline the main results of this paper, which hold under appropriate consistency and regularity conditions.

1. Asymptotic linearity of Ψˆ:

with D0 a known function of P0.

1. Ψˆ is an asymptotically efficient estimate of Ψ(P0).

2. One can obtain a consistent estimate σn2 for the variance of D0(O). An asymptotically valid 95% confidence intervals for Ψ(P0) given by Ψˆ(Pn)±1.96σn/n.

The upcoming sections give the consistency and regularity conditions which imply the above results.

4 Canonical gradient of the optimal R-C value

The pathwise derivative of Ψ will provide a key ingredient for analyzing the asymptotic properties of our estimator. We refer the reader to Pfanzagl [22] and Bickel et al. [23] for an overview of the crucial role that the pathwise derivative plays in semiparametric efficiency theory. We remind the reader that an estimator Φˆ is an asyptotically linear estimator of a parameter Φ(P0) with influence curve ICP0 provided that

If Φ is pathwise differentiable with canonical gradient ICP0, then Φˆ is RAL and asymptotically efficient (minimum variance) among all such RAL estimators of Φ(P0) [23, 22].

For o∈O, a deterministic rule d˜, and a real number τ, define

where gP(a|W)=ΔPrP(A=a|W). We will let g0=ΔgP0. We note that D1(d˜,P)+D2(d˜,P) is the efficient influence curve of the parameter Ψ˜d˜(P).

Let d be some stochastic rule. The canonical gradient of Ψd is given by

Define

For ease of reference, let D0=ΔD(d0,τ0,P0). The upcoming theorem makes use of the following assumptions.

• C1. g0 satisfies the positivity assumption: Pr0(0<g0(1|W)<1)=1.

• C2. Qˉb,0(W) has density f0 at η0, and 0<f0(η0)<∞.

• C3. S0 is continuous in a neighborhood of η0.

• C4. Pr0(Qˉb,0(V)=τ)=0 for all τ in a neighborhood of τ0.

We now present the canonical gradient of the optimal R-C value.

Note that C3) implies that Pr0(Qˉb,0(V)=τ0)=0. Thus d0 is (almost surely) deterministic and the expectation over PU in the definition of D0 is superfluous. Nonetheless, this representation will prove useful when we seek to show that our estimator solves the empirical estimating equation defined by an estimate of D(d0,τ0,P0).

When the resource constraint is active, i.e. τ0>0, the above theorem shows that Ψ has an additional component over the optimal value parameter when no resource constraints are present [11]. The additional component is τ0×EPU[d0(U,v)]−κ, and is the portion of the derivative that relies on the fact that d0 is estimated and falls on the edge of the parameter space. We note that it is possible that the variance of D0(O) is greater than the variance of ICd0(P0)(O). If τ0=0 then these two variances are the same, so suppose τ0>0. Then, provided that Pr0(Qˉb,0(V)=τ0)=0, we have that

For any κ∈(0,1), it is possible to exhibit a distribution P0 which satisfies the conditions of Theorem 3 and for which VarP0(D0(O))>VarP0(ICd0(P0)(O)). Perhaps more surprisingly, it is also possible to exhibit a distribution P0 which satisfies the conditions of Theorem 3 and for which VarP0(D0(O))<VarP0(ICd0(P0)(O)). We omit further the discussion here because the focus of this work is on considering the estimating the value from the optimization problem (2), rather than discussing how this procedure relates to the estimation of other parameters.

5 Results about the proposed estimator

We now show that Ψˆ is an asymptotically linear estimator for Ψ(P0) with influence curve D0 provided our estimates of the needed parts of P0 satisfy consistency and regularity conditions. Our theoretical results are presented in Section 5.1, and the conditions of our main theorem are discussed in Section 5.2.

6 Simulation methods

We simulated i.i.d. draws from two data generating distributions at sample sizes 100, 200, and 1,000. For each sample size and distribution we considered resource constraints κ=0.1 and κ=0.9. We ran 2,000 Monte Carlo draws of each simulation setting. All simulations were run in R [26].

We first present the two data generating distributions considered, and then present the estimation strategies used.

7 Simulation results

Figure 1:

Coverage of two-sided 95% confidence intervals. As expected, coverage increases with sample size. The coverage tends to be better for κ=0.1 than for κ=0.9, though the estimator performed well at the largest sample size (1,000) for all simulations and choices of κ. Error bars indicate 95% confidence intervals to account for uncertainty from the finite number of Monte Carlo draws.

The proposed estimation strategy performed well overall. Figure 1 demonstrates the coverage of 95% confidence intervals for the optimal R-C value. All methods performed well at all sample sizes for the highly constrained setting where κ=0.1. The results were more mixed for the resource constraint κ=0.9. All methods performed well at the largest sample size considered. This supports our theoretical results, which were all asymptotic in nature. For Simulation 1, in which the resource constraint was not active for either choice of V, the coverage dropped off at lower sample sizes. Coverage was approximately 90% in the two smaller sample sample sizes for V=W3, which may be expected for such an asymptotic method. For the more complex problem of estimating the optimal value when V=W1,...,W4 the coverage was somewhat lower (80% when n=100 and 84% when n=200). In Simulation 2, the coverage was better (>91%) for the smaller sample sizes. We note that the resource constraint was still active (τ0>0) when κ=0.9 for this simulation, and also that the estimation problem is easier because the baseline covariate was univariate.

We report the average confidence interval widths across the 2,000 Monte Carlo draws. For n=100, average confidence interval widths were between 0.25 and 0.26 across all simulations and choices of κ. For n=200, all average confidence interval widths were between 0.17 and 0.18. For n=1000, all average confidence interval widths were approximately 0.08. We note that the usefulness of such confidence intervals varies across simulations and choices of κ. When V=W3 and κ=0.1 in Simulation 1, the optimal R-C value is approximately 0.493, versus a baseline value μ0=EP0[Qˉ0(0,W)] of approximately 0.464. Thus here the confidence interval would give the investigator little information, even at a sample size of 1,000. In Simulation 2 with κ=0.9, on the other hand, the optimal R-C value is approximately 0.572, versus a baseline value of μ0≈0.3. Thus here all confidence intervals would likely be informative for investigators, even those made for data sets of size 100.

Figure 2:

Power of the α=0.025 level test of H0:Ψ(P0)=μ0 against H1:Ψ(P0)>μ0, where μ0=EP0[Qˉ0(0,W)] is treated as known. Power increases with sample size and κ. Error bars indicate 95% confidence intervals to account for uncertainty from the finite number of Monte Carlo draws.

Figure 2 gives the power of the α=0.025 level test H0:Ψ(P0)=μ0 against the alternative H1:Ψ(P0)>μ0. Overall our method appears to have reasonable power in this statistical test. We see that power increases with sample size, the key property of consistent statistical tests. We also see that power increases with κ, which is unsurprising given that Y is binary and g0(a|w) is 1/2 for all a,w. We note that power will not always increase with κ, for example if P0 is such that g0(1|w) is very small for individuals with covariate w who are treated at κ=0.9 but not at κ=0.1. This observation is not meant as a criticism to the estimation scheme that we have presented because we assume that κ will be chosen to reflect real resource constraints, rather than to maximize the power for a test H0′:Ψ(P0)=μ′ versus H1′:Ψ(P0)>μ′ for some fixed μ′.