Nonparametric estimation of a covariate-adjusted counterfactual treatment regimen response curve

Let $T$ be a univariate failure time, $\bm{W}\in\mathcal{W}\subset\mathbb{R}^{p}$ be a vector of all available baseline covariates measured before treatment $A\in\{0,1\}$, and $C$ be a censoring variable. Let $\bm{S}\in\mathcal{S}\subset\mathbb{R}^{q}$ be a subvector of $\bm{W}$ and $\bm{V}=\bm{W}\setminus\bm{S}$ where $\bm{V}\in\mathcal{V}\subset\mathbb{R}^{p-q}$. For simplicity of notation, we denote $p-q$ as $r$. A decision rule $d^{\theta}\in\mathcal{D}$ is defined as a function that maps values of $\bm{S}$ to a treatment option $\{0,1\}$. In our formulation $\mathcal{D}$ is a class of deterministic rules indexed by a vector of coefficients. Specifically, we consider $d^{\theta}=I\{\bm{S}^{\top}\theta>0\}$ with $\theta\in\Theta\subset\mathbb{R}^{q}$ where $I(.)$ is an indicator function. Let $T^{\theta}$ be the potential outcome that would have been observed under the decision rule $d^{\theta}$ for a given value of $\theta\in\Theta$. Suppose that we have the full data $X=\{(T^{\theta}:d^{\theta}\in\mathcal{D}),\bm{W}\}\sim P_{X}\in\mathcal{M}^{F}$ where $\mathcal{M}^{F}$ is a nonparametric full data model. Define $E_{P_{X}}\{I(T^{\theta}>t)\mid\bm{V}\}=\Psi_{P_{X}}(\theta,\bm{V})\in\mathcal{F}_{\lambda}$ where $\mathcal{F}_{\lambda}$ is a class of cadlag functions with sectional variation norm bounded by $\lambda$ and $t$ is a predetermined time point. The full-data risk of $\Psi$ is defined as $E_{P_{X}}\int_{\theta}\ell_{\theta}(\Psi,x)dF(\theta)$ where $\ell$ is the log-likelihood loss function, $\ell_{\theta}(\Psi)=I(T^{\theta}>t)\log\left(\frac{\Psi}{1-\Psi}\right)+\log(1-\Psi),$ and $dF(\theta)$ is a weight function/measure on the Euclidean set of all $\theta$-values (e.g., multivariate normal with possibly large variance). Define a V-specific optimal decision rule that maximizes the $t$-year survival probability as $\theta_{0}(\bm{V})=\operatorname*{argmax}_{\theta\in\Theta}\Psi_{P_{X}}(\theta,\bm{V})$. Hence, we allow $\theta$ to be a functional parameter mapping $\mathcal{V}$ onto $\mathbb{R}^{q}$.

The observed data $\mathcal{O}=\{Y=\min(T,C),\Delta=I(T<C),A,\bm{W}\}$ follow some probability distribution $P_{0}$ that lies in some nonparametric model $\mathcal{M}$. Suppose we observe $n$ independent, identically distributed trajectories of $\mathcal{O}$. Let $\Delta^{c}=\Delta+(1-\Delta)I\{Y>t\}$ and $G:\mathcal{M}\rightarrow\mathcal{G}$ be a functional nuisance parameter where $\mathcal{G}=\{G(P):P\in\mathcal{M}\}$. In our setting, $G(P)\{(A,\Delta^{c})\mid\bm{W}\}$ denotes the joint treatment and censoring mechanism under a distribution $P$. We refer to $G(P_{0})$ as $G_{0}$ and $G(P)$ as $G$. Under coarsening at random assumption, $G(P)\{(A,\Delta^{c})\mid{X}\}=G(P)\{(A,\Delta^{c})\mid\bm{W}\}$. Denote $G^{a}\equiv G^{a}(P)(A\mid\bm{W})$ and $G^{c}\equiv G^{c}(P)(\Delta^{c}\mid\bm{W},A)$. We denote the latter two functions under $P_{0}$ as $G^{a}_{0}$ and $G^{c}_{0}$. Also define $Q:\mathcal{M}\rightarrow\mathcal{Q}$ and $\mathcal{Q}=\{Q(P):P\in\mathcal{M}\}$ where $Q_{0}(\theta,\bm{W})=\mathbb{E}_{P_{0}}\{I(Y^{\theta}>t)\mid\bm{W}\}$. Let $\prod$ be a projection operator in the Hilbert space $L_{0}^{2}(P)$ with inner product $\langle h_{1},h_{2}\rangle=E_{P}(h_{1}h_{2})$.

The loss function $\int_{\theta}\ell_{\theta}(\Psi)dF(\theta)$ is based on the full-data which is not observable, and thus, cannot be used to define our full data functional parameter $\Psi_{P_{X}}$. However, under consistency ($T_{i}=T^{A_{i}}$) and no unmeasured confounder ($A\perp T^{a}\mid\bm{W}$ and $\Delta^{c}\perp T\mid A,\bm{W}$) assumptions, it can be identified using the observed data distribution $P_{0}$. Specifically, we define an inverse probability weighting mapping of the full-data loss as

Accordingly, define the regimen-response curve as $\Psi_{P_{0}}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda}}E_{P_{0}}L_{G_{0}}(\Psi)$. Following our convention, we denote $\Psi_{P_{0}}$ as $\Psi_{0}$. We define a highly adaptive estimator of $\Psi_{P_{0}}$ as $\Psi_{n}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda_{n}}}P_{n}L_{G_{n}}(\Psi)$ where $\lambda_{n}$ is a data adaptively selected $\lambda$ and $G_{n}$ is an estimator of $G_{0}$. The choice of function $\xi$ is inconsequential to the derivation of the efficient influence function of our target parameter but it may impact the level of undersmoothing needed in finite sample. One may consider $\xi(\theta,\bm{V})=1$ or $\xi(\theta,\bm{V})=G\{(A,\Delta^{c})\mid\bm{V}\}$. The latter choice corresponds to the stabilized weight.

The highly adaptive lasso is a nonparametric regression approach that can be used to estimate infinite dimensional functional parameters Benkeser and Van Der Laan (2016); van der Laan (2017); van der Laan and Bibaut (2017). It forms a linear combination of indicator basis functions to minimize the expected value of a loss function under the constraint that the constraint that the $L_{1}$-norm of the coefficient vector is bounded by a finite constant value.

Let $\mathcal{D}[0,\tau]$ be the Banach space of d-variate real valued cadlag functions (right-continuous with left-hand limits) on a cube $[0,\tau]\in\mathbb{R}^{d}$. For each function $f\in\mathcal{D}[0,\tau]$ define the supremum norm as $\|f\|_{\infty}=\sup_{w\in[0,\tau]}|f(w)|$. For any subset $s$ of $\{0,1,\ldots,d\}$, we partition $[0,\tau]$ in $\{0\}\{\cup_{s}(0_{s},\tau_{s}]\}$ and define the sectional variation norm of such $f$ as

where the sum is over all subsets of $\{0,1,\ldots,d\}$. For a given subset $s\subset\{0,1,\ldots,d\}$, define $u_{s}=(u_{j}:j\in s)$ and $u_{-s}$ as the complement of $u_{s}$. Then, $f_{s}:[0_{s},\tau_{s}]\rightarrow\mathbb{R}$ defined as $f_{s}(u_{s})=f(u_{s},0_{-s})$. Thus, $f_{s}(u_{s})$ is a section of $f$ that sets the components in the complement of $s$ equal to zero and varies only along the variables in $u_{r}$.

Let $G^{a}\equiv G(P)(A\mid W)$ and $G^{c}\equiv G(P)(\Delta^{c}\mid W,A)$ denote the treatment and censoring mechanism under an arbitrary distribution $P\in\mathcal{M}$. Assuming that our nuisance functional parameter $G^{a},G^{c}\in\mathcal{D}[0,\tau]$ has finite sectional variation norm, we can represent $\text{logit }G^{a}$ as (Gill, van der Laan and Wellner, 1993)

The representation (2) can be approximated using a discrete measure that puts mass on each observed $W_{s,i}$ denoted by $\alpha_{s,i}$. Let $\phi_{s,i}(c_{s})=I(w_{i,s}\leq c_{s})$ where $w_{i,s}$ are support points of $dG^{a}_{s}$. We then have

where $|\alpha_{0}|+\sum_{s\subset\{1,\ldots,d\}}\sum_{i=1}^{n}|\alpha_{s,i}|$ is an approximation of the sectional variation norm of $f$. The loss based highly adaptive lasso estimator $\beta_{n}$ is defined as

where $L(.)$ is a given loss function. Denote $G^{a}_{n,\lambda}\equiv G^{a}_{\alpha_{n}(\lambda)}$ as the highly adaptive lasso estimate of $G^{a}$. Similarly, we can define $G^{c}_{n,\lambda}$ as the highly adaptive lasso estimate $G^{c}$. When the nuisance parameter correspond to a conditional probability of a binary variable (e.g., propensity score with a binary treatment or a binary censoring indicator), log-likelihood loss can be used. In practice, $\lambda$ is unknown and must be obtained using data. We refer to $\lambda_{n}$ as a value of $\lambda$ that is determined using data.

The target functional parameter $\Psi_{0}$ is not a pathwise differentiable function which makes the inference challenging. The existing literature approximates this functional by pathwise differentiable function $\Psi_{\beta}$, where $\beta$ is a vector of unknown parameters van der Laan (2006); van der Laan and Petersen (2007); Murphy et al. (2001); Robins, Orellana and Rotnitzky (2008); Orellana, Rotnitzky and Robins (2010); Neugebauer and van der Laan (2007). The working model provides a user-specified summary of the unknown true function relating the expectation of the potential outcome under a strategy. Hence, the quality of the constructed decision rule depends on how well the working model approximates the true function. To mitigate this concern we propose a kernel smoothed highly adaptive lasso estimator of $\Psi_{0}(\bm{v}_{0},\theta)$. Specifically, we define an regimen-response curve estimator as

where $U(P_{n},\bm{v}_{0},\Psi_{n})=(P_{n}K_{h,\bm{v}_{0}})^{-1}\{K_{h,\bm{v}_{0}}\Psi_{n}(\theta)\}$ and $\Psi_{n}$ is a highly adaptive lasso estimate of $\Psi_{0}$ obtained based on the loss function (1). In this formulation, $K_{h,\bm{v}_{0}}(v)=h^{-r}K\{(v-\bm{v}_{0})h^{-1}\}$ with $K(v)=K(v_{1},\cdots,v_{r})$ where $K(v)$ is a kernel satisfying $\int K(v)dv=1$. Moreover, we assume that $K(v)=\prod_{j=1}^{r}K_{u}(v_{j})$ is a product kernel defined by a univariate kernel $K_{u}(.)$ and $K(v)$ is a J-orthogonal kernel such that $\int K_{u}(v)dv=1$ and $\int K_{u}(v)v^{j}dv=0$ for $j=1,\cdots,J$.

We study the asymptotic behaviour of our estimator relative to $\Psi_{0}(\bm{v_{0}},\theta)$ and a kernel smoothed parameter

where $h$ is a fixed bandwidth. In general, when the nuisance parameter $G_{0}$ is estimated using data adaptive regression techniques, the resulting regimen-response curve estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ won’t be asymptotically linear which makes inference challenging. For example, consider the target parameter $\Psi_{0h}(\bm{v}_{0},\theta)$, we have

The first term on the right-hand side of (4) is exactly linear and the third term is negligible under Donsker condition of $\Psi_{n}$. The problem arises because data adaptive regression techniques have a rate of convergence slower than the desired root-$n$ rate, and thus, the bias of $P_{0}\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\}$ diverges to infinity.

We show that when the nuisance parameters are properly undersmoothed, the asymptotic linearity of our estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ can be retrieved. Specifically, in the proof of Theorem 4.2, we show that

where $\tilde{D}^{*}_{h,\bm{v_{0}}}(P_{0})$ includes certain components of the canonical gradient of our statistical parameter $\Psi_{0h}(\bm{v_{0}},\theta)$ and thus, contributes to the efficient influence function of our estimator. The last three terms do not converge to zero at the appropriate rate and thus induce non-negligible bias. However, the terms (5), (6) and (7) resemble the form of score equations for $I(Y\geq t)$, $\Delta^{c}$ and $I(A=d^{\theta})$.

We first discuss the score equations for $G^{.}$ where $.\in\{a,c\}$. We can generate score $S_{g}(G_{n}^{.})$ for a nonparametric estimator $G_{n}^{.}$ by paths $\{G_{n,\epsilon}^{.g}(w)\}$, such that

where $g$ is a uniformly bounded function on $[0,\tau]^{p}$. Accordingly when a function is estimated using a highly adaptive lasso, the score functions can be generated by a path $\{1+\epsilon g(s,j)\}\beta_{n,s,j}$ for a bounded vector $g$ as

where $L(\cdot)$ is the log-likelihood loss function. For example, for the treatment indicator, $L(G^{a})=A\log\left(\frac{G^{a}}{1-G^{a}}\right)+\log(1-G^{a})$. The penalty in the highly adaptive lasso imposes a restriction on the $g$ function. Specifically, the path $\{1+\epsilon g(s,j)\}\beta_{n,s,j}$ can be generated using those $g$ functions that do not change the $L_{1}$-norm of the parameters. Let $\beta^{g}_{n,s,j}=\{1+\epsilon g(s,j)\}\beta_{n,s,j}$ denote the set of perturbed parameters. Then, for small enough $\epsilon$ such that $\{1+\epsilon g(s,j)\}>0$, $\|\beta^{g}_{n,s,j}\|_{1}=\|\beta^{g}_{n,s,j}\|_{1}\{1+\epsilon g(s,j)\}=\|\beta^{g}_{n,s,j}\|_{1}+\epsilon g(s,j)\|\beta_{n,s,j}\|_{1}$. Hence, the restriction is satisfied when the inner product of $g$ and the vector $|\beta|$ is zero (i.e., $<g.|\beta|>=0$).

We now provide a thought experiment on how undersmoothing the nuisance function estimates using a highly adaptive lasso can eliminate the terms (5), (6) and (7). We first ignore the $L_{1}$-norm restriction on the choice of function $g$ which is erroneous but contains the main idea of the approach. Without the restriction, one can choose $g(s_{0},j_{0})=I(s=s_{0},j=j_{0})$. The latter choice perturbs one coefficient at a time and corresponds to maximum likelihood estimators. Then, for any $s_{0}$ and $j_{0}$, the score function (19) becomes

Because, $\beta_{n,s_{0},j_{0}}$ is a finite constant for all $s_{0}$ and $j_{0}$, solving the above score equation is equivalent to solving

For example, for the treatment indicator, the score equation is given by $(A-G_{n,\lambda_{n}}^{a})\phi_{s_{0},j_{0}}$. As we undersmooth the fit, we solve more and more score equations (i.e., the number of score equations increases with the number of features included in the model) and any linear combination of those score equations will also be solved. Now let’s consider the term (7) and let $f=\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{G_{n}^{a}P_{n}K_{h,\bm{v}_{0}}}$. Assuming that $Q_{0}$, $\Psi_{0}$ and $G_{0}=G_{0}^{a}G_{0}^{c}$ are càdlàg with finite sectional variation norm (Assumption 1 in Section 4), Gill, van der Laan and Wellner (1993) showed that $f$ can be approximated as a linear combination of indicator basis functions. Therefore, if we sufficiently undersmooth $G_{n,\lambda}^{a}$ we will solve (7) up to a root-$n$ factor. That is

The same argument can be made to show that the other terms (5) and (6) can be made negligible when $\Psi_{n}$ and $G_{n}^{c}$ are properly undersmoothed. Hence, the challenging term will be asymptotically linear with

Now, we study the actual case where the choice function $h$ is restricted to those that do not change the $L_{1}-$norm of the coefficients. Lemma A.2 shows than when (23), (24) and (25) are satisfied, our proposed estimator achieves asymptotic linearity. As we undersmooth the fit, we start adding features with small coefficients into the model. Conditions (23), (24) and (25) imply that the $L_{1}-$norm must be increases until one of the corresponding score equations is solved to a precision of $o_{p}(n^{-1/2})$.

Here we provide details for the score equations corresponding to the treatment indicator where $S_{g}(G_{n,\lambda_{n}}^{a})=(A-G_{n,\lambda_{n}}^{a})\left\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\}$. Let $r(g,G_{n,\lambda_{n}}^{a})=\sum_{(s,j)}g(s,j)\lvert\beta_{n,s,j}\rvert$. For small enough $\epsilon$,

Hence, for any $g$ satisfying $r(g,G_{n,\lambda_{n}}^{a})=0$ (i.e., $h$ does not change the $L_{1}-$ norm of the coefficients), we have $P_{n}S_{g}(G_{n,\lambda_{n}}^{a})=0$.

Let $D\{f,G_{n,\lambda_{n}}^{a}\}=f\cdot(A-G_{n,\lambda_{n}}^{a})$, where $f$ is defined above, and let $\tilde{f}$ be an approximation of $f$ using the basis functions that satisfy condition (24). Then, $D\{\tilde{f},G_{n,\lambda_{n}}^{a}\}\in\{S_{g}(G_{n,\lambda_{n}}^{a}):\lVert g\rVert_{\infty}<\infty\}$. Thus, there exists $g^{\star}$ such that $D(\tilde{f},G_{n,\lambda_{n}}^{a})=S_{g^{\star}}(G_{n,\lambda_{n}}^{a})$; however, for this particular choice of $g^{\star}$, $r(g,G_{n,\lambda_{n}}^{a})$ may not be zero (i.e., the restriction on $g$ might be violated). Now, define $g$ such that $g(s,j)=g^{\star}(s,j)$ for $(s,j)\neq(s^{\star},j^{\star})$; $\tilde{g}(s^{\star},j^{\star})$ is defined such that

That is, $g$ matches $g^{\star}$ everywhere but for a single point $(s^{\star},j^{\star})$, where it is forced to take a value such that $r(g,G_{n,\lambda_{n}}^{a})=0$. As a result, for such a choice of $g$, $P_{n}S_{h}(G_{n,\lambda_{n}}^{a})=0$ by definition. Below, we show that $P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D(\tilde{f},G_{n,\lambda_{n}}^{a})=o_{p}(n^{-1/2})$ which then implies that $P_{n}D(\tilde{f},G_{n,\lambda_{n}}^{a})=o_{p}(n^{-1/2})$. We note that the choice of $(s^{\star},j^{\star})$ is inconsequential.

The details of the forth equality is given in the proof of Lemma A.2 and the last equality follows from the assumption that $\min_{(s,j)\in\mathcal{J}_{n}}\lVert P_{n}\frac{d}{d\text{logit }G^{a}_{n,\lambda_{n}}}L(\text{logit }G^{a}_{n,\lambda_{n}})(\phi_{s,j})\rVert=o_{p}(n^{-1/2})$ for $L(\cdot)$ being log-likelihood loss (i.e., condition (24)). As $P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$, it follows that $P_{n}D(\tilde{f},G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$. Using this result, under mild assumptions, we showed that $P_{n}D(f,G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$ indicating that the term (7) will be asymptotically negligible (see the proof of Lemma A.2 for details).

To improve the finite sample performance of our estimator we propose to estimate the nuisance parameters using cross-fitting  (Klaassen, 1987; Zheng and van der Laan, 2011; Chernozhukov et al., 2017). We split the data at random into $B$ mutually exclusive and exhaustive sets of size approximately $nB^{-1}$. Let $P_{n,b}^{0}$ and $P_{n,b}^{1}$ denote the the empirical distribution of a training and a validation sample, respectively. For a given $\lambda$ and $h$, exclude a single (validation) fold of data and fit the highly adaptive lasso estimator using data from the remaining $(B-1)$ folds; use this model to estimate the nuisance parameters for samples in the holdout (validation) fold. By repeating the process $B$ times, we will have estimates of the nuisance parameters for all sample units. Accordingly, we define the cross-fitted IPW estimator

where $U(P_{n,b}^{1},\bm{v}_{0},\Psi_{n,\lambda,b})=(G_{n,\lambda,b}P_{n,b}^{1}K_{h,\bm{v}_{0}})^{-1}\{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})\Psi_{n,\lambda,b}\}$, $\Psi_{n,\lambda,b}$ and $G_{n,\lambda,b}$ are highly adaptive lasso estimates of $\Psi_{0}$ and $G_{0}$ applied to the training sample for the b^th sample split for a given $\lambda$. The estimator uses a ($J-1$)-orthogonal kernel with bandwidth $h$ centered at $\bm{v_{0}}$.

The proposed estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ approaches to $\Psi_{0}(\bm{v}_{0},\theta)$ as $h_{n}\rightarrow 0$ at the cost of increasing the variance of the estimator. Mean square error (MSE) provides an ideal criterion for bias-variance trade-off. However, because the bias is unknown the latter cannot be used directly. Here we propose an approach that circumvent the need for knowing bias while targeting the optimal bias-variance trade-off. Let $\sigma_{nh}$ be a consistent estimator of $\sigma_{0}=[E\{D^{*2}_{\bm{v}_{0}}(P_{0})\}]^{1/2}$ where $D^{*2}_{\bm{v}_{0}}(P_{0})$ is the efficient influence function for $\Psi_{nh}(\bm{v}_{0},\theta)$ (see Theorem 4.2). The goal is to find an $h$ that minimizes the MSE or equivalently set the derivative of the MSE to zero. That is

We know that the optimal (up to a constant) bias-variance trade-off implies $\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\}\approx\sigma_{nh}/(nh^{r})^{1/2}$ (van der Laan, Bibaut and Luedtke, 2018). Hence, the optimal $h$ will also solve one of the following equations

where $\kappa$ is a positive constant. Consider a scenario where $\Psi_{nh}(\bm{v}_{0},\theta)$ shows an increasing trend as $h_{n}$ approaches zero. This implies that $\Psi_{0}(\bm{v}_{0},\theta)>\Psi_{nh}(\bm{v}_{0},\theta)-\kappa\sigma_{nh}/(nh^{r})^{1/2}$. We then define our optimal finite sample bandwidth selector as

Similarly, when $\Psi_{nh}(\bm{v}_{0},\theta)$ shows decreasing trend as $h_{n}$ approaches zero, we define $h_{n}=\operatorname*{argmin}_{h}\Psi_{nh}(\bm{v}_{0},\theta)+\kappa\sigma_{nh}/(nh^{r})^{1/2}$. A natural choice for the constant $\kappa$ is $(1-\alpha)$-quantile of the standard normal distribution (i.e., $\zeta_{1-\alpha}$). Hence, $h_{n}$ is the bandwidth value that minimizes the difference between the true value $\Psi_{0}(\bm{v}_{0},\theta)$ and the corresponding lower (upper) confidence bound.

The asymptotic linearity results of Theorems 4.2 and 4.3 rely on estimating the corresponding nuisance parameters using an undersmoothed highly adaptive lasso. Specifically, Theorem 4.2 requires that

• (a)

  $\left|P_{n}D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

• (b)

  $\left|P_{n}D_{CAR}^{c}(P_{n},G_{n},Q_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

• (c)

  $\left|P_{n}D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

where $D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)=\frac{K_{h,\bm{v}_{0}}}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{G_{n}^{a}-I(A=d^{\theta})}{G_{n}}\right)$, $D_{CAR}^{c}(P_{n},G_{n},Q_{n},\Psi_{n},\theta)=\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{G_{n}^{c}-\Delta^{c}}{G_{n}}\right)$ and $D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n},\theta)=\frac{K_{h,v}\Delta^{c}I(A=d^{\theta})}{G_{n}P_{n}K_{h,v}}\{\Psi_{n}-I(Y>t)\}$. Motivated by our theoretical results, we propose the following practical L1-norm bound selection criteria for $G_{n}^{a}$, $G_{n}^{c}$ and $\Psi_{n}$: