10 Continuous and Multivalued Treatments

Code

library(tidyverse)
library(WeightIt)
library(causaldrf)
library(lmtp)
library(SuperLearner)

So far most examples use a binary treatment. Many treatments are not binary: dose of a drug, hours in a program, dollars spent, or a treatment with several levels. In those cases “the treatment effect” is not one number. We usually want a dose-response curve,

\[ \mu(d) = \mathbb{E}[Y(d)], \quad d \in \mathcal{D}. \tag{10.1}\]

I use three approaches:

Generalized propensity score (GPS; Hirano-Imbens 2004), which models the conditional density of treatment given covariates.
Doubly-robust dose-response estimation (in the spirit of Kennedy et al. 2017), which combines a GPS with an outcome model.
Longitudinal modified treatment policies [LMTP; Dı́az et al. (2023)], which are useful when treatments are time-varying or when fixed treatment levels create positivity problems.

Related reading: The LMTP material below is condensed from the longer treatment in Longitudinal modified treatment policy (LMTP) in Topics on Econometrics and Causal Inference, which follows Nicholas Williams’s tutorials at beyondtheate.com.

10.1 The continuous-treatment setup

The aim is to recover a known dose-response curve from data in which the dose is chosen partly on covariates that also move the outcome.

We simulate $n = 2000$ units with two covariates, $X_1, X_2 \sim N(0,1)$, independent of each other. The dose is $D = 2 + 0.5X_1 + 0.5X_2 + \nu$ with $\nu \sim N(0,1)$. The dose-response curve is $\mu(d) = 1 + 2d - 0.15d^2$, and the outcome is $Y = \mu(D) + 0.3X_1 - 0.3X_2 + \varepsilon$ with $\varepsilon \sim N(0,1)$.

Code

set.seed(42)
n  <- 2000
X1 <- rnorm(n)
X2 <- rnorm(n)
# Continuous treatment, depends on X1, X2
D  <- 2 + 0.5 * X1 + 0.5 * X2 + rnorm(n, sd = 1.0)
# Dose-response: nonlinear in D
true_mu <- function(d) 1.0 + 2.0 * d - 0.15 * d^2
Y  <- true_mu(D) + 0.3 * X1 - 0.3 * X2 + rnorm(n)

df <- tibble(Y = Y, D = D, X1 = X1, X2 = X2)
cat(sprintf("Treatment range: [%.2f, %.2f]\n", min(D), max(D)))

Treatment range: [-1.82, 6.39]

Code

cat(sprintf("True μ(d=2) = %.2f, μ(d=4) = %.2f, μ(d=6) = %.2f\n",
            true_mu(2), true_mu(4), true_mu(6)))

True μ(d=2) = 4.40, μ(d=4) = 6.60, μ(d=6) = 7.60

The observed doses run from $-1.82$ to 6.39, and the true curve gives $\mu(2) = 4.40$, $\mu(4) = 6.60$ and $\mu(6) = 7.60$.

The dose-response curve is concave, with clearly diminishing returns. It is worth being precise about the shape rather than calling it hump-shaped, because the vertex sits at $d = 2/0.3 \approx 6.67$, just beyond the largest dose in the sample. Over the observed range the curve is therefore increasing throughout — its slope $2 - 0.3d$ is still positive (about 0.08) at the maximum observed dose — so the downturn is extrapolation, not something these data identify. The estimation problem is to recover the curvature that is in range when the observed dose is confounded by $X$.

10.2 Naive regression vs adjusted regression

We fit three models and compare their implied dose-response curves against the truth: $Y$ on $D$ alone, $Y$ on $D$ and the covariates, and $Y$ on a quadratic in $D$ and the covariates. All three are evaluated at $X_1 = X_2 = 0$ over $d \in [0,10]$.

One feature of this DGP is worth stating before reading the plot. $X_1$ and $X_2$ are confounders in the graph, since each raises the dose and each moves the outcome. But their linear bias cancels: each raises $D$ by 0.5, while they enter $Y$ with $+0.3$ and $-0.3$, so the omitted-variable term $\text{Cov}(D,\, 0.3X_1 - 0.3X_2)/\text{Var}(D)$ is $0.3(0.5) - 0.3(0.5) = 0$ exactly. The naive slope is 1.421 and the adjusted slope 1.409. Adjusting for $X$ changes almost nothing here, and the one problem left for the linear models is functional form.

Code

# Linear fit, ignoring covariates
linear_naive <- lm(Y ~ D, data = df)
# Linear fit with covariate adjustment
linear_adj   <- lm(Y ~ D + X1 + X2, data = df)
# Nonlinear fit with covariate adjustment
poly_adj     <- lm(Y ~ poly(D, 2) + X1 + X2, data = df)

ds <- seq(0, 10, by = 0.5)

predict_mu <- function(fit, ds) {
  data.frame(d = ds, X1 = 0, X2 = 0) |>
    rename(D = d) -> nd
  predict(fit, newdata = nd)
}

# True dose-response
df_true <- tibble(d = ds, mu = true_mu(ds))

df_fits <- bind_rows(
  tibble(d = ds, mu = predict_mu(linear_naive, ds), method = "Linear (naive)"),
  tibble(d = ds, mu = predict_mu(linear_adj,   ds), method = "Linear (adj X)"),
  tibble(d = ds, mu = predict_mu(poly_adj,     ds), method = "Quadratic (adj X)")
)

Code

ggplot() +
  geom_line(data = df_fits, aes(d, mu, colour = method), linewidth = 1.0) +
  geom_line(data = df_true, aes(d, mu), colour = "black",
            linetype = "dashed", linewidth = 1.0) +
  labs(x = "Dose D", y = "Expected outcome E[Y(d)]",
       colour = "Method",
       caption = "Black dashed line = true μ(d) = 1 + 2d - 0.15 d²") +
  theme_minimal()

Fitted vs true dose-response. Linear fits miss the curvature; the quadratic-adjusted fit comes closer.

The two linear fits lie almost on top of each other, because the omitted-variable term cancels exactly in this DGP, and both run straight through the curvature. At $d = 2$ they give 4.13 against a true 4.40; at $d = 6$ they give 9.82 and 9.77 against a true 7.60. Their common slope of about 1.41 is close to the average of the true slope $2 - 0.3d$ over the observed doses, which is 1.403 — a straight line through a concave curve recovers the average derivative, not the curve.

The quadratic adjusted fit gives 4.34, 6.61 and 7.83 at $d = 2, 4, 6$ against the true 4.40, 6.60 and 7.60. It is better here because the true curve is quadratic and we fitted a quadratic. In real applications the functional form is not known, which is what motivates the next two sections.

10.3 Generalized propensity score (GPS)

The generalized propensity score is the conditional density of the treatment given covariates:

\[ r(d, x) = f_{D \mid X}(d \mid x). \tag{10.2}\]

Two results from Hirano and Imbens (2004) are what make this object useful, and both are worth stating before we compute anything.

The first is the assumption they need, which is weaker than it looks. Weak unconfoundedness requires $Y(d) \perp D \mid X$ separately for each dose $d$. It does not require the whole path $\{Y(d) : d\}$ to be jointly independent of $D$ given $X$, which would be a much stronger demand at a continuous treatment.

The second is the balancing property. Conditional on the scalar $r(d, X)$, the covariates carry no further information about whether a unit received dose $d$:

\[ f_D\bigl(d \mid r(d, X)\bigr) = f_D\bigl(d \mid r(d, X),\, X\bigr). \tag{10.3}\]

So a one-dimensional score does the job of the full covariate vector, exactly as in the binary case. Put the two together and the dose-response function is

\[ \mu(d) = E[Y(d)] = E\Bigl[\, E\bigl(Y \mid D = d,\; r(d, X)\bigr) \Bigr], \tag{10.4}\]

where the outer expectation runs over the distribution of $r(d, X)$ in the whole sample. Read the argument of $r$ carefully, because it is where this estimator is usually got wrong: at each counterfactual dose $d$ we must recompute every unit’s score at that $d$, not reuse the score at the dose the unit happened to receive. The code below does this in estimate_dose(), and it is the reason the GPS has to be recomputed inside the loop over doses.

One more point of theory explains why the GPS enters as a regressor rather than as a weight. The GPS is a conditional density, not a probability, so it carries the units of $1/D$; weights of the form $1/\hat r$ are therefore not scale-free and change if the dose is rescaled. Multiplying by the marginal density $f_D(D)$ restores scale invariance, which is what the stabilised weights below do, and it is why they are reported here only as an overlap diagnostic. The Hirano-Imbens estimator sidesteps the issue by putting the score in the outcome model instead.

In practice we start with a normal model for $D \mid X$: $D \mid X \sim \mathcal{N}(\beta_0+\beta'X,\sigma^2)$. Fitting $D$ on $X_1$ and $X_2$ by OLS gives $\hat\beta$ and $\hat\sigma$, and the score at any dose $d$ for unit $i$ is the normal density evaluated there,

\[ \hat r(d, x_i) = \frac{1}{\hat\sigma\sqrt{2\pi}} \exp\!\left\{-\frac{\bigl(d - x_i'\hat\beta\bigr)^2}{2\hat\sigma^2}\right\}. \tag{10.5}\]

That is the whole calculation: one linear regression, then one density evaluation per unit per dose. Nothing here needs a package. lm fits the dose model, dnorm evaluates Equation 10.5, and density and approx from stats supply the marginal density $f_D$ used in the diagnostic weights, which are trimmed at the 99th percentile.

Code

# Fit GPS model: D ~ X with Normal residuals
gps_model <- lm(D ~ X1 + X2, data = df)
d_pred    <- predict(gps_model)
sigma_e   <- sigma(gps_model)

# GPS for each unit's observed treatment
gps <- dnorm(df$D, mean = d_pred, sd = sigma_e)
df$gps <- gps

# IPW-style weights for continuous treatment: 1 / GPS
# Stabilise by multiplying by the marginal density f_D(D)
# NOTE: sw_cont is computed only as an overlap/positivity diagnostic.
# The Hirano-Imbens estimator below does NOT weight by it: the GPS
# enters as a regressor in the outcome model.
marginal_D <- density(df$D)
fD <- approx(marginal_D$x, marginal_D$y, xout = df$D)$y
df$sw_cont <- pmin(fD / gps, quantile(fD / gps, 0.99))

cat(sprintf("Weight summary: min=%.3f, mean=%.3f, max=%.3f\n",
            min(df$sw_cont), mean(df$sw_cont), max(df$sw_cont)))

Weight summary: min=0.039, mean=0.947, max=4.390

The weights run from 0.039 to 4.390 with a mean of 0.947. No unit carries an extreme weight, so overlap along the dose is good and the estimators below are not resting on a handful of observations.

The Hirano-Imbens estimator has two steps. First regress $Y$ on flexible functions of $D$ and the estimated GPS. Then, for each dose $d$, average the predicted outcome over the sample. The plot compares the resulting curve against the truth.

Code

# Stage 2: regress Y on flexible function of (D, GPS)
hi_fit <- lm(Y ~ poly(D, 2) + poly(gps, 2) + I(D * gps), data = df)

# Dose-response at each d: average over the empirical distribution of GPS
# at that d.
estimate_dose <- function(d) {
  # Predict GPS at this d for each unit
  gps_at_d <- dnorm(d, mean = d_pred, sd = sigma_e)
  nd <- data.frame(D = d, gps = gps_at_d)
  mean(predict(hi_fit, newdata = nd))
}

ds <- seq(0, 10, by = 0.5)
hi_dose <- sapply(ds, estimate_dose)

df_hi <- tibble(d = ds, mu = hi_dose, method = "Hirano-Imbens GPS")

ggplot() +
  geom_line(data = df_hi, aes(d, mu, colour = method), linewidth = 1.2) +
  geom_line(data = df_true, aes(d, mu), colour = "black",
            linetype = "dashed", linewidth = 1.0) +
  labs(x = "Dose D", y = "E[Y(d)]",
       caption = "Dashed black = truth; coloured = Hirano-Imbens GPS estimator") +
  theme_minimal()

The GPS curve tracks the concavity that the linear fits missed. It is fitted with a quadratic in $D$ and in the GPS, so it is not assumption-free either, but the assumption is now about how the outcome depends on the dose and the score rather than about a straight line.

10.4 Doubly-robust dose-response estimation

A kernel-localized AIPW estimator — in the spirit of Kennedy et al. (2017), who develop a more elaborate pseudo-outcome version with formal guarantees — is the continuous-treatment analogue of AIPW. It combines:

An outcome model $\hat\mu(d, x) = \mathbb{E}[Y \mid D = d, X = x]$
A generalised propensity score $\hat r(d \mid x)$

Then the dose-response at $d$ is estimated by

\[ \hat\mu_{DR}(d) = \frac{1}{n} \sum_i \left[ \hat\mu(d, X_i) + K_h(D_i - d) \cdot \frac{Y_i - \hat\mu(D_i, X_i)}{\hat r(D_i \mid X_i)} \right], \tag{10.6}\]

where $K_h$ is a kernel around dose $d$. The bandwidth $h$ controls how much nearby observed doses contribute to the estimate. (The code below uses the weight-normalized — Hájek — version of the correction term, dividing the weights by their mean; the population limit is unchanged.)

Two features of this estimator are worth understanding, because they are what the comparison table at the end of the chapter is really about.

The first is why it is doubly robust, and the algebra is the same as in the binary case worked through in Section 8.2.1. The first term is a plug-in built from the outcome model. The second subtracts that model’s own residual, reweighted by the score. A first-order error in $\hat\mu$ therefore enters twice with opposite signs and cancels, and what survives involves the product of the error in $\hat\mu$ and the error in $\hat r$. If either nuisance is right, the product is zero.

The second is why a bandwidth appears here when none appeared for the ATE. The ATE is a single number, and it is estimable at the $\sqrt n$ rate. But $\mu(d)$ at a point $d$ is a regression function, and no $\sqrt n$ estimator of it exists without smoothness assumptions, because the data hold only finitely many observations near any given dose. The kernel is how we borrow from neighbouring doses, and $h$ sets the exchange rate: a wide bandwidth pulls in doses where the truth is different, a narrow one leaves too few observations to average over. Under the usual second-order smoothness the best trade is $h \propto n^{-1/5}$ and the rate is $n^{-2/5}$, slower than $\sqrt n$ by construction. That is why bandwidth sensitivity is listed below as a cost of the method rather than as a tuning detail.

Kennedy et al. (2017) handle the localisation more carefully. Instead of smoothing the correction term, they construct a pseudo-outcome whose conditional mean given $D = d$ is exactly $\mu(d)$, and then run a single nonparametric regression of that pseudo-outcome on $D$. All the smoothing sits in one place, which is what lets them state formal guarantees. The version used here keeps the doubly-robust structure and is easier to read, with weaker theory behind it.

Using the GPS already estimated above, the implementation is short. The outcome model is a quadratic in $D$ with both covariates and their interactions with $D$; the kernel is Gaussian with bandwidth $h = 0.5$. The plot puts this curve next to the GPS curve and the truth.

Code

# Outcome model μ̂(D, X): flexible polynomial
mu_fit <- lm(Y ~ D + I(D^2) + X1 + X2 + I(D * X1) + I(D * X2), data = df)
df$mu_at_obs <- predict(mu_fit)

# DR estimator at dose d with bandwidth h
dr_estimate <- function(d, h = 0.5) {
  # Counterfactual prediction at D = d for each unit
  mu_at_d <- predict(mu_fit, newdata = transform(df, D = d))
  # Kernel weight on observed treatment near d
  kernel <- dnorm(df$D - d, sd = h)
  weight <- kernel / df$gps
  weight <- weight / mean(weight)  # normalise
  # DR correction
  correction <- weight * (df$Y - df$mu_at_obs)
  mean(mu_at_d) + mean(correction)
}

ds <- seq(0, 10, by = 0.5)
dr_dose <- sapply(ds, dr_estimate)

df_drf <- tibble(d = ds, mu = dr_dose, method = "DR (kernel AIPW)")

ggplot() +
  geom_line(data = df_hi,  aes(d, mu, colour = method), linewidth = 1.0) +
  geom_line(data = df_drf, aes(d, mu, colour = method), linewidth = 1.0) +
  geom_line(data = df_true, aes(d, mu), colour = "black",
            linetype = "dashed", linewidth = 1.0) +
  labs(x = "Dose D", y = "E[Y(d)]",
       caption = "Black dashed = truth") +
  theme_minimal()

Read this as an illustrative kernel-AIPW construction rather than a turnkey doubly robust theorem. For continuous treatments, a formal doubly robust dose-response result – consistency if either the outcome model or the GPS is correct – requires additional bandwidth, nuisance-convergence-rate, smoothness, and orthogonalization conditions (see Kennedy et al. on continuous-treatment DR estimators). The self-normalized kernel correction shown here conveys the intuition but does not by itself establish those guarantees. As usual, the bandwidth trades bias against variance.

The causaldrf package provides ready-made versions of these estimators, including hi_est, iptw_est, aipwee_est, gam_est, nw_est, and bart_est.

10.5 Longitudinal Modified Treatment Policies (LMTP)

Fixed-dose interventions can be unrealistic. For example, there may be almost no observations near $D=10$, so estimating $Y(10)$ requires extrapolation. LMTP avoids this by defining interventions as changes to the observed treatment, such as “shift everyone’s dose down by one unit.” This keeps the counterfactual treatment closer to the observed support.

10.5.1 What LMTP identifies

Write the policy as a function $d$ that takes a unit’s history $H_t$ up to time $t$ and returns a treatment value: $d(H_t) = A_t - 1$ for a one-unit downward shift, or $d(H_t) = 0$ for the static “treat nobody” policy. The estimand is $E[Y(\bar d)]$, the mean outcome if every period’s treatment were replaced by the value the policy assigns.

Three assumptions identify it. Sequential exchangeability: at each period the potential outcome is independent of the current treatment given the observed history, $Y(\bar d) \perp A_t \mid H_t$. Consistency, in the usual sense. And positivity, which is the assumption that motivates the whole approach: the shifted value must have positive density given the history. A shift of one unit stays inside the observed support, while a fixed dose of $D = 10$ may sit where no comparable unit exists.

Given those three, the estimand has a sequential-regression form. Start at the last period and work backwards:

\[ m_T(H_T) = E\bigl[Y \mid H_T,\, A_T = d(H_T)\bigr], \qquad m_t(H_t) = E\bigl[m_{t+1} \mid H_t,\, A_t = d(H_t)\bigr], \tag{10.7}\]

with $E[Y(\bar d)] = E[m_1]$. In plain terms: predict the outcome as if the last treatment had followed the policy, treat that prediction as the outcome for the previous period and predict it the same way, and repeat back to the first period. This is the g-formula written as a sequence of regressions rather than as an integral, which is what makes it estimable at all. lmtp fits each regression, optionally with machine learning, and adds a targeting or weighting step at each stage so that the recursion is doubly robust rather than a plug-in.

10.5.2 A two-period example

We simulate $N = 1500$ units over two periods. The first covariate is $L_1 \sim U(0,1)$ and the first treatment is a coin flip, $A_1 \sim \text{Bernoulli}(0.5)$. The second covariate is $L_2 \sim N(0.25L_1,\, 0.5^2)$, and the second treatment depends on both the first treatment and the current covariate, $A_2 \sim \text{Bernoulli}(\Lambda(0.5 - 0.2A_1 + 0.1L_2))$. The outcome is $Y \sim N(A_2 + L_2,\, 0.3^2)$.

The truths follow from the DGP without simulation. $L_2$ does not depend on $A_1$, so setting both treatments to 0 leaves $E[Y] = E[L_2] = 0.25 \cdot E[L_1] = 0.125$. Setting both to 1 gives $1 + 0.125 = 1.125$. And $A_1$ affects $Y$ only through $A_2$, which the policy fixes, so the contrast between the two policies is exactly 1.

Code

set.seed(1)
N    <- 1500
L_1  <- runif(N, 0, 1)
A_1  <- rbinom(N, 1, 0.5)
L_2  <- rnorm(N, 0.25 * L_1, 0.5)
A_2  <- rbinom(N, 1, plogis(0.5 - 0.2 * A_1 + 0.1 * L_2))
Y    <- rnorm(N, A_2 + L_2, 0.3)

lmtp_data <- data.frame(L_1, A_1, L_2, A_2, Y)
head(lmtp_data)

        L_1 A_1        L_2 A_2           Y
1 0.2655087   0  0.4913989   0 -0.42997401
2 0.3721239   1 -0.3696255   0 -0.76396535
3 0.5728534   0  0.5900039   0  0.79312079
4 0.9082078   1 -0.2434529   1  1.23628221
5 0.2016819   0  0.3198965   0  0.02620015
6 0.8983897   0  0.1336102   1  1.00441306

For simplicity the worked example uses static policies (set every treatment to 0, then to 1) on a binary treatment — lmtp handles these in the same framework. A genuinely modified policy on a continuous dose would instead return a shift of the observed value, e.g. function(data, trt) data[[trt]] - 1, which is what keeps the counterfactual on the observed support. Here is the set-to-zero policy:

Code

shift_zero <- function(data, trt) {
  # Set treatment to 0 everywhere — for comparison to "never treat"
  a <- data[[trt]]
  rep(0, length(a))
}

Then fit the LMTP estimator with lmtp::lmtp_tmle:

Code

# TMLE-based LMTP estimator with super-learners for nuisance functions
# (using SL.glm to keep it fast; in practice use richer SL libraries)
lmtp_fit <- lmtp_tmle(
  data       = lmtp_data,
  trt        = c("A_1", "A_2"),
  outcome    = "Y",
  baseline   = NULL,
  time_vary  = list(c("L_1"), c("L_2")),
  shift      = shift_zero,
  outcome_type = "continuous",
  learners_outcome = c("SL.glm"),
  learners_trt     = c("SL.glm"),
  folds = 2
)

# tidy() gives the estimate, SE, and CI as a data frame
# (the styled console print does not render cleanly in the book)
tidy(lmtp_fit)

# A tibble: 1 × 4
  estimate std.error conf.low conf.high
     <dbl>     <dbl>    <dbl>     <dbl>
1    0.118    0.0252   0.0681     0.167

The set-to-zero policy gives 0.118 with a standard error of 0.0252 and a 95% interval of $[0.068, 0.167]$, which covers the true 0.125.

The useful part of lmtp is that it handles the time-ordering of covariates, treatments, and outcomes: $L_2$ is a post-treatment covariate for $A_1$ and a pre-treatment covariate for $A_2$, and the estimator has to respect that ordering rather than conditioning on $L_2$ throughout.

10.5.3 Contrasts between policies

We can also compare two policies, such as “set to 0” versus “set to 1”:

Code

shift_one <- function(data, trt) {
  a <- data[[trt]]
  rep(1, length(a))
}

lmtp_one <- lmtp_tmle(
  data       = lmtp_data,
  trt        = c("A_1", "A_2"),
  outcome    = "Y",
  baseline   = NULL,
  time_vary  = list(c("L_1"), c("L_2")),
  shift      = shift_one,
  outcome_type = "continuous",
  learners_outcome = c("SL.glm"),
  learners_trt     = c("SL.glm"),
  folds = 2
)

contrast <- lmtp_contrast(lmtp_one, ref = lmtp_fit)
as.data.frame(contrast$estimates)

     shift       ref  estimate  std.error  conf.low conf.high       p.value
1 1.103091 0.1176079 0.9854826 0.03489166 0.9170962  1.053869 1.683133e-175

The contrast is the difference between the two policy-specific mean outcomes. Treat-always gives 1.103 against a true 1.125, treat-never gives 0.118 against a true 0.125, and their difference is 0.9855 with a standard error of 0.0349 and a 95% interval of $[0.917, 1.054]$. The true contrast is exactly 1, and the interval covers it.

10.6 Choosing among the three

Method	Strengths	Weaknesses
GPS (Hirano-Imbens)	Conceptually transparent; widely understood	Parametric assumptions on the GPS and the outcome-GPS regression; sensitive to their misspecification
Doubly-robust (kernel AIPW)	Consistent if either nuisance is correct, and can attain nonparametric oracle rates, under the bandwidth, smoothness, and nuisance-rate conditions the formal estimator requires	More complex; the construction shown above conveys the intuition rather than establishing those guarantees; nuisance models need careful tuning; sensitive to bandwidth
LMTP	Handles continuous + multivalued + time-varying uniformly; respects positivity	More setup; requires defining a shift function

In practice:

For a simple cross-sectional continuous treatment with good overlap, GPS or the Kennedy DR estimator is a natural start.
For time-varying treatments or policy shifts, use LMTP.
For multivalued treatments, use a multinomial propensity model and the same adjustment logic.

10.7 When to use these methods

These methods are useful when:

Treatment has more than two values.
The research question is about the shape of the dose-response curve.
The intervention is better described as shifting treatment than setting treatment to a fixed value.

For strictly binary, single-time-period problems, the standard Estimation and Heterogeneous Effects methods are sufficient.

10.8 Summary

Continuous and multivalued treatments usually require a dose-response curve, not a single coefficient.
GPS extends propensity-score logic to continuous treatments by modeling the density of treatment given covariates.
Kennedy’s DR estimator adds an outcome model to the GPS.
LMTP is useful when the intervention is a policy shift or when treatment is time-varying.

# Continuous and Multivalued Treatments ```{r} #| include: false library(tidyverse) library(WeightIt) library(causaldrf) library(lmtp) library(SuperLearner) library(broom) ``` ```{r} #| eval: false library(tidyverse) library(WeightIt) library(causaldrf) library(lmtp) library(SuperLearner) ``` So far most examples use a binary treatment. Many treatments are not binary: dose of a drug, hours in a program, dollars spent, or a treatment with several levels. In those cases "the treatment effect" is not one number. We usually want a dose-response curve, $$ \mu(d) = \mathbb{E}[Y(d)], \quad d \in \mathcal{D}. $$ {#eq-continuous-treatments-1} I use three approaches: 1. **Generalized propensity score** (GPS; Hirano-Imbens 2004), which models the conditional density of treatment given covariates. 2. **Doubly-robust dose-response estimation** [in the spirit of @kennedy-2017-continuous], which combines a GPS with an outcome model. 3. **Longitudinal modified treatment policies** [LMTP; @diaz-2023-lmtp], which are useful when treatments are time-varying or when fixed treatment levels create positivity problems. > **Related reading**: The LMTP material below is condensed from the > longer treatment in [Longitudinal modified treatment policy (LMTP)](https://xiangao.github.io/blog_book/lmtp.html) > in *Topics on Econometrics and Causal Inference*, which follows > Nicholas Williams's tutorials at > [beyondtheate.com](https://beyondtheate.com/). ## The continuous-treatment setup The aim is to recover a known dose-response curve from data in which the dose is chosen partly on covariates that also move the outcome. We simulate $n = 2000$ units with two covariates, $X_1, X_2 \sim N(0,1)$, independent of each other. The dose is $D = 2 + 0.5X_1 + 0.5X_2 + \nu$ with $\nu \sim N(0,1)$. The dose-response curve is $\mu(d) = 1 + 2d - 0.15d^2$, and the outcome is $Y = \mu(D) + 0.3X_1 - 0.3X_2 + \varepsilon$ with $\varepsilon \sim N(0,1)$. ```{r} #| label: dgp #| cache: true set.seed(42) n <- 2000 X1 <- rnorm(n) X2 <- rnorm(n) # Continuous treatment, depends on X1, X2 D <- 2 + 0.5 * X1 + 0.5 * X2 + rnorm(n, sd = 1.0) # Dose-response: nonlinear in D true_mu <- function(d) 1.0 + 2.0 * d - 0.15 * d^2 Y <- true_mu(D) + 0.3 * X1 - 0.3 * X2 + rnorm(n) df <- tibble(Y = Y, D = D, X1 = X1, X2 = X2) cat(sprintf("Treatment range: [%.2f, %.2f]\n", min(D), max(D))) cat(sprintf("True μ(d=2) = %.2f, μ(d=4) = %.2f, μ(d=6) = %.2f\n", true_mu(2), true_mu(4), true_mu(6))) ``` The observed doses run from $-1.82$ to 6.39, and the true curve gives $\mu(2) = 4.40$, $\mu(4) = 6.60$ and $\mu(6) = 7.60$. The dose-response curve is concave, with clearly diminishing returns. It is worth being precise about the shape rather than calling it hump-shaped, because the vertex sits at $d = 2/0.3 \approx 6.67$, just *beyond* the largest dose in the sample. Over the observed range the curve is therefore increasing throughout --- its slope $2 - 0.3d$ is still positive (about 0.08) at the maximum observed dose --- so the downturn is extrapolation, not something these data identify. The estimation problem is to recover the curvature that is in range when the observed dose is confounded by $X$. ## Naive regression vs adjusted regression We fit three models and compare their implied dose-response curves against the truth: $Y$ on $D$ alone, $Y$ on $D$ and the covariates, and $Y$ on a quadratic in $D$ and the covariates. All three are evaluated at $X_1 = X_2 = 0$ over $d \in [0,10]$. One feature of this DGP is worth stating before reading the plot. $X_1$ and $X_2$ are confounders in the graph, since each raises the dose and each moves the outcome. But their linear bias cancels: each raises $D$ by 0.5, while they enter $Y$ with $+0.3$ and $-0.3$, so the omitted-variable term $\text{Cov}(D,\, 0.3X_1 - 0.3X_2)/\text{Var}(D)$ is $0.3(0.5) - 0.3(0.5) = 0$ exactly. The naive slope is 1.421 and the adjusted slope 1.409. Adjusting for $X$ changes almost nothing here, and the one problem left for the linear models is functional form. ```{r} #| label: naive #| cache: true # Linear fit, ignoring covariates linear_naive <- lm(Y ~ D, data = df) # Linear fit with covariate adjustment linear_adj <- lm(Y ~ D + X1 + X2, data = df) # Nonlinear fit with covariate adjustment poly_adj <- lm(Y ~ poly(D, 2) + X1 + X2, data = df) ds <- seq(0, 10, by = 0.5) predict_mu <- function(fit, ds) { data.frame(d = ds, X1 = 0, X2 = 0) |> rename(D = d) -> nd predict(fit, newdata = nd) } # True dose-response df_true <- tibble(d = ds, mu = true_mu(ds)) df_fits <- bind_rows( tibble(d = ds, mu = predict_mu(linear_naive, ds), method = "Linear (naive)"), tibble(d = ds, mu = predict_mu(linear_adj, ds), method = "Linear (adj X)"), tibble(d = ds, mu = predict_mu(poly_adj, ds), method = "Quadratic (adj X)") ) ``` ```{r} #| label: naive-plot #| fig-cap: "Fitted vs true dose-response. Linear fits miss the curvature; the quadratic-adjusted fit comes closer." #| fig-width: 7 #| fig-height: 4 ggplot() + geom_line(data = df_fits, aes(d, mu, colour = method), linewidth = 1.0) + geom_line(data = df_true, aes(d, mu), colour = "black", linetype = "dashed", linewidth = 1.0) + labs(x = "Dose D", y = "Expected outcome E[Y(d)]", colour = "Method", caption = "Black dashed line = true μ(d) = 1 + 2d - 0.15 d²") + theme_minimal() ``` The two linear fits lie almost on top of each other, because the omitted-variable term cancels exactly in this DGP, and both run straight through the curvature. At $d = 2$ they give 4.13 against a true 4.40; at $d = 6$ they give 9.82 and 9.77 against a true 7.60. Their common slope of about 1.41 is close to the average of the true slope $2 - 0.3d$ over the observed doses, which is 1.403 --- a straight line through a concave curve recovers the average derivative, not the curve. The quadratic adjusted fit gives 4.34, 6.61 and 7.83 at $d = 2, 4, 6$ against the true 4.40, 6.60 and 7.60. It is better here because the true curve is quadratic and we fitted a quadratic. In real applications the functional form is not known, which is what motivates the next two sections. ## Generalized propensity score (GPS) The **generalized propensity score** is the conditional density of the treatment given covariates: $$ r(d, x) = f_{D \mid X}(d \mid x). $$ {#eq-continuous-treatments-2} Two results from @hirano-imbens-2004-gps are what make this object useful, and both are worth stating before we compute anything. The first is the assumption they need, which is weaker than it looks. *Weak unconfoundedness* requires $Y(d) \perp D \mid X$ separately for each dose $d$. It does not require the whole path $\{Y(d) : d\}$ to be jointly independent of $D$ given $X$, which would be a much stronger demand at a continuous treatment. The second is the balancing property. Conditional on the scalar $r(d, X)$, the covariates carry no further information about whether a unit received dose $d$: $$ f_D\bigl(d \mid r(d, X)\bigr) = f_D\bigl(d \mid r(d, X),\, X\bigr). $$ {#eq-continuous-treatments-balance} So a one-dimensional score does the job of the full covariate vector, exactly as in the binary case. Put the two together and the dose-response function is $$ \mu(d) = E[Y(d)] = E\Bigl[\, E\bigl(Y \mid D = d,\; r(d, X)\bigr) \Bigr], $$ {#eq-continuous-treatments-doseresp} where the outer expectation runs over the distribution of $r(d, X)$ in the whole sample. Read the argument of $r$ carefully, because it is where this estimator is usually got wrong: at each counterfactual dose $d$ we must recompute every unit's score *at that* $d$, not reuse the score at the dose the unit happened to receive. The code below does this in `estimate_dose()`, and it is the reason the GPS has to be recomputed inside the loop over doses. One more point of theory explains why the GPS enters as a regressor rather than as a weight. The GPS is a conditional *density*, not a probability, so it carries the units of $1/D$; weights of the form $1/\hat r$ are therefore not scale-free and change if the dose is rescaled. Multiplying by the marginal density $f_D(D)$ restores scale invariance, which is what the stabilised weights below do, and it is why they are reported here only as an overlap diagnostic. The Hirano-Imbens estimator sidesteps the issue by putting the score in the outcome model instead. In practice we start with a normal model for $D \mid X$: $D \mid X \sim \mathcal{N}(\beta_0+\beta'X,\sigma^2)$. Fitting $D$ on $X_1$ and $X_2$ by OLS gives $\hat\beta$ and $\hat\sigma$, and the score at any dose $d$ for unit $i$ is the normal density evaluated there, $$ \hat r(d, x_i) = \frac{1}{\hat\sigma\sqrt{2\pi}} \exp\!\left\{-\frac{\bigl(d - x_i'\hat\beta\bigr)^2}{2\hat\sigma^2}\right\}. $$ {#eq-continuous-treatments-gpshat} That is the whole calculation: one linear regression, then one density evaluation per unit per dose. Nothing here needs a package. `lm` fits the dose model, `dnorm` evaluates @eq-continuous-treatments-gpshat, and `density` and `approx` from `stats` supply the marginal density $f_D$ used in the diagnostic weights, which are trimmed at the 99th percentile. ```{r} #| label: gps #| cache: true # Fit GPS model: D ~ X with Normal residuals gps_model <- lm(D ~ X1 + X2, data = df) d_pred <- predict(gps_model) sigma_e <- sigma(gps_model) # GPS for each unit's observed treatment gps <- dnorm(df$D, mean = d_pred, sd = sigma_e) df$gps <- gps # IPW-style weights for continuous treatment: 1 / GPS # Stabilise by multiplying by the marginal density f_D(D) # NOTE: sw_cont is computed only as an overlap/positivity diagnostic. # The Hirano-Imbens estimator below does NOT weight by it: the GPS # enters as a regressor in the outcome model. marginal_D <- density(df$D) fD <- approx(marginal_D$x, marginal_D$y, xout = df$D)$y df$sw_cont <- pmin(fD / gps, quantile(fD / gps, 0.99)) cat(sprintf("Weight summary: min=%.3f, mean=%.3f, max=%.3f\n", min(df$sw_cont), mean(df$sw_cont), max(df$sw_cont))) ``` The weights run from 0.039 to 4.390 with a mean of 0.947. No unit carries an extreme weight, so overlap along the dose is good and the estimators below are not resting on a handful of observations. The Hirano-Imbens estimator has two steps. First regress $Y$ on flexible functions of $D$ and the estimated GPS. Then, for each dose $d$, average the predicted outcome over the sample. The plot compares the resulting curve against the truth. ```{r} #| label: hi-est #| cache: true # Stage 2: regress Y on flexible function of (D, GPS) hi_fit <- lm(Y ~ poly(D, 2) + poly(gps, 2) + I(D * gps), data = df) # Dose-response at each d: average over the empirical distribution of GPS # at that d. estimate_dose <- function(d) { # Predict GPS at this d for each unit gps_at_d <- dnorm(d, mean = d_pred, sd = sigma_e) nd <- data.frame(D = d, gps = gps_at_d) mean(predict(hi_fit, newdata = nd)) } ds <- seq(0, 10, by = 0.5) hi_dose <- sapply(ds, estimate_dose) df_hi <- tibble(d = ds, mu = hi_dose, method = "Hirano-Imbens GPS") ggplot() + geom_line(data = df_hi, aes(d, mu, colour = method), linewidth = 1.2) + geom_line(data = df_true, aes(d, mu), colour = "black", linetype = "dashed", linewidth = 1.0) + labs(x = "Dose D", y = "E[Y(d)]", caption = "Dashed black = truth; coloured = Hirano-Imbens GPS estimator") + theme_minimal() ``` The GPS curve tracks the concavity that the linear fits missed. It is fitted with a quadratic in $D$ and in the GPS, so it is not assumption-free either, but the assumption is now about how the outcome depends on the dose and the score rather than about a straight line. ## Doubly-robust dose-response estimation A kernel-localized AIPW estimator — in the spirit of @kennedy-2017-continuous, who develop a more elaborate pseudo-outcome version with formal guarantees — is the continuous-treatment analogue of AIPW. It combines: - An outcome model $\hat\mu(d, x) = \mathbb{E}[Y \mid D = d, X = x]$ - A generalised propensity score $\hat r(d \mid x)$ Then the dose-response at $d$ is estimated by $$ \hat\mu_{DR}(d) = \frac{1}{n} \sum_i \left[ \hat\mu(d, X_i) + K_h(D_i - d) \cdot \frac{Y_i - \hat\mu(D_i, X_i)}{\hat r(D_i \mid X_i)} \right], $$ {#eq-continuous-treatments-3} where $K_h$ is a kernel around dose $d$. The bandwidth $h$ controls how much nearby observed doses contribute to the estimate. (The code below uses the weight-normalized — Hájek — version of the correction term, dividing the weights by their mean; the population limit is unchanged.) Two features of this estimator are worth understanding, because they are what the comparison table at the end of the chapter is really about. The first is why it is doubly robust, and the algebra is the same as in the binary case worked through in @sec-plug-in-bias. The first term is a plug-in built from the outcome model. The second subtracts that model's own residual, reweighted by the score. A first-order error in $\hat\mu$ therefore enters twice with opposite signs and cancels, and what survives involves the product of the error in $\hat\mu$ and the error in $\hat r$. If either nuisance is right, the product is zero. The second is why a bandwidth appears here when none appeared for the ATE. The ATE is a single number, and it is estimable at the $\sqrt n$ rate. But $\mu(d)$ at a *point* $d$ is a regression function, and no $\sqrt n$ estimator of it exists without smoothness assumptions, because the data hold only finitely many observations near any given dose. The kernel is how we borrow from neighbouring doses, and $h$ sets the exchange rate: a wide bandwidth pulls in doses where the truth is different, a narrow one leaves too few observations to average over. Under the usual second-order smoothness the best trade is $h \propto n^{-1/5}$ and the rate is $n^{-2/5}$, slower than $\sqrt n$ by construction. That is why bandwidth sensitivity is listed below as a cost of the method rather than as a tuning detail. @kennedy-2017-continuous handle the localisation more carefully. Instead of smoothing the correction term, they construct a pseudo-outcome whose conditional mean given $D = d$ is exactly $\mu(d)$, and then run a single nonparametric regression of that pseudo-outcome on $D$. All the smoothing sits in one place, which is what lets them state formal guarantees. The version used here keeps the doubly-robust structure and is easier to read, with weaker theory behind it. Using the GPS already estimated above, the implementation is short. The outcome model is a quadratic in $D$ with both covariates and their interactions with $D$; the kernel is Gaussian with bandwidth $h = 0.5$. The plot puts this curve next to the GPS curve and the truth. ```{r} #| label: drf #| cache: true # Outcome model μ̂(D, X): flexible polynomial mu_fit <- lm(Y ~ D + I(D^2) + X1 + X2 + I(D * X1) + I(D * X2), data = df) df$mu_at_obs <- predict(mu_fit) # DR estimator at dose d with bandwidth h dr_estimate <- function(d, h = 0.5) { # Counterfactual prediction at D = d for each unit mu_at_d <- predict(mu_fit, newdata = transform(df, D = d)) # Kernel weight on observed treatment near d kernel <- dnorm(df$D - d, sd = h) weight <- kernel / df$gps weight <- weight / mean(weight) # normalise # DR correction correction <- weight * (df$Y - df$mu_at_obs) mean(mu_at_d) + mean(correction) } ds <- seq(0, 10, by = 0.5) dr_dose <- sapply(ds, dr_estimate) df_drf <- tibble(d = ds, mu = dr_dose, method = "DR (kernel AIPW)") ggplot() + geom_line(data = df_hi, aes(d, mu, colour = method), linewidth = 1.0) + geom_line(data = df_drf, aes(d, mu, colour = method), linewidth = 1.0) + geom_line(data = df_true, aes(d, mu), colour = "black", linetype = "dashed", linewidth = 1.0) + labs(x = "Dose D", y = "E[Y(d)]", caption = "Black dashed = truth") + theme_minimal() ``` Read this as an illustrative kernel-AIPW construction rather than a turnkey doubly robust theorem. For continuous treatments, a formal doubly robust dose-response result -- consistency if either the outcome model or the GPS is correct -- requires additional bandwidth, nuisance-convergence-rate, smoothness, and orthogonalization conditions (see Kennedy et al. on continuous-treatment DR estimators). The self-normalized kernel correction shown here conveys the intuition but does not by itself establish those guarantees. As usual, the bandwidth trades bias against variance. The `causaldrf` package provides ready-made versions of these estimators, including `hi_est`, `iptw_est`, `aipwee_est`, `gam_est`, `nw_est`, and `bart_est`. ## Longitudinal Modified Treatment Policies (LMTP) Fixed-dose interventions can be unrealistic. For example, there may be almost no observations near $D=10$, so estimating $Y(10)$ requires extrapolation. LMTP avoids this by defining interventions as changes to the observed treatment, such as "shift everyone's dose down by one unit." This keeps the counterfactual treatment closer to the observed support. ### What LMTP identifies Write the policy as a function $d$ that takes a unit's history $H_t$ up to time $t$ and returns a treatment value: $d(H_t) = A_t - 1$ for a one-unit downward shift, or $d(H_t) = 0$ for the static "treat nobody" policy. The estimand is $E[Y(\bar d)]$, the mean outcome if every period's treatment were replaced by the value the policy assigns. Three assumptions identify it. Sequential exchangeability: at each period the potential outcome is independent of the current treatment given the observed history, $Y(\bar d) \perp A_t \mid H_t$. Consistency, in the usual sense. And positivity, which is the assumption that motivates the whole approach: the *shifted* value must have positive density given the history. A shift of one unit stays inside the observed support, while a fixed dose of $D = 10$ may sit where no comparable unit exists. Given those three, the estimand has a sequential-regression form. Start at the last period and work backwards: $$ m_T(H_T) = E\bigl[Y \mid H_T,\, A_T = d(H_T)\bigr], \qquad m_t(H_t) = E\bigl[m_{t+1} \mid H_t,\, A_t = d(H_t)\bigr], $$ {#eq-continuous-treatments-lmtp} with $E[Y(\bar d)] = E[m_1]$. In plain terms: predict the outcome as if the last treatment had followed the policy, treat that prediction as the outcome for the previous period and predict it the same way, and repeat back to the first period. This is the g-formula written as a sequence of regressions rather than as an integral, which is what makes it estimable at all. `lmtp` fits each regression, optionally with machine learning, and adds a targeting or weighting step at each stage so that the recursion is doubly robust rather than a plug-in. ### A two-period example We simulate $N = 1500$ units over two periods. The first covariate is $L_1 \sim U(0,1)$ and the first treatment is a coin flip, $A_1 \sim \text{Bernoulli}(0.5)$. The second covariate is $L_2 \sim N(0.25L_1,\, 0.5^2)$, and the second treatment depends on both the first treatment and the current covariate, $A_2 \sim \text{Bernoulli}(\Lambda(0.5 - 0.2A_1 + 0.1L_2))$. The outcome is $Y \sim N(A_2 + L_2,\, 0.3^2)$. The truths follow from the DGP without simulation. $L_2$ does not depend on $A_1$, so setting both treatments to 0 leaves $E[Y] = E[L_2] = 0.25 \cdot E[L_1] = 0.125$. Setting both to 1 gives $1 + 0.125 = 1.125$. And $A_1$ affects $Y$ only through $A_2$, which the policy fixes, so the contrast between the two policies is exactly 1. ```{r} #| label: lmtp-sim #| cache: true set.seed(1) N <- 1500 L_1 <- runif(N, 0, 1) A_1 <- rbinom(N, 1, 0.5) L_2 <- rnorm(N, 0.25 * L_1, 0.5) A_2 <- rbinom(N, 1, plogis(0.5 - 0.2 * A_1 + 0.1 * L_2)) Y <- rnorm(N, A_2 + L_2, 0.3) lmtp_data <- data.frame(L_1, A_1, L_2, A_2, Y) head(lmtp_data) ``` For simplicity the worked example uses *static* policies (set every treatment to 0, then to 1) on a binary treatment — `lmtp` handles these in the same framework. A genuinely *modified* policy on a continuous dose would instead return a shift of the observed value, e.g. `function(data, trt) data[[trt]] - 1`, which is what keeps the counterfactual on the observed support. Here is the set-to-zero policy: ```{r} #| label: shift-fn shift_zero <- function(data, trt) { # Set treatment to 0 everywhere — for comparison to "never treat" a <- data[[trt]] rep(0, length(a)) } ``` Then fit the LMTP estimator with `lmtp::lmtp_tmle`: ```{r} #| label: lmtp-fit #| cache: true #| warning: false # TMLE-based LMTP estimator with super-learners for nuisance functions # (using SL.glm to keep it fast; in practice use richer SL libraries) lmtp_fit <- lmtp_tmle( data = lmtp_data, trt = c("A_1", "A_2"), outcome = "Y", baseline = NULL, time_vary = list(c("L_1"), c("L_2")), shift = shift_zero, outcome_type = "continuous", learners_outcome = c("SL.glm"), learners_trt = c("SL.glm"), folds = 2 ) # tidy() gives the estimate, SE, and CI as a data frame # (the styled console print does not render cleanly in the book) tidy(lmtp_fit) ``` The set-to-zero policy gives 0.118 with a standard error of 0.0252 and a 95% interval of $[0.068, 0.167]$, which covers the true 0.125. The useful part of `lmtp` is that it handles the time-ordering of covariates, treatments, and outcomes: $L_2$ is a post-treatment covariate for $A_1$ and a pre-treatment covariate for $A_2$, and the estimator has to respect that ordering rather than conditioning on $L_2$ throughout. ### Contrasts between policies We can also compare two policies, such as "set to 0" versus "set to 1": ```{r} #| label: lmtp-contrast #| cache: true #| warning: false shift_one <- function(data, trt) { a <- data[[trt]] rep(1, length(a)) } lmtp_one <- lmtp_tmle( data = lmtp_data, trt = c("A_1", "A_2"), outcome = "Y", baseline = NULL, time_vary = list(c("L_1"), c("L_2")), shift = shift_one, outcome_type = "continuous", learners_outcome = c("SL.glm"), learners_trt = c("SL.glm"), folds = 2 ) contrast <- lmtp_contrast(lmtp_one, ref = lmtp_fit) as.data.frame(contrast$estimates) ``` The contrast is the difference between the two policy-specific mean outcomes. Treat-always gives 1.103 against a true 1.125, treat-never gives 0.118 against a true 0.125, and their difference is 0.9855 with a standard error of 0.0349 and a 95% interval of $[0.917, 1.054]$. The true contrast is exactly 1, and the interval covers it. ## Choosing among the three | Method | Strengths | Weaknesses | |---|---|---| | **GPS (Hirano-Imbens)** | Conceptually transparent; widely understood | Parametric assumptions on the GPS and the outcome-GPS regression; sensitive to their misspecification | | **Doubly-robust** (kernel AIPW) | Consistent if either nuisance is correct, and can attain nonparametric oracle rates, under the bandwidth, smoothness, and nuisance-rate conditions the formal estimator requires | More complex; the construction shown above conveys the intuition rather than establishing those guarantees; nuisance models need careful tuning; sensitive to bandwidth | | **LMTP** | Handles continuous + multivalued + time-varying uniformly; respects positivity | More setup; requires defining a shift function | In practice: - For a simple cross-sectional continuous treatment with good overlap, GPS or the Kennedy DR estimator is a natural start. - For time-varying treatments or policy shifts, use LMTP. - For multivalued treatments, use a multinomial propensity model and the same adjustment logic. ## When to use these methods These methods are useful when: 1. Treatment has more than two values. 2. The research question is about the shape of the dose-response curve. 3. The intervention is better described as shifting treatment than setting treatment to a fixed value. For strictly binary, single-time-period problems, the standard [Estimation](estimation.qmd) and [Heterogeneous Effects](heterogeneous-effects.qmd) methods are sufficient. ## Summary - Continuous and multivalued treatments usually require a dose-response curve, not a single coefficient. - GPS extends propensity-score logic to continuous treatments by modeling the density of treatment given covariates. - Kennedy's DR estimator adds an outcome model to the GPS. - LMTP is useful when the intervention is a policy shift or when treatment is time-varying.