Targeted Learning with Stochastic Treatment Regimes


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:


Data and Notation

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 Data

# 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


Estimating Stochastic Interventions Effects with Generalized Linear Models

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 ~ .")

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 ~ .")

Interlude: Constructing Optimal Stacked Regressions with 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 and

# 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())

Estimating Stochastic Interventions Effects with Stacked Regressions

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 =
                           Q_fit_args = list(fit_type = "sl",
                                             sl_learners = sl_learner)

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 =
                           Q_fit_args = list(fit_type = "sl",
                                             sl_learners = sl_learner)

Statistical Inference for Targeted Maximum Likelihood Estimates

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:

  1. If $D(\bar{Q}_n^{\star}, g_n)$ converges to $D(P_0)$ in $L_2(P_0)$ norm, and
  2. 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}$ whp, 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)$.

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))

Advanced Usage: User-Specified Regressions

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 <-, gn_noshift, gn_upshift,
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 <-, 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)


Try the txshift package in your browser

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

txshift documentation built on Oct. 23, 2020, 8:27 p.m.