Regression Adjustment

The idea

The simplest way to estimate a causal effect from observational data: run a regression of the outcome on the treatment indicator and the confounders.

\[Y_i = \alpha + \tau D_i + \beta X_i + \varepsilon_i\]

The coefficient \(\hat{\tau}\) is your treatment effect estimate, “adjusted” for \(X\). OLS compares treated and control units within levels of X — holding the confounder constant to isolate the effect of treatment.

The gap between the two regression lines (treated vs control) is \(\hat{\tau}\).

Assumptions

1. Selection on observables (CIA): all confounders are included in \(X\). If an unobserved variable drives both treatment and outcome, the estimate is biased — same as every other tool under SOO.
2. Correct functional form: the relationship between \(X\) and \(Y\) is correctly specified. If you fit a linear model but the truth is quadratic, the adjustment is wrong and \(\hat{\tau}\) is biased. This is the key vulnerability.
3. Overlap: treated and control groups share common support in \(X\).

The vulnerability: functional form

Regression adjustment models the outcome — it writes down a specific equation for how \(Y\) relates to \(X\) and \(D\). That equation is the functional form:

• Linear: \(Y = \alpha + \tau D + \beta X\) — assumes \(Y\) changes at a constant rate with \(X\). A straight line.
• Quadratic: \(Y = \alpha + \tau D + \beta_1 X + \beta_2 X^2\) — allows curvature.
• Log: \(Y = \alpha + \tau D + \beta \log(X)\) — assumes diminishing returns.

When you run lm(y ~ treat + x), you’ve chosen the linear form. If the truth is quadratic but you fit a line, the line can’t capture the curvature — and that misfit contaminates your estimate of \(\hat{\tau}\). The model is wrong, so the estimate is biased, even when all confounders are observed. This is an estimation bias, not an identification failure.

Compare this to IPW, which models the treatment assignment instead — it writes down an equation for \(P(D = 1 \mid X)\), not for \(Y\). IPW never touches the outcome. It just reweights the raw \(Y\) values using the propensity score. So if the outcome-\(X\) relationship is nonlinear, IPW doesn’t care — as long as the propensity score model is right.

	What it models	What can go wrong
Regression adjustment	\(E[Y \mid X, D]\) (the outcome)	Wrong functional form for \(Y\)
IPW	\(P(D=1 \mid X)\) (the treatment)	Wrong propensity score model
Doubly robust	Both	Fails only if both are wrong

Each method bets on getting one model right. Doubly robust estimation hedges by combining both — it only needs one to be correct.

#| standalone: true
#| viewerHeight: 620

library(shiny)

ui <- fluidPage(
  tags$head(tags$style(HTML("
    .stats-box {
      background: #f0f4f8; border-radius: 6px; padding: 14px;
      margin-top: 12px; font-size: 14px; line-height: 1.9;
    }
    .stats-box b { color: #2c3e50; }
    .good { color: #27ae60; font-weight: bold; }
    .bad  { color: #e74c3c; font-weight: bold; }
  "))),

  sidebarLayout(
    sidebarPanel(
      width = 3,

      sliderInput("n", "Sample size:",
                  min = 200, max = 2000, value = 500, step = 100),

      sliderInput("ate", "True ATE:",
                  min = 0, max = 5, value = 2, step = 0.5),

      sliderInput("confounding", "Confounding strength:",
                  min = 0, max = 3, value = 1.5, step = 0.25),

      selectInput("shape", "True outcome model:",
                  choices = c("Linear (Y ~ X)"      = "linear",
                              "Quadratic (Y ~ X\u00b2)" = "quadratic")),

      actionButton("go", "New draw",
                   class = "btn-primary", width = "100%"),

      uiOutput("results")
    ),

    mainPanel(
      width = 9,
      fluidRow(
        column(6, plotOutput("scatter_plot", height = "400px")),
        column(6, plotOutput("compare_plot", height = "400px"))
      )
    )
  )
)

server <- function(input, output, session) {

  dat <- reactive({
    input$go
    n    <- input$n
    ate  <- input$ate
    conf <- input$confounding
    is_linear <- input$shape == "linear"

    # Confounder
    x <- rnorm(n)

    # Treatment depends on x
    p_true <- pnorm(conf * x)
    treat  <- rbinom(n, 1, p_true)

    # Outcome depends on x (linearly or quadratically) + treatment
    if (is_linear) {
      y <- 1 + 2 * x + ate * treat + rnorm(n)
    } else {
      y <- 1 + 2 * x^2 + ate * treat + rnorm(n)
    }

    # Naive (no controls)
    naive <- mean(y[treat == 1]) - mean(y[treat == 0])

    # Regression adjustment (always linear: Y ~ D + X)
    fit <- lm(y ~ treat + x)
    reg_est <- coef(fit)["treat"]

    # IPW for comparison
    ps <- fitted(glm(treat ~ x, family = binomial))
    w <- ifelse(treat == 1, 1 / ps, 1 / (1 - ps))
    ipw_est <- weighted.mean(y[treat == 1], w[treat == 1]) -
               weighted.mean(y[treat == 0], w[treat == 0])

    list(x = x, treat = treat, y = y, fit = fit,
         naive = naive, reg_est = reg_est, ipw_est = ipw_est,
         ate = ate, is_linear = is_linear)
  })

  output$scatter_plot <- renderPlot({
    d <- dat()
    par(mar = c(4.5, 4.5, 3, 1))

    plot(d$x, d$y, pch = 16, cex = 0.5,
         col = ifelse(d$treat == 1, "#3498db50", "#e74c3c50"),
         xlab = "X (confounder)", ylab = "Y (outcome)",
         main = "Data + Regression Fit")

    # Fitted regression lines (always linear)
    x_seq <- seq(min(d$x), max(d$x), length.out = 200)
    pred_t <- predict(d$fit,
      newdata = data.frame(treat = 1, x = x_seq))
    pred_c <- predict(d$fit,
      newdata = data.frame(treat = 0, x = x_seq))
    lines(x_seq, pred_t, col = "#3498db", lwd = 2.5)
    lines(x_seq, pred_c, col = "#e74c3c", lwd = 2.5)

    # True conditional means (dashed)
    if (d$is_linear) {
      true_t <- 1 + 2 * x_seq + d$ate
      true_c <- 1 + 2 * x_seq
    } else {
      true_t <- 1 + 2 * x_seq^2 + d$ate
      true_c <- 1 + 2 * x_seq^2
    }
    lines(x_seq, true_t, col = "#3498db", lwd = 1.5, lty = 3)
    lines(x_seq, true_c, col = "#e74c3c", lwd = 1.5, lty = 3)

    # Gap arrow at x = 0
    mid_t <- predict(d$fit,
      newdata = data.frame(treat = 1, x = 0))
    mid_c <- predict(d$fit,
      newdata = data.frame(treat = 0, x = 0))
    arrows(0, mid_c, 0, mid_t, code = 3,
           col = "#2c3e50", lwd = 2, length = 0.08)
    text(0.3, (mid_t + mid_c) / 2,
         bquote(hat(tau) == .(round(coef(d$fit)["treat"], 2))),
         cex = 0.95, font = 2)

    legend("topleft", bty = "n", cex = 0.8,
           legend = c("Treated", "Control",
                      "Regression fit", "True E[Y|X,D]"),
           col = c("#3498db", "#e74c3c", "gray40", "gray40"),
           pch = c(16, 16, NA, NA),
           lwd = c(NA, NA, 2.5, 1.5),
           lty = c(NA, NA, 1, 3))
  })

  output$compare_plot <- renderPlot({
    d <- dat()
    par(mar = c(4.5, 9, 3, 2))

    ests <- c(d$ipw_est, d$reg_est, d$naive)
    labels <- c("IPW", "Regression adj.", "Naive")
    cols <- c("#27ae60", "#3498db", "#e74c3c")

    xlim <- range(c(ests, d$ate))
    pad  <- max(diff(xlim) * 0.4, 0.5)
    xlim <- xlim + c(-pad, pad)

    plot(ests, 1:3, pch = 19, cex = 2, col = cols,
         xlim = xlim, ylim = c(0.5, 3.5),
         yaxt = "n", xlab = "Estimated treatment effect",
         ylab = "", main = "Estimator Comparison")
    axis(2, at = 1:3, labels = labels,
         las = 1, cex.axis = 0.9)

    abline(v = d$ate, lty = 2, col = "#2c3e50", lwd = 2)
    text(d$ate, 3.45,
         paste0("True ATE = ", d$ate),
         cex = 0.85, font = 2, col = "#2c3e50")

    # Lines from true to estimate
    segments(d$ate, 1:3, ests, 1:3,
             col = cols, lwd = 2, lty = 2)
  })

  output$results <- renderUI({
    d <- dat()
    b_naive <- d$naive - d$ate
    b_reg   <- d$reg_est - d$ate
    b_ipw   <- d$ipw_est - d$ate
    reg_ok  <- abs(b_reg) < abs(b_naive) * 0.5
    ipw_ok  <- abs(b_ipw) < abs(b_naive) * 0.5

    tags$div(class = "stats-box",
      HTML(paste0(
        "<b>True ATE:</b> ", d$ate, "<br>",
        "<b>Naive:</b> <span class='bad'>",
        round(d$naive, 3), "</span>",
        " (bias: ", round(b_naive, 3), ")<br>",
        "<b>Reg. adj:</b> <span class='",
        ifelse(reg_ok, "good", "bad"), "'>",
        round(d$reg_est, 3), "</span>",
        " (bias: ", round(b_reg, 3), ")<br>",
        "<b>IPW:</b> <span class='",
        ifelse(ipw_ok, "good", "bad"), "'>",
        round(d$ipw_est, 3), "</span>",
        " (bias: ", round(b_ipw, 3), ")"
      ))
    )
  })
}

shinyApp(ui, server)

Things to try

• Linear outcome, confounding = 1.5: regression adjustment works well. The fitted lines (solid) match the true curves (dotted). The estimate is close to the true ATE.
• Switch to Quadratic: now the true outcome is \(Y = 1 + 2X^2 + \tau D\), but regression still fits a line. The solid lines miss the curvature (compare to the dotted true curves). The regression estimate is biased — it’s a functional form problem, not a confounding problem.
• Quadratic outcome, compare to IPW: IPW doesn’t model the outcome, so it still works. The green dot (IPW) is close to the true ATE; the blue dot (regression) is not.
• Confounding = 0: treatment is random. All three estimators agree, and functional form doesn’t matter (no confounding to adjust for).

The bottom line

Regression adjustment is the simplest and most common tool. It works great when the outcome model is correctly specified. But it’s fragile to functional form misspecification — fitting a line when the truth is curved introduces bias even when all confounders are observed.

This is why alternatives like IPW and entropy balancing exist: they avoid modeling the outcome entirely. And doubly robust methods hedge both bets.

---

In Stata

* Basic regression adjustment
reg outcome treatment x1 x2

* Stata's teffects version (potential outcomes framework)
teffects ra (outcome x1 x2) (treatment)

* Lin (2013) fix: interact treatment with demeaned covariates
* (robust to functional form misspecification)
foreach var of varlist x1 x2 {
    sum `var', meanonly
    gen `var'_dm = `var' - r(mean)
}
reg outcome treatment c.treatment#c.(x1_dm x2_dm) x1 x2

Plain reg outcome treatment x1 x2 and teffects ra give the same point estimate under linearity. The Lin (2013) interaction approach is safer when you’re unsure about functional form.

---

Did you know?

• Regression adjustment is so natural that most applied researchers use it without thinking of it as a “method.” But Freedman (2008) showed that even in randomized experiments, regression adjustment can be biased in finite samples if the model is wrong — the linear adjustment can introduce bias that wasn’t there. Lin (2013) showed how to fix this: interact the treatment indicator with all covariates, and the resulting estimator is asymptotically unbiased regardless of the true functional form.

• The vulnerability of regression to functional form is a key motivation for semiparametric and nonparametric methods. Instead of assuming \(E[Y \mid X] = \alpha + \beta X\) (which might be wrong), methods like kernel regression or series estimation let the data determine the shape.