```{r set-options, echo=FALSE, cache=FALSE} options(width = 1000) knitr::opts_chunk$set(fig.width=5, fig.height=5, out.width = "40%", fig.align = 'center') knitr::opts_chunk$set(class.source = "fold-show") knitr::opts_chunk$set(collapse=TRUE) ``` ## Overview We have now seen several examples of estimating the optimal treatment regime. In this section, our goal is to explore several issues and topics around this problem. Hopefully they can provide some guidelines for you to think about when you are working on your own project. ## Outcome Weighted Learning Outcome Weighted Learning [@zhao2012estimating] is probably one of the most popular direct learning approaches. Formally, let's be precise about our target of estimation. Following [@qian2011performance], let us define a quantity called the __value function__ of a decision rule \( \pi(x) \). In fact, this is also a concept in reinforcement learning, and it is more common to use \( R \) (reward) to denote the outcome rather than \(Y \). Hence we will follow the notation in the literature. Let \(R(a)\) denote the potential reward under treatment \(a\). Under consistency, conditional exchangeability \(R(a) \perp A \mid X\), and positivity, the value function can be identified as follows: \[ \begin{aligned} {\cV}(\pi) =& \E^\pi(R)\\ =& \E\left[R\{\pi(X)\}\right]\\ =& \E_X\big[ \E_R\{R \mid X, A = \pi(X)\} \big]. \end{aligned} \] This means that we first (inner expectation \( E_R \)) let everyone take the treatment label suggested by the decision rule \( \pi \), i.e., \( A = \pi(X) \), and obtain the expected reward, and then average over the entire population (\(E_X\)). Hence, our goal is to maximize ${\cV}(\pi)$ by searching for the best \(\pi(X)\): \[ \pi_\text{opt} = \underset{\pi}{\arg\max} \,\, {\cV}(\pi) \] Note that this is the same as our previous definition of the optimal treatment regime, since if a rule maximizes the outcome for each individual, it would also maximize the population average. However, a tricky problem here is that we do not observe the distribution under which everyone follows the suggested $\pi(X)$ treatment. Instead, our observational study is collected under its own mechanism, or distribution, denoted as $(R, X, A) \sim \mu$. This is called a behavior policy. On the other hand, our target is to estimate the value function under a restricted policy that gives $(R, X, A = \pi(X)) \sim \mu^\pi$. Let $g(a \mid x) = \Pr(A=a \mid X=x)$ denote the behavior propensity. A tool we could utilize is the [Radon-Nikodym theorem](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem)[^rn], which suggests that \[ \begin{aligned} \cV(\pi) =& \int R d\mu^\pi \\ =& \int R \frac{d\mu^\pi}{d\mu}d\mu \\ =& \int \frac{R \cdot \mathbf{1}\{A = \pi(X)\}}{g(A \mid X)} d\mu. \end{aligned} \] Now, since $d\mu$ is observed, we can use the sample (empirical) estimation of this value function. \[ \widehat{\pi}_\text{opt} = \underset{\pi}{\arg\max} \,\, \frac{1}{n} \sum_{i=1}^n \frac{R_i \cdot \mathbf{1}\{A_i = \pi(X_i)\} }{\widehat g(A_i \mid X_i)}. \] We can simply change the equal sign to unequal sign and this would become a minimization problem: \[ \widehat{\pi}_\text{opt} = \underset{\pi}{\arg\min} \frac{1}{n} \sum_{i=1}^n \frac{R_i \cdot \mathbf{1}\{A_i \neq \pi(X_i)\} }{\widehat g(A_i \mid X_i)}. \] The interesting view of this formulation is that, this becomes a classification problem, where * $\pi(X_i)$ is the decision rule * $A_i$ is the observed class label * $\frac{R_i}{\widehat g(A_i \mid X_i)}$ is the subject-specific weight^[In some cases, the observed outcome can be negative. At the population level, the optimal treatment rule is invariant to a location shift of \(R\). With the correct behavior propensity, we can verify that \( \E\left[ \frac{\alpha \cdot \mathbf{1}\{A \neq \pi(X)\}}{g(A \mid X)}\right] = \alpha \) for any constant \( \alpha \). Hence, a common implementation adds a constant to all \( R_i \) to make the classification weights positive. In a finite sample, especially with estimated propensities and regularization, the fitted rule can still be affected by this shift.]. Hence, our goal is to fit a weighted classification model. We want to encourage the model to do a better job for those with high rewards (large weights). For the binary-treatment implementation below, we encode the treatment as \(A \in \{-1,+1\}\). The following code visually demonstrates how this estimator works. ```{r fig.dim=c(5, 5), out.width='40%'} set.seed(1) n = 800 p = 2 x1 <- runif(n, 0, 1) x2 <- runif(n, 0, 1) a = rbinom(n, 1, 0.5)*2-1 side = sign(x2 - sin(2*pi*x1)/3-0.5) R <- rnorm(n, mean = ifelse(side==a, 1.5, 0.5), sd = 0.5) # the observed data, with class labels par(mar=rep(1,4)) plot(x1, x2, pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = ifelse(a==1, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) ``` Now let us add the reward as weights, indicated by the point size. It is clear that on the lower side, treatment $-1$ is more dominant, while on the upper side, treatment $1$ is more dominant. The true decision boundary is also given below. ```{r fig.dim=c(10, 5), out.width='80%'} par(mfrow = c(1, 2)) par(mar=rep(1,4)) plot(x1, x2, pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", cex = (R+2)/3, col = ifelse(a==1, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) plot(x1, x2, pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", cex = (R+2)/3, col = ifelse(a==1, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) x1_grid = sort(x1) lines(x1_grid, sin(2*pi*x1_grid)/3+0.5, lwd = 3) ``` ### Choosing the Functional Class The outcome weighted learning framework permits a wide range of functional classes for estimating \(\pi\). For example, we can use a simple linear model, a function of reproducing kernel Hilbert space, a decision tree, a random forest, or a neural network. But also be aware that the functional class cannot be too flexible, otherwise it would exactly fit the observed data with \(\widehat\pi(X_i) = A_i\). The following code demonstrates the idea using a simple linear model. This can fit reasonably well but we know that a linear function would lead to a large bias for our problem because the true optimal treatment regime is a sign function. ```{r warning=FALSE} mydata = data.frame("x1" = x1, "x2" = x2, "R" = R, "a" = a) # estimate Pr(A = 1 | X), although we know this is a randomized trial pscore = glm(as.factor(a) ~ x1 + x2, data = mydata, family = "binomial")$fitted.values # convert to the probability of each subject's observed treatment pscore.observed = ifelse(a == 1, pscore, 1 - pscore) # add a constant to the reward to make the classification weights positive W = (mydata$R + 2)/pscore.observed # fit a weighted classification model to estimate the OTR fit = glm(as.factor(a) ~ x1 + x2, data = mydata, family = "binomial", weights = W) # the decision rule table(fit$fitted.values > 0.5, side) mean((fit$fitted.values > 0.5) == (side == 1)) ``` Since the linear model does not perform well, we now use a more flexible functional class based on the **Reproducing Kernel Hilbert Space (RKHS)**. Let \( f \in \mathcal{H} \) be a real-valued function, and define the decision rule as \(\pi(x) = \mathrm{sign}\{f(x)\}\). A common approach is to use the weighted support vector machine (SVM) framework, which replaces the nonconvex misclassification loss \(1\{A_i \neq \pi(X_i)\}\) with the convex hinge loss (see [this note](https://teazrq.github.io/SMLR/support-vector-machines.html#svm-as-a-penalized-model)). The optimization problem becomes \[ \widehat{f} = \arg\min_{f \in \mathcal{H}} \frac{1}{n} \sum_{i=1}^n \underbrace{\frac{R_i^{+}}{\widehat g(A_i \mid X_i)}}_{\text{subject weight}} \, \big[1 - A_i\, f(X_i)\big]_+ + \lambda \| f \|_{\mathcal{H}}^2, \qquad \widehat{\pi}(x) = \mathrm{sign}\{\widehat{f}(x)\}, \] where \(A_i \in \{-1, +1\}\), \([u]_+ = \max(u, 0)\), and \(R_i^{+} = R_i + c\) for any constant \(c\) chosen so that \(R_i^{+} \ge 0\). The term \(\| f \|_{\mathcal{H}}^2\) controls the smoothness or complexity of the function. By the [Representer Theorem](https://teazrq.github.io/SMLR/reproducing-kernel-hilbert-space.html), the solution takes the form \[ \widehat{f}(x) = \sum_{i=1}^n \alpha_i K(x, X_i), \] where \(K(\cdot, \cdot)\) is a kernel function and \(\bK\) is the corresponding kernel matrix with entries \(K(X_i, X_j)\). The weighted SVM can be efficiently implemented using the `WeightSVM` package by specifying the subject-specific weights. The following code demonstrates this idea using a radial basis function (RBF) kernel. ```{r fig.dim=c(5, 5), out.width='40%'} library(WeightSVM) owl.fit = wsvm(a ~ x1 + x2, data = mydata, kernel = "radial", gamma = 0.15, weight = W) table(owl.fit$fitted > 0, side) mean((owl.fit$fitted > 0) == (side == 1)) par(mfrow = c(1, 1)) par(mar=rep(1,4)) plot(x1, x2, pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = ifelse(owl.fit$fitted > 0, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) ``` ### Example: `lalonde` Data Let's use the `lalonde` data again. Recall that the outcome is `re78`, the treatment `treat` is whether the subject received a job training program. For this example, we will use the `DTRlearn2` package. ```{r} library(Matching) data(lalonde) head(lalonde) library(DTRlearn2) X = data.frame("age" = scale(lalonde$age), "educ" = scale(lalonde$educ), "black" = lalonde$black, "hisp" = lalonde$hisp, "married" = lalonde$married, "nodegr" = lalonde$nodegr, "re74" = lalonde$re74/max(lalonde$re74), "re75" = lalonde$re75/max(lalonde$re75)) A = 2*(as.numeric(lalonde$treat) - 0.5) R = lalonde$re78 # fit the OWL model set.seed(1) OWL.fit = owl(H = data.matrix(X), AA = A, RR = R, n = nrow(X), K = 1, # for our problem, it is only one stage of treatment so we use K = 1 kernel = "rbf", # radial basis function kernel sigma = 0.5) # bandwidth for each covariate # the estimated value function OWL.fit$valuefun ``` The package reports an estimated value function for the learned rule. Because the same data are used to learn and evaluate the rule here, this is an in-sample value. It should not by itself be used to conclude that individualization improves outcomes relative to a __one-size-fits-all__ rule. Such a comparison requires held-out or cross-fitted policy evaluation, together with an assessment of uncertainty. ## Continuous Treatment with Stochastic Policy View In the previous sections, we have mainly focused on binary treatment. However, in many cases, the treatment can be continuous. The most common example is the dosage of a drug. In this case, we need to modify the methods we have learned to accommodate continuous treatment. There are several simple approaches. For example, we can discretize the treatment into several levels and then apply the methods we have learned. However, this would greatly reduce the information we have for each treatment level. We can also use a linear regression and add interaction terms, possibly even second-order terms of \(A\), between the treatment and the covariates. When we need to choose the best treatment, we can either solve for the optimal treatment directly or use a grid search. However, a linear model can be biased if its functional form is misspecified. Hence some dedicated methods are needed. A natural idea is to utilize kernel smoothing to modify the outcome weighted learning approach. This leads to the method proposed in @chen2016personalized. Let's recall that the outcome weighted learning approach for binary/multi-categorical treatment is to solve the following optimization problem: \[ \widehat{\pi}_\text{opt} = \underset{\pi}{\arg\max} \frac{1}{n} \sum_{i=1}^n \frac{R_i \cdot \mathbf{1}\{A_i = \pi(X_i)\} }{\widehat g(A_i \mid X_i)} \] When the treatment is continuous, we need to modify the indicator function since the probability that \( A_i = \pi(X_i) \) is zero. It might be helpful to revisit our definition of the value function and the Radon-Nikodym theorem, which converted the observed behavior-policy distribution to the target-policy distribution. A deterministic target policy with \(A=\pi(X)\) is degenerate relative to a continuous behavior density. A natural alternative is to let the target treatment follow a distribution centered at \(\pi(X)\). This is connected with the reinforcement learning literature, where we use a __stochastic target policy__^[Note that in @chen2016personalized, the authors started from a different perspective by replacing exact agreement with a small tolerance while still viewing the problem as a deterministic policy with \(A = \pi(X)\). The resulting computational framework also aligns with this stochastic-policy view.]. For example, let the stochastic target policy be Uniform on \( [\pi(X) - h, \pi(X) + h] \), where \(h\) is a bandwidth. Its conditional density is \[ q_{\pi,h}(a \mid x) = \frac{1}{2h} \cdot \mathbf{1}\left\{|a-\pi(x)| \le h\right\}. \] Then if we revisit the definition of the value function, we would have \[ \begin{aligned} \cV_h(\pi) =& \int R d\mu^{\pi,h} \\ =& \int R \frac{d\mu^{\pi,h}}{d\mu}d\mu \\ =& \int \frac{R \cdot q_{\pi,h}(A \mid X)}{g(A \mid X)} d\mu \\ =& \int \frac{R \cdot \mathbf{1}\{|A-\pi(X)| \le h\}}{2h \cdot g(A \mid X)} d\mu. \\ \end{aligned} \] This leads to the sample version of the estimator, since we observe samples from \( \mu \): \[ \widehat{\pi}_{h,\text{opt}} = \underset{\pi}{\arg\max} \frac{1}{n} \sum_{i=1}^n \frac{R_i \cdot \mathbf{1}\left\{|A_i-\pi(X_i)| \le h\right\} }{2h\,\widehat g(A_i \mid X_i)}. \] This permits a variety of choices for the kernel density. Instead of using a Uniform kernel, we can in general write \[ \widehat{\pi}_{h,\text{opt}} = \underset{\pi}{\arg\max} \frac{1}{n} \sum_{i=1}^n \frac{R_i \cdot K_h\{A_i-\pi(X_i)\} }{\widehat g(A_i \mid X_i)}. \] where \(K_h(u) = K(u/h)/h\) is a kernel density. If \(K\) is nonnegative and integrates to one, this corresponds to the stochastic target-policy density \(q_{\pi,h}(a \mid x) = K_h\{a-\pi(x)\}\). In @chen2016personalized, however, the Uniform kernel is used because the resulting empirical loss can be written as a difference of convex functions and solved using the DC algorithm [@tao1998dc]. This makes the computational task easier[^dc]. ## Random Forest for Outcome Weighted Learning Many of our previously mentioned methods utilize the Reproducing Kernel Hilbert Space (RKHS) to estimate the optimal treatment regime. However, the RKHS method can be computationally expensive when the number of observations is large, may perform poorly when the treatment rule is discontinuous, and may also suffer from the curse of dimensionality. In addition, the choice of kernel function can be difficult. Random forests provide some useful advantages because they can naturally handle high-dimensional data and adapt to model sparsity. Furthermore, a random forest decision rule can handle nonlinear and discontinuous treatment rules well. Using random forests to solve the optimal treatment regime under the outcome weighted learning framework would be relatively straightforward. We can use a random forest to estimate the propensity score (using the out-of-bag prediction) and also follow the subject-weighted classification objective function to fit each tree. The following code demonstrates an example using the `RLT` package from [GitHub](https://github.com/teazrq/RLT) to implement this idea. A tutorial is provided [here](https://teazrq.github.io/RLT/index.html). Be careful that these tasks (weighted classification) using random forests can be sensitive to tuning parameters and the construction of weights. The following results were obtained after some tuning. ```{r fig.dim=c(5, 5), out.width='40%'} # install.packages("devtools") # devtools::install_github("teazrq/RLT") library(RLT) set.seed(1) n = 800 p = 5 x = matrix(runif(n*p), n, p) ps = exp(2*x[,3]-1)/(1+exp(2*x[,3]-1)) a = rbinom(n, 1, ps)*2-1 side = sign(x[,2] - sin(2*pi*x[,1])/3-0.5) R <- rnorm(n, mean = x[, 3] + x[, 4]^2 + x[, 5] + a*(x[,2] - sin(2*pi*x[,1])/3-0.5), sd = 0.5) # Extract the OOB probability for the second treatment class get.oob.class.prob = function(fit, x, class = 2) { tree.prob = vapply(predict(fit, x)$AllPrediction, function(z) z[, class], numeric(nrow(x))) oob = fit$ObsTrack == 0 rowSums(tree.prob * oob) / rowSums(oob) } # fit the propensity score pscore.fit <- RLT(x = x[, 1:p], y = as.factor(a), model = "classification", track.obs = TRUE, ntrees = 5000, mtry = 2, nmin = 2*sqrt(n)) pscore.oob = get.oob.class.prob(pscore.fit, x[, 1:p]) # We can check how well it fits the propensity score par(mar=rep(2, 4)) plot(ps, pscore.oob, xlab = "True Propensity Score", ylab = "Estimated Propensity Score", xlim = c(0.2, 0.8), ylim = c(0.2, 0.8)) abline(0, 1, col = "red", lwd = 2) ``` ```{r} # calculate subject weight by forcing them to be positive pscore.observed = ifelse(a == 1, pscore.oob, 1 - pscore.oob) W = (R - min(R)) / pscore.observed # Fit random forest with subject weight owl.fit <- RLT(x = x[, 1:p], y = as.factor(a), model = "classification", subject.weight = W, track.obs = TRUE, alpha = 0.1, ntrees = 1000, mtry = 5, nmin = 3*sqrt(n)) owl.oob = get.oob.class.prob(owl.fit, x[, 1:p]) oobrule = owl.oob > 0.5 table(oobrule, side) mean(oobrule == (side == 1)) # plot the estimated treatment regime par(mar=rep(2, 4)) plot(x[,1], x[,2], pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = ifelse(oobrule, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) ``` ## Variance Reduction of OWL One potential issue with OWL is that it may be sensitive to the scale or location of the outcome weight \( R \). Subjects with much larger weights can dominate the learning process. To mitigate this issue, we notice that, at the population level with the correct behavior propensity, subtracting any function of \( X \) from \( R \) does not change the maximizer of the problem: \[ \begin{aligned} & \quad \,\, \E \left[ \frac{ \left\{ R - m(X) \right\} \mathbf{1}\{A = \pi(X)\} }{g(A \mid X)} \right] \\ &= \E \left[ \frac{ R \mathbf{1}\{A = \pi(X)\} }{g(A \mid X)} \right] - \E \left[ \frac{ m(X) \mathbf{1}\{A = \pi(X)\} }{g(A \mid X)} \right] \\ &= \cV(\pi) - \E_X\left[ m(X) \E_A \left\{ \frac{\mathbf{1}\{A = \pi(X)\}}{g(A \mid X)} \Biggm| X \right\} \right] \\ &= \cV(\pi) - \E_X\left[ m(X) \sum_{a \in \mathcal A} \frac{\mathbf{1}\{\pi(X)=a\}g(a \mid X)}{g(a \mid X)} \right] \\ &= \cV(\pi) - \E_X\left[ m(X) \sum_{a \in \mathcal A}\mathbf{1}\{\pi(X)=a\} \right] \\ &= \cV(\pi) - \E_X\left[ m(X) \right]. \end{aligned} \] Hence, this quantity does not depend on the choice of the decision rule \( \pi \). A well-chosen \(m(X)\) can reduce the variability of the weights, although an arbitrary or poorly estimated \(m(X)\) need not reduce variance. A convenient choice used in @zhou2017residual is \[ m(X) = \E\left[ \frac{R}{2g(A \mid X)} \Biggm| X \right] = \frac{\E[R \mid X, A = 1] + \E[R \mid X, A = 0]}{2}. \] In @zhu2017greedy, the choice is the same as the one used in causal random forests: \[ m(X) = \E[ R \mid X] = \E[ R \mid X, A = 1] g(1 \mid X) + \E[ R \mid X, A = 0] g(0 \mid X). \] Both choices are intended to reduce variability and can be especially helpful for random forest splitting. The realized reduction depends on how well \(m(X)\) is estimated. The following is an example using the same dataset. ```{r} # Fit an outcome model r.fit <- RLT(x = x[, 1:p], y = R, model = "regression", track.obs = TRUE, alpha = 0.1, ntrees = 1000, mtry = 5, nmin = 3*sqrt(n)) # calculate the subject weight after adjusting the mean r.tree.pred = predict(r.fit, x[, 1:p])$AllPrediction r.oob = rowSums(r.tree.pred * (r.fit$ObsTrack == 0)) / rowSums(r.fit$ObsTrack == 0) residual = R - r.oob # compare the variance var(R) var(residual) # define the new observed-action weight W = (residual - min(residual)) / pscore.observed # Fit random forest with subject weight owl.fit <- RLT(x = x[, 1:p], y = as.factor(a), model = "classification", subject.weight = W, track.obs = TRUE, alpha = 0.1, ntrees = 1000, mtry = 5, nmin = 3*sqrt(n)) owl.oob = get.oob.class.prob(owl.fit, x[, 1:p]) oobrule = owl.oob > 0.5 table(oobrule, side) mean(oobrule == (side == 1)) # exploratory sensitivity check only: # squaring the weights changes the OWL objective and no longer # implements the same value-based criterion owl.fit <- RLT(x = x[, 1:p], y = as.factor(a), model = "classification", subject.weight = W^2, track.obs = TRUE, alpha = 0.1, ntrees = 1000, mtry = 5, nmin = 3*sqrt(n)) owl.oob = get.oob.class.prob(owl.fit, x[, 1:p]) oobrule = owl.oob > 0.5 table(oobrule, side) mean(oobrule == (side == 1)) # plot the estimated treatment regime par(mar=rep(2, 4)) plot(x[,1], x[,2], pch = 19, xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = ifelse(oobrule, "deepskyblue", "darkorange")) legend("topright", c("Treatment = 1", "Treatment = -1"), pch = 19, col = c("deepskyblue", "darkorange"), cex = 1.3) ``` ```{r} # Compare with the grf causal forest method for estimating CATE # this also includes an internal centering step to model m(X) = E[R | X] library(grf) a01 = as.numeric(a == 1) owl.grf.fit <- causal_forest(X = x[, 1:p], Y = R, W = a01, num.trees = 1000, alpha = 0.1, min.node.size = 10, mtry = 5) pred = owl.grf.fit$predictions oobrule.grf = pred > 0 table(oobrule.grf, side) mean(oobrule.grf == (side == 1)) ``` *** [^rn]: Let \((\Omega,\mathcal{F})\) be a measurable space and let \(\nu\) and \(\mu\) be \(\sigma\)-finite measures with \(\nu \ll \mu\) (i.e., whenever \(\mu(E)=0\) then \(\nu(E)=0\)). The RN theorem guarantees the existence of a measurable function \(\frac{d\nu}{d\mu}:\Omega\to[0,\infty)\) (the RN derivative, or “density” of \(\nu\) w.r.t. \(\mu\)) such that, for every integrable \(f\), \[ \int f\,d\nu \;=\; \int f\,\frac{d\nu}{d\mu}\,d\mu. \] In our notation, \(\mu\) is the observed/behavior data-generating distribution of \((R,X,A)\) and \(\mu^\pi\) is the (counterfactual) distribution induced by the target policy \(\pi\). Under **positivity** (\(\mu^\pi \ll \mu\)), the value can be rewritten as \[ \mathcal V(\pi) \;=\; \int R\, d\mu^\pi \;=\; \int R\,\frac{d\mu^\pi}{d\mu}\, d\mu. \] When actions are discrete and we factor measures by the behavior propensity \(g(a \mid x)\), the RN derivative reduces to the importance weight: \[ \frac{d\mu^\pi}{d\mu}(R,X,A) \;=\; \frac{\Pr_\pi(A\mid X)}{g(A\mid X)} \quad\Longrightarrow\quad \mathcal V(\pi) \;=\; \mathbb E\!\left[\,R\,\frac{\Pr_\pi(A\mid X)}{g(A\mid X)}\,\right]. \] [^dc]: The Difference of Convex (DC) algorithm is a general optimization framework for problems whose objective can be expressed as the difference between two convex functions, i.e., \(f(\theta) = g(\theta) - h(\theta)\) where both \(g\) and \(h\) are convex. The key idea is to linearize the concave part \( -h(\theta) \) at the current iterate \(\theta^{(t)}\) and solve the resulting convex subproblem iteratively: \[ \theta^{(t+1)} = \arg\min_{\theta} \; g(\theta) - \langle \nabla h(\theta^{(t)}), \theta \rangle. \] This yields a sequence of convex optimizations that decreases the original nonconvex objective and seeks a local minimum. In the formulation of @chen2016personalized, the Uniform-kernel construction leads to the truncated absolute loss \[ L(u) = \min(|u|, 1), \qquad u = \frac{a - \pi(x)}{h}. \] This loss is zero at exact agreement, increases linearly within the bandwidth, and remains at one for larger deviations. Thus it measures disagreement rather than approximating the inside-window indicator directly. Using the identity \[ \min(|u|,1)=|u|-(|u|-1)_+, \] the loss can be expressed as the difference of two convex functions. Hence, the empirical risk in continuous-treatment OWL can be optimized using the DC algorithm.