28  Policy Evaluation and Dynamic Treatment Strategies

Status: Draft

v0.5

28.1 Introduction

This chapter focuses on Doing in the Observable stratum, how can we evaluate policies and treatment strategies? A dynamic treatment strategy (policy) is a rule that maps history to treatment, and we need methods to evaluate these policies using causal models.

28.2 Dynamic Treatment Strategies

28.2.1 Definition

A dynamic treatment strategy (policy) is a rule that maps history to treatment:

\[ A_t = \pi(H_t) \]

where \(H_t = (Y_{1:t}, A_{1:t-1}, L_{1:t})\) is the observed history.

28.2.2 Examples

  • “Treat if severity > threshold”: \(A_t = \mathbb{1}(L_t > \tau)\)
  • “Stop if recovered”: \(A_t = \mathbb{1}(Y_{t-1} = 0)\)
  • “Adaptive dosing”: \(A_t = f(Y_{1:t})\)

28.3 Evaluating Policies

28.3.1 The Policy Evaluation Problem

Question: What is the expected outcome under policy \(\pi\)?

\[ \mathbb{E}[Y^{do(\pi)}] = \mathbb{E}[Y \mid \text{policy } \pi] \]

28.3.2 Answer Using CDM

From a fitted CDM, compute: \[ \mathbb{E}^{do(\pi)}[Y_T] = \int Y_T \, P^{do(\pi)}(Y_T \mid Y_1) \, dY_T \]

Computation:

  1. Simulate trajectories under \(do(\pi)\)
  2. Average over trajectories
  3. Account for uncertainty

28.3.3 Methods for Policy Evaluation

G-computation: Simulate outcomes under policy

  • Fit outcome and confounder models
  • Simulate trajectories under policy
  • Average over simulated outcomes

IPTW: Weight observations by inverse probability of following policy

  • Fit treatment models
  • Compute weights: \(w_i = \prod_{t=1}^{T} \frac{1}{P(A_t = \pi(H_t) \mid H_t)}\)
  • Fit outcome model on weighted data

TMLE: Targeted estimation of policy effects

  • Fit initial outcome model
  • Fit treatment model for weights
  • Targeted update to estimate policy effect
  • Provides valid inference

28.3.4 Implementation: Policy Evaluation

Here’s an example evaluating different policies using G-computation:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))

@auto_using Random Distributions OrdinaryDiffEq CairoMakie

Random.seed!(42)

# Example: Disease management policies
# Policy 1: Treat if severity > 0.6
# Policy 2: Treat if severity > 0.4
# Policy 3: Always treat

β = 0.2  # Disease progression
α = 0.3  # Treatment effectiveness

# Policy definitions
policy1 = (X, t) -> X > 0.6 ? 1.0 : 0.0  # Treat if severe
policy2 = (X, t) -> X > 0.4 ? 1.0 : 0.0  # Treat if moderate
policy3 = (X, t) -> 1.0  # Always treat

function disease_policy!(du, u, p, t)
    """Disease model with policy function: treatment determined by policy(X, t)."""
    X = u[1]
    policy = p[1]
    A = policy(X, t)
    du[1] = β * X - α * A * X
end

u0 = [0.3]
tspan = (0.0, 20.0)

# Evaluate each policy
policies = [("Treat if X > 0.6", policy1), 
            ("Treat if X > 0.4", policy2),
            ("Always treat", policy3)]
outcomes = Float64[]
trajectories = []

for (name, policy) in policies
    prob = ODEProblem(disease_policy!, u0, tspan, (policy,))  # Using tuple for better performance
    sol = solve(prob, Tsit5())
    push!(outcomes, sol.u[end][1])
    push!(trajectories, (name=name, sol=sol))
end

println("Policy evaluation results:")
for (i, (name, _)) in enumerate(policies)
    println("  ", name, ": Final severity = ", round(outcomes[i], digits=3))
end
best_policy = argmin(outcomes)
println("\nBest policy: ", policies[best_policy][1], " (lowest final severity)")
Policy evaluation results:
  Treat if X > 0.6: Final severity = 0.603
  Treat if X > 0.4: Final severity = 0.402
  Always treat: Final severity = 0.041

Best policy: Always treat (lowest final severity)
Figure 28.1: Policy evaluation: comparing different treatment strategies

28.3.5 Value Function

The value function of policy \(\pi\) is: \[ V^\pi = \mathbb{E}^{do(\pi)}\left[\sum_{t=1}^T \gamma^t R_t\right] \]

where \(R_t\) is the reward at time \(t\) and \(\gamma\) is the discount factor.

28.3.6 Comparing Policies

Compare expected outcomes under different policies:

\[ \mathbb{E}[Y^{do(\pi_1)}] \quad \text{vs} \quad \mathbb{E}[Y^{do(\pi_2)}] \]

This allows us to choose the best policy for a given objective.

28.4 Optimal Policy

28.4.1 Definition

The optimal policy maximises expected outcome:

\[ \pi^* = \arg\max_\pi \mathbb{E}^{do(\pi)}[Y_T] \]

28.4.2 Finding Optimal Policies

Methods (Murphy 2003; Schulam and Saria 2017):

28.5 Off-Policy Evaluation

28.5.1 The Problem

Evaluate policy \(\pi\) using data collected under different policy \(\pi_0\).

28.5.2 Methods (Tennenholtz et al. 2020; Sutton and Barto 2018):

  • Importance sampling: Reweight observations by \(\frac{\pi(A \mid H)}{\pi_0(A \mid H)}\)
  • Doubly robust methods: Combine model-based and importance sampling (e.g., TMLE, augmented IPTW)
  • Causal methods: Use CDM to simulate under \(\pi\)

28.6 Off-Policy Evaluation

28.6.1 The Problem

Evaluate policy \(\pi\) using data collected under different policy \(\pi_0\).

28.6.2 Methods (Tennenholtz et al. 2020; Sutton and Barto 2018):

  • Importance sampling: Reweight observations by \(\frac{\pi(A \mid H)}{\pi_0(A \mid H)}\)
  • Doubly robust methods: Combine model-based and importance sampling (e.g., TMLE, augmented IPTW)
  • Causal methods: Use CDM to simulate under \(\pi\)

28.6.3 Targeted Maximum Likelihood Estimation (TMLE) for Off-Policy Evaluation

TMLE provides a robust approach to off-policy evaluation that combines the strengths of importance sampling and model-based methods (Laan and Rubin 2006; Laan and Rose 2011).

28.6.3.1 TMLE Algorithm for Off-Policy Evaluation

Goal: Estimate \(E[Y^{\pi}]\) (expected outcome under policy \(\pi\)) using data collected under policy \(\pi_0\).

Steps:

  1. Fit initial outcome model: \(Q_0(A, H) = E[Y \mid A, H]\) using data from \(\pi_0\)
  2. Fit behaviour policy model: \(g_0(A \mid H) = \pi_0(A \mid H)\) (known or estimated)
  3. Targeted update: Update outcome model to target \(E[Y^{\pi}]\):
  • Compute clever covariate: \(H(A, H) = \frac{\pi(A \mid H)}{g_0(A \mid H)}\)
  • Fit logistic regression: \(\text{logit}(Q_1(A, H)) = \text{logit}(Q_0(A, H)) + \epsilon H(A, H)\)
  • Update: \(Q_1(A, H) = Q_0(A, H) + \epsilon H(A, H)\)
  1. Compute policy value: \[ \hat{E}[Y^{\pi}] = \frac{1}{n} \sum_{i=1}^n \sum_a \pi(a \mid H_i) Q_1(a, H_i) \]

Advantages:

  • Double robust: Consistent if either outcome model or behaviour policy model is correct
  • Semiparametric efficient: Achieves optimal variance
  • Valid inference: Provides confidence intervals
  • Handles policy mismatch: Works even when \(\pi\) and \(\pi_0\) differ substantially

Limitations:

  • Requires overlap: \(\pi(a \mid H) > 0\) implies \(\pi_0(a \mid H) > 0\) (positivity)
  • More complex than simple importance sampling

28.6.4 Implementation: Off-Policy Evaluation with Importance Sampling

Here’s an example of off-policy evaluation using importance sampling:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using Random Distributions GLM DataFrames CairoMakie

Random.seed!(42)

# Simulate data collected under behaviour policy π₀
# π₀: Treat if severity > 0.5
n = 500
L = rand(Uniform(0, 1), n)  # Severity
A_behavior = Int.(L .> 0.5)  # Behavior policy: treat if L > 0.5
Y = 0.5 .* A_behavior .- 0.3 .* L .+ rand(Normal(0, 0.1), n)  # Outcome

df = DataFrame(L = L, A = A_behavior, Y = Y)

# Target policy π: Treat if severity > 0.3 (more aggressive)
π_target = (L) -> L > 0.3 ? 1 : 0
π_behavior = (L) -> L > 0.5 ? 1 : 0

# Fit behaviour policy model g₀(A | L) = π₀(A | L)
g0_model = glm(@formula(A ~ L), df, Binomial(), LogitLink())

# Compute importance sampling weights
# w = π(A | L) / π₀(A | L)
weights = Float64[]
for i in 1:n
    L_i = L[i]
    A_i = A_behavior[i]
    
    # π₀(A | L) from fitted model
    p_behavior = predict(g0_model, DataFrame(L = [L_i]))[1]
    π0_A = A_i == 1 ? p_behavior : (1 - p_behavior)
    
    # π(A | L) from target policy
    # Target policy is deterministic: π_target(L) returns 0 or 1
    # Probability of observed action A_i under target policy
    π_target_prob = π_target(L_i)  # This is 0 or 1 (deterministic policy)
    π_A = A_i == 1 ? π_target_prob : (1 - π_target_prob)
    
    # Weight (with stabilisation to avoid extreme values)
    # Avoid division by zero
    if π0_A > 1e-10
        w = π_A / π0_A
    else
        w = 0.0  # If behaviour policy assigns zero probability, weight is zero
    end
    push!(weights, w)
end

# Truncate extreme weights
weights = min.(weights, 10.0)  # Cap at 10

# Estimate policy value: E[Y^π] = E[w * Y]
policy_value = sum(weights .* Y) / sum(weights)

# Compare to naive estimate (ignoring policy mismatch)
naive_value = mean(Y)

println("Off-policy evaluation:")
println("  Behavior policy π₀: Treat if L > 0.5")
println("  Target policy π: Treat if L > 0.3")
println("  Naive estimate (ignores policy): ", round(naive_value, digits=3))
println("  IS-weighted estimate: ", round(policy_value, digits=3))
Off-policy evaluation:
  Behavior policy π₀: Treat if L > 0.5
  Target policy π: Treat if L > 0.3
  Naive estimate (ignores policy): 0.109
  IS-weighted estimate: 0.158
Figure 28.2: Off-policy evaluation: evaluating a new policy using data from a different policy

28.6.5 Advantages of Causal Approach

CDMs enable:

  • Interventional simulation: Simulate under any policy
  • No positivity issues: Don’t need policy overlap (can simulate policies not observed)
  • Uncertainty quantification: Full uncertainty propagation
  • Mechanistic understanding: Understand why policies work, not just that they work

28.7 Continuous modified treatment policies (LMTP)

When the exposure is continuous, a scientifically useful intervention is often a modified treatment policy (a shift of the natural value \(A\) (e.g. add one standardised unit, subject to quantile clamps)) rather than \(do(A=a^*)\) for an arbitrary \(a^*\) (Díaz Muñoz and Laan 2012; Díaz et al. 2023). Estimation of δ-indexed total-effect curves (cross-fitted Super Learner nuisances, optional TMLE fluctuation) is illustrated with CausalTargeted; identification certificates come from the graph layer (Section 25.4).

NoteFrom CausalTargeted.jl

The scientific object in continuous MTP estimation is the shift policy, not a point \(do(A=a^*)\). Additive z-shifts are the default; multiplicative and threshold policies share the same ShiftPolicy type.

# packages/CausalTargeted.jl/src/shift_policies.jl (excerpt)
additive_shift_policy(; scale = "z", lower_q = 0.01, upper_q = 0.99) =
    ShiftPolicy(scale, lower_q, upper_q)
# Realised map (additive): A ↦ clamp(A + δ · s, L, U)

multiplicative_shift_policy(; lower_q = 0.01, upper_q = 0.99) =
    ShiftPolicy("multiplicative", lower_q, upper_q)
# A ↦ clamp(A * (1 + δ), L, U)

threshold_shift_policy(; lower_q = 0.01, upper_q = 0.99) =
    ShiftPolicy("threshold", lower_q, upper_q)
# Raise low exposures toward the median scale (clamped)

run_lmtp_grid estimates a δ-indexed TE curve under such a policy (with cross-fitted nuisances). Identification of the adjustment set remains upstream in CausalDynamics.

TipSmall-n checklist (CausalTargeted)

For continuous MTP and mediation with tens to low hundreds of units (Díaz et al. 2023; Laan and Rose 2011):

  1. Start from recommend_run_options(n; engine, n_mediators) (lean Super Learner when \(n < 80\); parallel=false by default).
  2. Inspect positivity / support (positivity_report or grid positivity=true) (Petersen et al. 2012).
  3. For mediation, sweep nested MC (mediation_n_mc_sweep) until signs and SEs stabilise (Liu et al. 2024).
  4. Report tipping-point / partial-\(R^2\) sensitivity (sensitivity_report) (Cinelli and Hazlett 2020).
  5. Treat discovery graphs as sensitivity only (merge_discovery_sensitivity!), never as silent DAG replacement.

Package prose and DOIs: CausalTargeted methods · references. Worked δ-grid: Chapter 23 (Policy Evaluation).

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using CausalTargeted DataFrames StableRNGs

rng = StableRNG(23)
df, _ = simulate_linear_mtp(120; rng = rng)
opts = recommend_run_options(nrow(df); engine = :lmtp)
println("Recommended learners (lean default): ", opts.learners_outcome)
println("Package DEFAULT_SL_LEARNERS: ", DEFAULT_SL_LEARNERS)

lmtp_grid = run_lmtp_grid(
    df, :A, :Y;
    baseline = [:W],
    deltas = [-1.0, -0.5, 0.0, 0.5, 1.0],
    folds = opts.folds,
    learners_outcome = opts.learners_outcome,
    parallel = false,
    positivity = true,
)
pos = positivity_report(df, :A; deltas = [-1.0, 0.0, 1.0])
sens = sensitivity_report(lmtp_grid.est[end], lmtp_grid.se[end]; n = nrow(df))
println("LMTP TE at δ = ", lmtp_grid.delta[end], ": ",
    round(lmtp_grid.est[end]; digits = 3), " (SE ",
    round(lmtp_grid.se[end]; digits = 3), ")")
println("Positivity rows: ", nrow(pos), "; sensitivity kinds: ", unique(sens.kind))
Recommended learners (lean default): (:glm, :mean)
Package DEFAULT_SL_LEARNERS: (:glm, :mean)
LMTP TE at δ = 1.0: 0.438 (SE 0.038)
Positivity rows: 3; sensitivity kinds: ["tipping_point", "partial_r2"]
NoteLean vs rich SuperLearner libraries

Production grids default to DEFAULT_SL_LEARNERS = (:glm, :mean) so Quarto renders stay fast and dependencies stay lean. For recovery / larger \(n\), pass RICH_SL_LEARNERS (GLM expansions, :glmnet* via MLJLinearModels, :evotree* when EvoTrees is loaded). Activate MLJ candidates with using MLJ, MLJLinearModels before the grid.

Parallel LMTP TE at δ = 0.5: 0.216
RICH_SL_LEARNERS (opt-in): (:glm, :glm_interact, :glm_quad, :glmnet, :glmnet_lasso, :glmnet_ridge, :evotree, :evotree_deep, :mean)
Figure 28.3: Continuous MTP total-effect curve under additive z-shifts (CausalTargeted run_lmtp_grid). Vertical bars are ±1.96 SE. (Díaz et al. 2023)

Figure 28.3 shows a small synthetic MTP curve with positivity and Cinelli–Hazlett-style sensitivity helpers attached (Cinelli and Hazlett 2020; Petersen et al. 2012). For multi-time exposures, use SequentialPolicy / run_sequential_lmtp with a CausalDynamics TemporalEffectQuery certificate.

28.8 Stratum context

This chapter addresses Doing in the Observable stratum: how can we evaluate policies? Policy evaluation uses observable data to reason about what will happen under different treatment strategies, bridging the Observable stratum (what we observe) with the Structural stratum (what would happen under interventions).

28.9 Key Takeaways

  1. Dynamic treatment strategies: Policies map history to treatment
  2. Policy evaluation: Estimate expected outcomes under policies
  3. Methods: G-computation, IPTW, TMLE, and continuous LMTP (CausalTargeted) all apply to policy evaluation
  4. Comparing policies: Choose best policy for given objective; check positivity before interpreting MTP curves

28.10 Further Reading