Stochastic treatment regimes present a relatively simple manner in which to assess the effects of continuous treatments by way of parameters that examine the effects induced by the counterfactual shifting of the observed values of a treatment of interest. Here, we present an implementation of a new algorithm for computing targeted minimum loss-based estimates of treatment shift parameters defined based on a shifting function $d(A,W)$. For a technical presentation of the algorithm, the interested reader is invited to consult @diaz2018stochastic. For additional background on Targeted Learning and previous work on stochastic treatment regimes, please consider consulting @vdl2011targeted, @vdl2018targeted, and @diaz2012population.

To start, let's load the packages we'll use and set a seed for simulation:

library(data.table) library(haldensify) library(txshift) set.seed(429153)

Consider $n$ observed units $O_1, \ldots, O_n$, where each random variable $O =
(W, A, Y)$ corresponds to a single observational unit. Let $W$ denote baseline
covariates (e.g., age, sex, education level), $A$ an intervention variable of
interest (e.g., nutritional supplements), and $Y$ an outcome of interest (e.g.,
disease status). Though it need not be the case, let $A$ be continuous-valued,
i.e. $A \in \mathbb{R}$. Let $O_i \sim \mathcal{P} \in \mathcal{M}$, where
$\mathcal{M}$ is the nonparametric statistical model defined as the set of
continuous densities on $O$ with respect to some dominating measure. To
formalize the definition of stochastic interventions and their corresponding
causal effects, we introduce a nonparametric structural equation model (NPSEM),
based on @pearl2000causality, to define how the system changes under posited
interventions:
\begin{align*}\label{eqn:npsem}
W &= f_W(U_W) \ A &= f_A(W, U_A) \ Y &= f_Y(A, W, U_Y),
\end{align*}
We denote the observed data structure $O = (W, A, Y)$

Letting $A$ denote a continuous-valued treatment, we assume that the
distribution of $A$ conditional on $W = w$ has support in the interval
$(l(w), u(w))$ -- for convenience, let this support be *a.e.* That is, the
minimum natural value of treatment $A$ for an individual with covariates
$W = w$ is $l(w)$; similarly, the maximum is $u(w)$. Then, a simple stochastic
intervention, based on a shift $\delta$, may be defined
\begin{equation}\label{eqn:shift}
d(a, w) =
\begin{cases}
a - \delta & \text{if } a > l(w) + \delta \
a & \text{if } a \leq l(w) + \delta,
\end{cases}
\end{equation}
where $0 \leq \delta \leq u(w)$ is an arbitrary pre-specified value that
defines the degree to which the observed value $A$ is to be shifted, where
possible. For the purpose of using such a shift in practice, the present
software provides estimators for a shift function that assumes that the density
of treatment $A$, conditional on the covariates $W$, has support *a.e.*

# simulate simple data for tmle-shift sketch n_obs <- 1000 # number of observations n_w <- 1 # number of baseline covariates tx_mult <- 2 # multiplier for the effect of W = 1 on the treatment ## baseline covariate -- simple, binary W <- as.numeric(replicate(n_w, rbinom(n_obs, 1, 0.5))) ## create treatment based on baseline W A <- as.numeric(rnorm(n_obs, mean = tx_mult * W, sd = 1)) # create outcome as a linear function of A, W + white noise Y <- A + W + rnorm(n_obs, mean = 0, sd = 1) # shift parameter delta <- 0.5

The simplest way to compute the TML estimator for the stochastic treatment
regime is to fit each of the major constituent parts (nuisance parameters) of
the estimator with generalized linear models. To do this, one may merely call
the `txshift`

function, providing the standard inputs listed below.

tmle_hal_shift_1 <- txshift(W = W, A = A, Y = Y, delta = delta, fluctuation = "standard", g_fit_args = list(fit_type = "hal", n_bins = 5, grid_type = "equal_range", lambda_seq = exp(seq(-1, -12, length = 500))), Q_fit_args = list(fit_type = "glm", glm_formula = "Y ~ .") ) summary(tmle_hal_shift_1)

When computing any such TML estimator, we may, of course, vary the regressions
used in fitting the nuisance parameters; however, an even simpler variation is
to fit the step for the fluctuation submodels with a *weighted* method, simply
weighting each observation by an auxiliary covariate (often denoted $H_n$, and
sometimes called the "clever covariate" in the literature) rather than using
such a covariate directly in the submodel regression fit.

tmle_hal_shift_2 <- txshift(W = W, A = A, Y = Y, delta = delta, fluctuation = "weighted", g_fit_args = list(fit_type = "hal", n_bins = 5, grid_type = "equal_range", lambda_seq = exp(seq(-1, -12, length = 500))), Q_fit_args = list(fit_type = "glm", glm_formula = "Y ~ .") ) summary(tmle_hal_shift_2)

`sl3`

To easily incorporate ensemble machine learning into the estimation procedure,
we rely on the facilities provided in the `sl3`

R
package [@coyle2020sl3]. For a complete guide on
using the `sl3`

R package, consider consulting https://tlverse.org/sl3/ and
https://tlverse.org/tlverse-handbook/sl3.html.

# SL learners to be used for most fits (e.g., IPCW, outcome regression) mean_learner <- Lrnr_mean$new() glm_learner <- Lrnr_glm$new() rf_learner <- Lrnr_ranger$new() Q_lib <- Stack$new(mean_learner, glm_learner, rf_learner) sl_learner <- Lrnr_sl$new(learners = Q_lib, metalearner = Lrnr_nnls$new()) # SL learners for fitting the generalized propensity score fit hse_learner <- make_learner(Lrnr_density_semiparametric, mean_learner = glm_learner ) mvd_learner <- make_learner(Lrnr_density_semiparametric, mean_learner = rf_learner, var_learner = glm_learner ) g_lib <- Stack$new(hse_learner, mvd_learner) sl_learner_density <- Lrnr_sl$new(learners = g_lib, metalearner = Lrnr_solnp_density$new())

Using the framework provided by the `sl3`

package,
the nuisance parameters of the TML estimator may be fit with ensemble learning,
using the cross-validation framework of the Super Learner algorithm of
@vdl2007super.

tmle_sl_shift_1 <- txshift(W = W, A = A, Y = Y, delta = delta, fluctuation = "standard", g_fit_args = list(fit_type = "sl", sl_learners_density = sl_learner_density), Q_fit_args = list(fit_type = "sl", sl_learners = sl_learner) ) summary(tmle_sl_shift_1)

As before, we may vary the regression for the submodel fluctuation procedure by weighting each observation by the value of the so-called clever covariate rather than using such an auxiliary covariate directly in the regression procedure:

tmle_sl_shift_2 <- txshift(W = W, A = A, Y = Y, delta = delta, fluctuation = "weighted", g_fit_args = list(fit_type = "sl", sl_learners_density = sl_learner_density), Q_fit_args = list(fit_type = "sl", sl_learners = sl_learner) ) summary(tmle_sl_shift_2)

Recall that the asymptotic distribution of TML estimators has been studied
thoroughly:
$$\psi_n - \psi_0 = (P_n - P_0) \cdot D(\bar{Q}_n^*, g_n) + R(\hat{P}^*, P_0),$$
which, provided the following two conditions:

- If $D(\bar{Q}_n^{\star}, g_n)$ converges to $D(P_0)$ in $L_2(P_0)$ norm, and
- the size of the class of functions considered for estimation of
$\bar{Q}_n^{\star}$ and $g_n$ is bounded (technically, $\exists \mathcal{F}$
st $D(\bar{Q}_n^{\star}, g_n) \in \mathcal{F}$
, where $\mathcal{F}$ is a Donsker class), readily admits the conclusion that $\psi_n - \psi_0 = (P_n - P_0) \cdot D(P_0) + R(\hat{P}^{\star}, P_0)$.**whp**

Under the additional condition that the remainder term $R(\hat{P}^{\star}, P_0)$ decays as $o_P \left( \frac{1}{\sqrt{n}} \right),$ we have that $$\psi_n - \psi_0 = (P_n - P_0) \cdot D(P_0) + o_P \left( \frac{1}{\sqrt{n}} \right),$$ which, by a central limit theorem, establishes a Gaussian limiting distribution for the estimator:

$$\sqrt{n}(\psi_n - \psi) \to N(0, V(D(P_0))),$$ where $V(D(P_0))$ is the variance of the efficient influence curve (canonical gradient) when $\psi$ admits an asymptotically linear representation.

The above implies that $\psi_n$ is a $\sqrt{n}$-consistent estimator of $\psi$, that it is asymptotically normal (as given above), and that it is locally efficient. This allows us to build Wald-type confidence intervals in a straightforward manner:

$$\psi_n \pm z_{\alpha} \cdot \frac{\sigma_n}{\sqrt{n}},$$ where $\sigma_n^2$ is an estimator of $V(D(P_0))$. The estimator $\sigma_n^2$ may be obtained using the bootstrap or computed directly via the following

$$\sigma_n^2 = \frac{1}{n} \sum_{i = 1}^{n} D^2(\bar{Q}_n^{\star}, g_n)(O_i)$$

(ci_param <- confint(tmle_hal_shift_1))

In some special cases it may be useful for the experienced user to compute the
treatment mechanism and outcome regressions separately (i.e., outside of the
`tmle_txshift`

wrapper function), instead applying this user-facing wrapper only
to invoke the *targeting* steps involved in computing the TML estimator for the
treatment shift parameter. In such cases, the optional arguments `gn_fit_ext`

and `Qn_fit_ext`

may be utilized. We present a case of using these here:

# compute treatment mechanism (propensity score) externally ## first, produce the down-shifted treatment data gn_downshift <- dnorm(A - delta, mean = tx_mult * W, sd = 1) ## next, initialize and produce the up-shifted treatment data gn_upshift <- dnorm(A + delta, mean = tx_mult * W, sd = 1) ## now, initialize and produce the up-up-shifted (2 * delta) treatment data gn_upupshift <- dnorm(A + 2 * delta, mean = tx_mult * W, sd = 1) ## then, initialize and produce the un-shifted treatment data gn_noshift <- dnorm(A, mean = tx_mult * W, sd = 1) ## finally, put it all together into an object like what's produced internally gn_out <- as.data.table(cbind(gn_downshift, gn_noshift, gn_upshift, gn_upupshift)) setnames(gn_out, c("downshift", "noshift", "upshift", "upupshift")) # compute outcome regression externally Qn_noshift <- (W + A - min(Y)) / diff(range(Y)) Qn_upshift <- (W + A + delta - min(Y)) / diff(range(Y)) Qn_noshift[Qn_noshift < 0] <- .Machine$double.neg.eps Qn_noshift[Qn_noshift > 1] <- 1 - .Machine$double.neg.eps Qn_upshift[Qn_upshift < 0] <- .Machine$double.neg.eps Qn_upshift[Qn_upshift > 1] <- 1 - .Machine$double.neg.eps Qn_out <- as.data.table(cbind(Qn_noshift, Qn_upshift)) setnames(Qn_out, c("noshift", "upshift")) # invoke the wrapper function only to apply the targeting step tmle_shift_spec <- txshift(W = W, A = A, Y = Y, delta = delta, fluctuation = "standard", ipcw_fit_args = NULL, g_fit_args = list(fit_type = "external"), Q_fit_args = list(fit_type = "external"), gn_fit_ext = gn_out, Qn_fit_ext = Qn_out) summary(tmle_shift_spec)

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.