Dose-Response Estimation with HAL: Simulated and Real Data

Data Download

1. Introduction

This document demonstrates how to use the HAL-based plugin estimator to estimate the marginal causal dose-response curve in both simulated and real data settings. The estimator uses Highly Adaptive Lasso (HAL) with optional undersmoothing to flexibly model the conditional outcome, enabling valid inference on marginal intervention contrasts.

We provide two simulated examples and one real data example: - A binary outcome setting - A continuous outcome setting - A real-data application using RHC study data

---

2. Simulated Data Examples

In this section, we simulate hypothetical datasets to validate the HAL-based dose-response estimator. These examples serve to: - illustrate the mechanics of generating potential outcomes under continuous interventions, - evaluate estimator performance under known data-generating mechanisms, - and visualize the agreement between estimated and true dose-response curves.

2.1 Binary Outcome Example

This example simulates a binary outcome \(Y\), influenced by a continuous treatment \(A\) and a baseline confounder \(W\). The goal is to estimate the causal effect curve \(\mathbb{E}[Y | do(A = a)]\) using the plugin estimator with HAL under a binomial outcome model.

We simulate data where: - \(W\) is a baseline confounder - \(A\) is a continuous treatment - \(Y\) is a binary outcome depending on both \(A\) and \(W\)

2.1.1 Generate Simulated Data

set.seed(145)

generate_data <- function(n, a = NA){
  U_W <- rnorm(n)
  U_A <- rnorm(n, sd = 2)
  U_Y <- runif(n)
  W <- U_W

  if (is.na(a)) {
    A <- 2 - 0.5 * W + U_A
    A <- pmin(pmax(A, 0), 5)
  } else {
    A <- rep(a, n)
  }

  Y <- as.numeric(U_Y < plogis(-3 + 0.5*W + 1.25*A - 0.5*W*A))
  data.frame(W, A, Y)
}

obs <- generate_data(n = 500)
head(obs)

##           W         A Y
## 1 0.6869129 3.7840125 0
## 2 1.0663631 3.2232627 1
## 3 0.5367006 0.5392545 1
## 4 1.9060287 3.4649305 0
## 5 1.0631596 0.0000000 0
## 6 1.3703436 0.8576644 0

2.1.2 Estimate Dose-Response Curve

DRC_fit_UAdaptive <- fit_UHAL_DRC(dat = obs, y_var_name = "Y", trt_var_name = "A", family = "binomial")

2.1.3 Visualize Estimated Curve

df_ests <- DRC_fit_UAdaptive$curve_est

ggplot(df_ests, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  labs(title = "Estimated Dose-Response (Binary Outcome)",
       x = "Treatment (A)",
       y = "E[Y | do(A = a)]") +
  theme_bw()

2.1.4 Compare to True Curve

a.vals <- seq(0, 5, length.out = 20)
psi0_a <- sapply(a.vals, function(a) mean(generate_data(n = 1e6, a = a)$Y))
psi0 <- data.frame(a = a.vals, psi0 = psi0_a)

ggplot(df_ests, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  geom_point(data = psi0, aes(x = a, y = psi0), color = "red") +
  geom_line(data = psi0, aes(x = a, y = psi0), color = "red") +
  labs(title = "Estimated vs. True Dose-Response (Binary Y)",
       y = "E[Y | do(A = a)]") +
  theme_bw()

---

2.2 Continuous Outcome Example

This simulation example uses a continuous outcome \(Y\), which is a nonlinear function of treatment \(A\) and baseline covariate \(W\). This showcases how HAL handles nonlinear functional relationships and continuous outcomes. We again compare the estimated dose-response curve with the true curve obtained by Monte Carlo integration.

2.2.1 Simulate Continuous Outcome Data

generate_data_2 <- function(n, a = NA) {
  U_W <- rnorm(n)
  U_A <- rnorm(n, sd = 0.8)
  U_Y <- rnorm(n, sd = 0.5)

  W <- U_W
  A <- if (is.na(a)) 0.5 * W + U_A else rep(a, n)
  logit_Y <- 0.5 * A^2 + A + 0.3 * W + U_Y
  Y <- 1 / (1 + exp(-logit_Y))
  data.frame(W, A, Y)
}

obs <- generate_data_2(n = 500)

2.2.2 Visualize Data

par(mfrow = c(1, 2))
plot(obs$W, obs$A, main = "W vs A", pch = 19, col = rgb(0, 0, 1, 0.4))
plot(obs$A, obs$Y, main = "A vs Y", pch = 19, col = rgb(1, 0, 0, 0.4))

2.2.3 Estimate and Plot Dose-Response Curve

DRC_fit_UAdaptive <- fit_UHAL_DRC(obs, "Y", "A", "gaussian")
df_ests <- DRC_fit_UAdaptive$curve_est

ggplot(df_ests, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  labs(title = "Estimated Dose-Response (Continuous Outcome)",
       x = "Treatment (A)",
       y = "E[Y | do(A = a)]") +
  theme_bw()

2.2.4 Compare to True Curve

a.vals <- seq(min(obs$A), max(obs$A), length.out = 100)
psi0 <- data.frame(
  a = a.vals,
  psi0 = sapply(a.vals, function(a) mean(generate_data_2(n = 1e6, a = a)$Y))
)

ggplot(df_ests, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  geom_line(data = psi0, aes(x = a, y = psi0), color = "red") +
  labs(title = "Estimated vs. True Dose-Response (Continuous Y)",
       y = "E[Y | do(A = a)]") +
  theme_bw()

---

3. Real Data Illustration: Right Heart Catheterization (RHC) Study

In this section, we demonstrate the use of HAL-based marginal dose-response estimation on observational data from the RHC study. This dataset contains detailed clinical information on critically ill patients. We investigate the causal effect of age on two outcomes: in-hospital mortality and length of hospital stay, adjusting for demographic and clinical covariates such as sex and severity of illness.

3.1 Common Setup

We begin by loading the dataset, processing variables, and defining a reusable function for fitting and visualizing dose-response curves.

# Load packages
library(dplyr)
library(ggplot2)
library(halDRC)

# Load RHC dataset
ObsData <- read.csv("https://hbiostat.org/data/repo/rhc.csv", header = TRUE)

# Clean outcome: length of stay (Y)
ObsData$Y_los <- ObsData$dschdte - ObsData$sadmdte
ObsData$Y_los[is.na(ObsData$Y_los)] <- ObsData$dthdte[is.na(ObsData$Y_los)] - ObsData$sadmdte[is.na(ObsData$Y_los)]

# Clean outcome: in-hospital death (binary)
ObsData$Y_death <- as.integer(ObsData$dthdte == ObsData$dschdte)
ObsData$Y_death[is.na(ObsData$Y_death)] <- 0

# Rename common confounders
ObsData <- ObsData %>%
  rename(
    APACHE.score = aps1,
    DASIndex = das2d3pc,
    Creatinine = crea1,
    Bilirubin = bili1,
    WBC = wblc1,
    Heart.rate = hrt1
  ) %>%
  mutate(
    sex_m = ifelse(sex == 'Male', 1, 0),
    ProcedureCount = round(runif(n(), min = 0, max = 10))
  )


# ℹ️ For computational efficiency in this demo, we randomly selected 500 observations from the full cleaned dataset for this HAL analysis.
set.seed(2025)
ObsData <- ObsData %>% sample_n(500)

3.2 Causal Relationship 1: Effect of Age on Length of Stay

We treat age as the exposure of interest and estimate its causal effect on hospital length of stay, controlling for baseline severity and other confounders.

drc_data1 <- ObsData %>%
  select(
    Y = Y_los,
    A = age,
    sex_m, APACHE.score, DASIndex, Creatinine, Bilirubin, WBC, Heart.rate, Y_death
  ) %>%
  na.omit()

drc_fit1 <- fit_UHAL_DRC(dat = drc_data1, y_var_name = "Y", trt_var_name = "A", family = "gaussian")
drc_est1 <- drc_fit1$curve_est

ggplot(drc_est1, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  labs(title = "Dose-Response: Age vs Length of Stay", 
       x = "Age (A)", 
       y = "Expected Length of Stay (Y)") +
  theme_bw()

3.2 Causal Relationship 2: Effect of Age on In-Hospital Mortality

Here, we again consider age as the treatment variable but now estimate its effect on the probability of in-hospital death, while adjusting for mortality and clinical severity indicators.

drc_data2 <- ObsData %>%
  select(
    Y = Y_death,
    A = age,
    sex_m, APACHE.score, DASIndex, Creatinine, Bilirubin, WBC, Heart.rate, Y_los
  ) %>%
  na.omit()

drc_fit2 <- fit_UHAL_DRC(dat = drc_data2, y_var_name = "Y", trt_var_name = "A", family = "binomial")
drc_est2 <- drc_fit2$curve_est

ggplot(drc_est2, aes(x = a, y = y_hat)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.3) +
  labs(title = "Dose-Response: Age vs In-Hospital Mortality", 
       x = "Age (A)", 
         y = "Expected In-Hospital Death (Y)") +
  theme_bw()

4. Summary

This unified demonstration shows how to apply the HAL-based plugin estimator to estimate causal dose-response curves for both simulated and real-world observational data. The method enables flexible functional modeling and valid inference, with clear application to diverse outcomes and settings.