Dose-Response Estimation for Continuous Treatments

Estimates the dose-response function (average potential outcome as a function of treatment dose) for continuous treatments using a generalized propensity score (GPS) approach or regression-based methods. Produces smooth dose-response curves with confidence bands.

Usage

dose_response(
  data,
  outcome,
  treatment,
  covariates = NULL,
  method = c("gps", "regression", "loess"),
  n_grid = 50,
  poly_degree = 2,
  conf_level = 0.95,
  boot_reps = 200,
  seed = 42,
  trim_percentile = 0.01,
  plot = TRUE
)

Arguments

data

A data frame.

outcome

Character. Outcome variable name.

treatment

Character. Continuous treatment variable name.

covariates

Character vector. Covariate names for GPS model.

method

Character. Estimation method: "gps" (generalized propensity score, default), "regression" (polynomial OLS), or "loess" (nonparametric LOESS).

n_grid

Integer. Number of grid points for the dose-response curve. Default 50.

poly_degree

Integer. Polynomial degree for GPS or regression methods. Default 2.

conf_level

Numeric. Confidence level. Default 0.95.

boot_reps

Integer. Bootstrap replicates for CI. Default 200.

seed

Integer. Random seed. Default 42.

trim_percentile

Numeric. Trim extreme treatment values at this percentile (both tails). Default 0.01.

plot

Logical. Return a dose-response plot. Default TRUE.

Value

A list with:

dose_response

Data frame: dose grid, estimated outcome, SE, CI bounds.

ape

Numeric. Average partial effect (slope at mean dose).

plot

ggplot2 dose-response curve.

References

Hirano, K., & Imbens, G. W. (2004). The propensity score with continuous treatments. In A. Gelman & X.-L. Meng (Eds.), Applied Bayesian Modeling and Causal Inference from Incomplete- Data Perspectives (pp. 73–84). Wiley.

Examples

set.seed(42)
n <- 500
df <- data.frame(
  X1 = rnorm(n), X2 = rnorm(n),
  D  = runif(n, 0, 10)
)
df$D  <- df$D + df$X1  # confounded dose
df$Y  <- 0.5 * df$D - 0.03 * df$D^2 + df$X1 + rnorm(n)

result <- dose_response(
  data       = df,
  outcome    = "Y",
  treatment  = "D",
  covariates = c("X1", "X2"),
  method     = "gps"
)
result$plot