Standard Iv
Executive Summary

Executive Summary

Ordinary versus instrumented effect of Years Of Schooling on Hourly Earnings

Observations
1200
OLS Estimate
2.962
IV Estimate
1.739
Confounding Gap
-1.223
First-Stage F
154.9
Instrument Strength
strong enough
Across 1,200 rows, an ordinary least squares regression puts the effect of Years Of Schooling on Hourly Earnings at +2.96 per unit (95% CI 2.83 to 3.09). Instrumenting Years Of Schooling with Miles To College moves that to +1.74 (95% CI 1.3 to 2.18, p < 0.001), a difference of -1.22. That difference is the story: it is the size of what the confounding was adding to the ordinary estimate, and the endogeneity test confirms the two estimates differ by more than sampling noise. The instrument's first-stage F is 155, at or above the conventional threshold of 10. With exactly one instrument the exclusion restriction is untestable — no statistic in this analysis can check whether Miles To College reaches Hourly Earnings by some route other than Years Of Schooling. This rests on an assumption the data cannot check: that Miles To College affects Hourly Earnings only by changing Years Of Schooling, and is otherwise unrelated to whatever else moves Hourly Earnings. With exactly one instrument that assumption is untestable — no statistic in this analysis can check it. Nothing here establishes causation on its own. The instrumented estimate describes the subgroup whose Years Of Schooling actually moved with Miles To College, not everyone in the data.
Suggested Interpretation

Accounting for confounding via instrumental variables, the causal return to schooling is 1.74 dollars per hour per additional year (95% CI 1.30 to 2.18, p < 0.001), compared to 2.96 from ordinary regression—a confounding gap of −1.22. The first-stage F of 155 clears the threshold of 10, making the instrumented estimate usable. However, with exactly one instrument the exclusion restriction—that Miles To College reaches Hourly Earnings only through Years Of Schooling—is untestable and must be assumed. The estimate applies only to the subgroup whose schooling actually responded to proximity to college.

Overview

Analysis Overview

Two-stage least squares of Hourly Earnings on Years Of Schooling, instrumented by Miles To College, across 1,200 observations.

N Observations1200
N Instruments1
N Controls0
First Stage F154.9
Suggested Interpretation

Ordinary regression conflates Years Of Schooling with unmeasured traits that raise both schooling and earnings. This analysis isolates the causal effect using Miles To College as an instrument: the first stage explains 11.4 percent of Years Of Schooling variation through the instrument, and the second stage uses only that instrumented portion to estimate the effect on Hourly Earnings. The instrumented estimate answers the causal question because anything unmeasured that the instrument is unrelated to drops out. Standard errors are rebuilt from residuals against actual Years Of Schooling, not the second stage's predicted values, correcting for two-stage least squares variance.

Data Preparation

Data Quality

Row accounting, column typing, and which mapped columns were excluded.

Initial Rows1200
Final Rows1200
Rows Removed0
N Instruments1
Suggested Interpretation

All 1,200 rows were complete on Hourly Earnings, Years Of Schooling, and Miles To College; no rows were removed. Numeric columns had missing values filled with the column median before completeness checks. The analysis uses 1 instrument (Miles To College) and no control columns. All mapped columns were usable for estimation.

Data Table

First Stage

First-stage coefficients and the F statistic on the excluded instruments.

TermRoleEstimateStd ErrorT StatP Value
Miles To CollegeInstrument-0.08770.007-12.44< 0.001
Suggested Interpretation

The first-stage F of 155 on 1 and 1,198 degrees of freedom (p < 0.001) clears the conventional threshold of 10, establishing that Miles To College is a strong instrument. Miles To College carries a coefficient of −0.0877 (t = −12.44, p < 0.001), explaining 11.4 percent of the remaining variation in Years Of Schooling beyond controls. A strong first stage confirms relevance—the instrument moves the treatment—but says nothing about validity. The exclusion restriction, which requires Miles To College to affect Hourly Earnings only through Years Of Schooling, cannot be tested from the first stage and remains an untestable assumption.

Visualization

Instrument versus Treatment

The first-stage relationship the instrumented estimate runs through.

Suggested Interpretation

The scatter of Miles To College against Years Of Schooling shows a visible downward slope: as distance to college increases, years of schooling tend to decrease. The pattern is moderately tight around the trend, consistent with the first-stage F of 155. This slope is the mechanism through which the instrumented estimate operates. A tight band indicates the instrument pins down a substantial share of treatment variation; a diffuse cloud would mean lower precision. However, this plot cannot distinguish a valid from an invalid instrument—a strong first stage and an excluded variable that violates the exclusion restriction look identical visually.

Visualization

Ordinary versus Instrumented

The confounded estimate next to the instrumented one, with 95 percent intervals.

Suggested Interpretation

The ordinary least squares estimate (+2.96, 95% CI 2.83 to 3.09) is larger than the instrumented estimate (+1.74, 95% CI 1.30 to 2.18), a gap of −1.22 representing the confounding removed by the instrument. The ordinary estimate credits Years Of Schooling with unmeasured traits that raise both schooling and earnings; the instrumented estimate uses only the variation in schooling driven by proximity to college. The wider instrumented confidence interval reflects the cost of using only 11.4 percent of schooling variation. Both intervals are clearly above zero, but the exclusion restriction—that Miles To College affects earnings only through schooling—cannot be verified and must be assumed.

Data Table

Effect Estimates

Ordinary and instrumented estimates with intervals and p-values.

Estimate TypeEstimateStd ErrorCI LowCI HighP ValueInterpretation
Ordinary least squares2.9620.06692.8313.093< 0.001Confounded: mixes the effect of Years Of Schooling with anything unmeasured that moves both it and Hourly Earnings.
Instrumental variables (2SLS)1.7390.22351.32.177< 0.001Uses only the variation in Years Of Schooling driven by the instrument; valid only if the exclusion restriction holds.
Suggested Interpretation

The short answer

The instrumented effect is 1.739 per year (standard error 0.2235, 95% CI 1.3 to 2.177, p < 0.001), compared to the ordinary estimate of 2.962 (standard error 0.0669). The instrumented standard error is wider because it uses only 11.4 percent of schooling variation, and it is rebuilt from residuals against actual schooling, not the second stage's predicted value.

The detail

Ordinary least squares: estimate 2.962, standard error 0.0669, 95% CI 2.831 to 3.093, p < 0.001. Instrumental variables: estimate 1.739, standard error 0.2235, 95% CI 1.3 to 2.177, p < 0.001. The instrumented standard error is the two-stage least squares standard error, measured against actual Years Of Schooling. The second stage's own standard error would have been 0.317, a factor of 0.705 away from the correct figure. Both estimates are highly significant.

What this can't tell you

This rests on the assumption that Miles To College affects Hourly Earnings only by changing Years Of Schooling and is otherwise unrelated to whatever else moves earnings. With exactly one instrument that assumption is untestable. The instrumented estimate describes only the subgroup whose schooling actually moved with distance to college.

Data Table

Assumptions & Diagnostics

What was tested, what was assumed, and what the estimate actually applies to.

CheckResultInterpretation
Instrument strength (first-stage F on the excluded instruments)F = 155 on 1 and 1,198 degrees of freedom, p < 0.001Above the conventional threshold of 10, so Miles To College moves Years Of Schooling strongly enough for the instrumented estimate to be read.
Instrument relevance (extra variation explained)11.4 percent of the variation in Years Of Schooling left by the controlsThe share of Years Of Schooling that the instrument explains beyond the control columns. This is the only part of Years Of Schooling the instrumented estimate uses.
Over-identification (Sargan test of the exclusion restriction)Not testable with one instrumentWith exactly one instrument and one confounded regressor the model is exactly identified, and the exclusion restriction is UNTESTABLE — no statistic in this or any other analysis can check it. It is assumed, not established.
Endogeneity (Wu-Hausman test of ordinary vs instrumented)t = 6.7, p < 0.001The two estimates differ by more than sampling noise, which is evidence that Years Of Schooling is confounded and the ordinary least squares estimate is biased.
Exclusion restrictionAssumed, never testedThe analysis assumes Miles To College affects Hourly Earnings ONLY by changing Years Of Schooling, and is unrelated to whatever else moves Hourly Earnings. No amount of data can verify that. It is a claim about the world that you make, and the whole estimate rests on it.
What the estimate applies toA local average treatment effectTwo-stage least squares recovers the effect of Years Of Schooling on Hourly Earnings for the subgroup whose Years Of Schooling actually responded to Miles To College — not the average effect across everyone in the data. If the effect differs between people who respond to the instrument and people who do not, this number does not describe the second group. It also assumes the instrument pushes every unit in the same direction.
Suggested Interpretation

Relevance is testable and confirmed: the first-stage F of 155 exceeds the threshold of 10. Over-identification is untestable because the model is exactly identified with one instrument and one confounded regressor; the exclusion restriction cannot be checked statistically and must be assumed. Endogeneity is testable: the Wu-Hausman test (t = 6.7, p < 0.001) shows the ordinary and instrumented estimates differ significantly, confirming that Years Of Schooling is confounded. The estimate applies to a local average treatment effect—the subgroup whose schooling responded to Miles To College—and assumes the instrument pushes all units in the same direction. No statistic can verify either the exclusion restriction or monotonicity; both are claims about the world you are making.

Methodology

Methodology

Statistical methodology and diagnostics for Instrumental Variables

Statistical Method

Instrumental Variables

Standard-library analysis: what is the causal effect when the treatment is confounded? Map an outcome, a confounded treatment, and at least one instrument — a variable that shifts the treatment but has no other route to the outcome — and get two-stage least squares built stage by stage: the first-stage F on the excluded instruments and the instrument-strength verdict it implies, the ordinary least squares estimate and the instrumented estimate side by side with the difference between them made explicit, an over-identification (Sargan) test when more than one instrument is supplied, an endogeneity (Wu-Hausman) test, and a plain statement of the exclusion restriction the whole estimate rests on.

Data
N = 1200 observations
Assumptions
  • Relevance: the instrument moves the treatment (this IS tested — the first-stage F on the excluded instruments)
  • Exclusion restriction: the instrument affects the outcome ONLY through the treatment (this is NOT testable and is assumed)
  • Independence: the instrument is unrelated to whatever else moves the outcome (not testable)
  • Monotonicity: the instrument pushes every unit's treatment in the same direction (not testable)
  • The relationship between treatment and outcome is linear in the treatment
  • Observations are independent; the standard errors are the classical homoskedastic two-stage least squares standard errors
Limitations
  • The exclusion restriction cannot be verified by any statistic in this or any other analysis — the estimate is only as good as that assumption
  • With exactly one instrument the model is exactly identified and over-identification is untestable, full stop
  • The Sargan test only detects DISAGREEMENT between instruments; instruments that are all invalid in the same direction pass it
  • Two-stage least squares recovers a local average treatment effect — the effect for the subgroup whose treatment responded to the instrument, not the average effect across everyone
Software & Citation
MCP Analytics · mcpanalytics.ai
Code Appendix

Analysis Code

Complete R source code for this analysis

Instrumental Variables — Two-Stage Least Squares

Estimates the effect of a treatment on an outcome when the treatment is confounded: something unmeasured moves both, so an ordinary least squares comparison is biased. An instrument — a variable that shifts the treatment but has no other path to the outcome — is used to isolate the part of the treatment's variation that is unrelated to the confounder.

Why This Method?

Two-stage least squares regresses the treatment on the instrument (plus any controls) to build a predicted treatment, then regresses the outcome on that predicted treatment. The ordinary estimate and the instrumented estimate are reported side by side: the distance between them is the size of the confounding the instrument removed.

What This Analysis Covers

  • The first stage: does the instrument actually move the treatment?

The F statistic on the excluded instruments and the verdict it implies.

  • The ordinary least squares estimate and the instrumented estimate, with

the difference between them made explicit.

  • Over-identification (Sargan) when more than one instrument is supplied —

and an explicit statement that with exactly one instrument the exclusion restriction cannot be tested at all.

  • An endogeneity (Wu-Hausman) test of whether the two estimates differ by

more than sampling noise.

Standard Library

Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {outcome, treatment, instrument_1..N, covariate_1..M}. All narrative is derived from the user's own column names and computed values. The exclusion restriction is an ASSUMPTION the data cannot check, and every summary says so.

Implementation Note

Both stages are fitted with base lm. The second stage's own standard errors are WRONG for two-stage least squares, because lm computes them from residuals taken against the FITTED treatment. The correct two-stage least squares covariance rebuilds the residual variance from the STRUCTURAL equation — the outcome minus the fitted coefficients applied to the ACTUAL treatment — and then scales the second-stage cross-product inverse by it. That is done explicitly below.

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

Step 1: Semantic column discovery

initial_rows <- nrow(df)
  if (!"outcome" %in% names(df)) {
    stop("column_mapping must map an &#x27;outcome' column (the numeric result you want to explain)")
  }
  if (!"treatment" %in% names(df)) {
    stop("column_mapping must map a &#x27;treatment' column (the confounded cause whose effect you want)")
  }
  inst_cols <- grep("^instrument_[0-9]+$", names(df), value = TRUE)
  inst_cols <- inst_cols[order(as.integer(sub("^instrument_", "", inst_cols)))]
  if (length(inst_cols) == 0) {
    stop("column_mapping must map at least one instrument column(instrument_1) — a variable that shifts the treatment but has no other route to the outcome")
  }
  cov_cols <- grep("^covariate_[0-9]+$", names(df), value = TRUE)
  cov_cols <- cov_cols[order(as.integer(sub("^covariate_", "", cov_cols)))]

  outcome_name   <- humanize_semantic("outcome", col_map)
  treatment_name <- humanize_semantic("treatment", col_map)
  inst_names <- setNames(humanize_semantic(inst_cols, col_map), inst_cols)
  cov_names  <- if (length(cov_cols) > 0) {
    setNames(humanize_semantic(cov_cols, col_map), cov_cols)
  } else character(0)

Step 2: Outcome and treatment must be numeric (95% coercion rule)

coerce_required <- function(v, label, role_hint) {
    if (is.numeric(v)) return(v)
    conv <- suppressWarnings(as.numeric(as.character(v)))
    n_orig <- sum(!is.na(v) & trimws(as.character(v)) != "")
    if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) return(conv)
    stop(sprintf(
      "The %s column(&#x27;%s') must be numeric — %s.", role_hint, label,
      if (role_hint == "outcome")
        "the result you want to explain, such as earnings, spend, or a score"
      else
        "the cause whose effect you want, such as years of schooling, dose, or price"))
  }
  df$outcome   <- coerce_required(df$outcome, outcome_name, "outcome")
  df$treatment <- coerce_required(df$treatment, treatment_name, "treatment")

Step 3: Type every instrument and covariate.

Numeric when at least 95% of non-blank values convert; otherwise a factor with blanks as "Missing" and rare levels lumped into "Other". Near-unique text columns are identifiers, not variables — excluded.

prep_var <- function(cc, dropped) {
    x <- df[[cc]]
    if (!is.numeric(x)) {
      conv <- suppressWarnings(as.numeric(as.character(x)))
      n_orig <- sum(!is.na(x) & trimws(as.character(x)) != "")
      if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) df[[cc]] <<- conv
    }
    x <- df[[cc]]
    if (is.numeric(x)) {
      med <- median(x, na.rm = TRUE)
      if (is.na(med)) return(c(dropped, cc))
      x[is.na(x)] <- med
      df[[cc]] <<- x
      if (is.na(var(x)) || isTRUE(var(x) == 0)) return(c(dropped, cc))
    } else {
      x <- as.character(x)
      x[is.na(x) | trimws(x) == ""] <- "Missing"

Decide identifier-vs-category on the RAW level count, BEFORE lumping — lumping a 1,200-value ID column into 13 levels would otherwise hide it.

n_levels_raw <- length(unique(x))
      if (n_levels_raw > nrow(df) / 2 || n_levels_raw <= 1) {
        return(c(dropped, cc))
      }
      tab <- sort(table(x), decreasing = TRUE)
      if (length(tab) > 12) {
        keep_lv <- names(tab)[1:12]
        x[!(x %in% keep_lv)] <- "Other"
      }
      df[[cc]] <<- factor(x)
    }
    dropped
  }

  dropped_insts <- character(0)
  for (cc in inst_cols) dropped_insts <- prep_var(cc, dropped_insts)
  dropped_covs <- character(0)
  for (cc in cov_cols) dropped_covs <- prep_var(cc, dropped_covs)

  model_insts <- setdiff(inst_cols, dropped_insts)
  model_covs  <- setdiff(cov_cols, dropped_covs)
  if (length(model_insts) == 0) {
    stop(sprintf(
      "None of the mapped instrument columns(%s) is usable — each was constant, empty, or an identifier. An instrument must vary across rows.",
      paste(inst_names[inst_cols], collapse = ", ")))
  }

Step 4: Complete cases across everything the models need

keep_cols <- c("outcome", "treatment", model_insts, model_covs)
  df_clean <- df[, keep_cols, drop = FALSE]
  df_clean <- df_clean[complete.cases(df_clean), , drop = FALSE]
  final_rows <- nrow(df_clean)
  rows_removed <- initial_rows - final_rows

Step 5: Size guards, named in the user's own columns

if (final_rows < MIN_ROWS) {
    stop(sprintf(
      "Only %d rows have usable values in %s, %s, and the instrument column(s) %s. At least %d rows are required for a two-stage least squares estimate.",
      final_rows, outcome_name, treatment_name,
      paste(inst_names[model_insts], collapse = ", "), MIN_ROWS))
  }
  if (is.na(var(df_clean$treatment)) || isTRUE(var(df_clean$treatment) == 0)) {
    stop(sprintf("The treatment column(&#x27;%s') is constant — there is no variation whose effect could be estimated.",
                 treatment_name))
  }
  if (is.na(var(df_clean$outcome)) || isTRUE(var(df_clean$outcome) == 0)) {
    stop(sprintf("The outcome column(&#x27;%s') is constant — there is nothing to explain.",
                 outcome_name))
  }

  cov_rhs <- if (length(model_covs) > 0) paste(model_covs, collapse = " + ") else NULL

Step 6: FIRST STAGE — treatment on instruments (plus controls),

and the restricted fit that EXCLUDES the instruments. The difference in residual sum of squares gives the F statistic on the excluded instruments, which is the number the weak-instrument rule is about.

f_full <- as.formula(paste("treatment ~",
                             paste(c(model_insts, model_covs), collapse = " + ")))
  f_rest <- as.formula(paste("treatment ~", cov_rhs %||% "1"))
  fs      <- lm(f_full, data = df_clean)
  fs_rest <- lm(f_rest,  data = df_clean)

  n_inst_params <- fs$rank - fs_rest$rank
  if (n_inst_params < 1) {
    stop(sprintf(
      "The mapped instrument column(s) (%s) add no independent information beyond the control column(s) — they are perfectly explained by the controls, so no instrumented estimate is possible.",
      paste(inst_names[model_insts], collapse = ", ")))
  }
  rss_u <- sum(residuals(fs)^2)
  rss_r <- sum(residuals(fs_rest)^2)
  f_df1 <- n_inst_params
  f_df2 <- final_rows - fs$rank
  if (f_df2 < 10) {
    stop(sprintf(
      "Only %d rows remain against %d model terms — too few observations to estimate the effect of %s on %s with these instruments and controls.",
      final_rows, fs$rank, treatment_name, outcome_name))
  }
  f_stat <- ((rss_r - rss_u) / f_df1) / (rss_u / f_df2)
  f_p <- if (is.finite(f_stat)) pf(f_stat, f_df1, f_df2, lower.tail = FALSE) else NA_real_
  partial_r2 <- if (rss_r > 0) max(0, (rss_r - rss_u) / rss_r) else NA_real_

  strength_verdict <- if (!is.finite(f_stat)) "weak"
    else if (f_stat >= 10) "strong"
    else if (f_stat >= 5) "borderline"
    else "weak"
  usable <- is.finite(f_stat) && f_stat >= WEAK_F

Step 7: SECOND STAGE — outcome on the FITTED treatment (plus controls)

df_clean$treatment_fitted <- as.numeric(fitted(fs))
  f_ss <- as.formula(paste("outcome ~ treatment_fitted",
                           if (!is.null(cov_rhs)) paste("+", cov_rhs) else ""))
  ss <- lm(f_ss, data = df_clean)

  b <- coef(ss)
  keep_b <- !is.na(b)
  b <- b[keep_b]
  if (!("treatment_fitted" %in% names(b))) {
    stop(sprintf(
      "The predicted %s could not be separated from the control column(s) in the second stage — the instruments carry no variation that the controls do not already contain.",
      treatment_name))
  }
  Xhat <- model.matrix(ss)[, names(b), drop = FALSE]

Step 8: THE CORRECT TWO-STAGE LEAST SQUARES STANDARD ERRORS.

summary(ss) is WRONG here. Its residual variance comes from outcome - b'[1, fitted treatment, controls] but the standard error of a two-stage least squares estimate must use the STRUCTURAL residual outcome - b'[1, ACTUAL treatment, controls]. Same coefficients, different residuals. Everything else is identical: Var(b) = sigma^2_structural * (Xhat' Xhat)^-1 with sigma^2_structural = sum(u^2) / (n - k).

Xact <- Xhat
  Xact[, "treatment_fitted"] <- df_clean$treatment
  b_vec <- as.numeric(b)
  u_struct <- as.numeric(df_clean$outcome - Xact %*% b_vec)
  k_ss <- length(b)
  sigma2_2sls <- sum(u_struct^2) / (final_rows - k_ss)

  XtX <- crossprod(Xhat)
  XtX_inv <- tryCatch(chol2inv(chol(XtX)),
                      error = function(e) tryCatch(solve(XtX), error = function(e2) NULL))
  if (is.null(XtX_inv)) {
    stop("The second-stage design is rank deficient — the control column(s) and the predicted treatment carry the same information, so no standard error can be formed.")
  }
  V_2sls <- sigma2_2sls * XtX_inv
  se_all <- sqrt(pmax(0, diag(V_2sls)))
  names(se_all) <- names(b)

  iv_pos   <- which(names(b) == "treatment_fitted")
  iv_beta  <- as.numeric(b[iv_pos])
  iv_se    <- as.numeric(se_all[iv_pos])
  iv_t     <- if (is.finite(iv_se) && iv_se > 0) iv_beta / iv_se else NA_real_
  iv_p     <- if (is.na(iv_t)) NA_real_ else 2 * pt(-abs(iv_t), df = final_rows - k_ss)
  iv_crit  <- qt(0.975, df = final_rows - k_ss)
  iv_ci    <- iv_beta + c(-1, 1) * iv_crit * iv_se

Step 9: The ordinary least squares estimate — the biased comparison

f_ols <- as.formula(paste("outcome ~ treatment",
                            if (!is.null(cov_rhs)) paste("+", cov_rhs) else ""))
  ols <- lm(f_ols, data = df_clean)
  ols_coef <- summary(ols)$coefficients
  ols_beta <- as.numeric(ols_coef["treatment", "Estimate"])
  ols_se   <- as.numeric(ols_coef["treatment", "Std. Error"])
  ols_p    <- as.numeric(ols_coef["treatment", "Pr(>|t|)"])
  ols_ci   <- as.numeric(confint(ols)["treatment", ])
  gap <- iv_beta - ols_beta

Step 10: Over-identification (Sargan). Only defined when there are

MORE excluded instruments than endogenous regressors. With exactly one instrument the exclusion restriction is untestable — full stop.

overid_testable <- n_inst_params > 1
  sargan_stat <- NA_real_; sargan_df <- NA_integer_; sargan_p <- NA_real_
  if (overid_testable) {
    df_clean$iv_resid <- u_struct
    f_sar <- as.formula(paste("iv_resid ~",
                              paste(c(model_insts, model_covs), collapse = " + ")))
    sar <- tryCatch(lm(f_sar, data = df_clean), error = function(e) NULL)
    if (!is.null(sar)) {
      r2_sar <- summary(sar)$r.squared
      if (is.finite(r2_sar)) {
        sargan_stat <- final_rows * r2_sar
        sargan_df <- as.integer(n_inst_params - 1L)
        sargan_p <- pchisq(sargan_stat, df = sargan_df, lower.tail = FALSE)
      }
    }
  }

Step 11: Endogeneity (Wu-Hausman). Add the first-stage residual to

the ordinary regression: if its coefficient is distinguishable from zero, the treatment is endogenous and the ordinary estimate is biased.

df_clean$fs_resid <- as.numeric(residuals(fs))
  f_haus <- as.formula(paste("outcome ~ treatment + fs_resid",
                             if (!is.null(cov_rhs)) paste("+", cov_rhs) else ""))
  haus <- tryCatch(lm(f_haus, data = df_clean), error = function(e) NULL)
  hausman_t <- NA_real_; hausman_p <- NA_real_
  if (!is.null(haus)) {
    hc <- summary(haus)$coefficients
    if ("fs_resid" %in% rownames(hc)) {
      hausman_t <- as.numeric(hc["fs_resid", "t value"])
      hausman_p <- as.numeric(hc["fs_resid", "Pr(>|t|)"])
    }
  }

Step 12: First-stage coefficient table

fs_coef <- summary(fs)$coefficients
  sem_vars <- c(model_insts, model_covs)
  human_map <- as.list(c(inst_names[model_insts], cov_names[model_covs]))
  fs_terms <- rownames(fs_coef)
  fs_terms <- fs_terms[fs_terms != "(Intercept)"]
  first_stage_df <- data.frame(
    term = sapply(fs_terms, function(tm) humanize_term(tm, sem_vars, human_map),
                  USE.NAMES = FALSE),
    role = sapply(fs_terms, function(tm) {
      hits <- sem_vars[startsWith(tm, sem_vars)]
      v <- if (length(hits)) hits[which.max(nchar(hits))] else ""
      if (v %in% model_insts) "Instrument" else "Control"
    }, USE.NAMES = FALSE),
    estimate  = signif(as.numeric(fs_coef[fs_terms, "Estimate"]), 4),
    std_error = signif(as.numeric(fs_coef[fs_terms, "Std. Error"]), 4),
    t_stat    = signif(as.numeric(fs_coef[fs_terms, "t value"]), 4),
    p_value   = sapply(as.numeric(fs_coef[fs_terms, "Pr(>|t|)"]), fmt_p),
    stringsAsFactors = FALSE
  )
  rownames(first_stage_df) <- NULL

Step 13: First-stage picture. Prefer the numeric instrument with the

largest absolute t statistic; if no instrument is numeric, fall back to the fitted treatment. LAT-1445 guard: filter NA before taking a maximum.

numeric_insts <- model_insts[sapply(model_insts, function(cc) is.numeric(df_clean[[cc]]))]
  fs_plot_kind <- "fitted"
  plot_inst <- NULL
  if (length(numeric_insts) > 0) {
    tvals <- sapply(numeric_insts, function(cc) {
      if (cc %in% rownames(fs_coef)) abs(as.numeric(fs_coef[cc, "t value"])) else NA_real_
    })
    ok <- which(is.finite(tvals))
    if (length(ok) > 0) {
      plot_inst <- numeric_insts[ok[which.max(tvals[ok])]]
      fs_plot_kind <- "instrument"
    }
  }
  x_vals <- if (fs_plot_kind == "instrument") df_clean[[plot_inst]] else df_clean$treatment_fitted
  set.seed(42)
  sidx <- if (final_rows > 1000) sort(sample(final_rows, 1000)) else seq_len(final_rows)
  first_stage_points <- data.frame(
    instrument_value = round(as.numeric(x_vals[sidx]), 4),
    treatment_value  = round(as.numeric(df_clean$treatment[sidx]), 4),
    stringsAsFactors = FALSE
  )
  first_stage_points <- first_stage_points[order(first_stage_points$instrument_value), ]
  rownames(first_stage_points) <- NULL
  fs_plot_label <- if (fs_plot_kind == "instrument") inst_names[[plot_inst]]
                   else paste0("Predicted ", treatment_name)

Step 14: Estimate tables

ols_sig <- is.finite(ols_ci[1]) && is.finite(ols_ci[2]) &&
    (ols_ci[1] > 0 || ols_ci[2] < 0)
  iv_sig <- is.finite(iv_ci[1]) && is.finite(iv_ci[2]) &&
    (iv_ci[1] > 0 || iv_ci[2] < 0)

  estimate_compare_df <- data.frame(
    estimate_type = c("Ordinary least squares", "Instrumental variables(2SLS)"),
    estimate = signif(c(ols_beta, iv_beta), 4),
    ci_low   = signif(c(ols_ci[1], iv_ci[1]), 4),
    ci_high  = signif(c(ols_ci[2], iv_ci[2]), 4),
    stringsAsFactors = FALSE
  )

  estimate_results_df <- data.frame(
    estimate_type = c("Ordinary least squares", "Instrumental variables(2SLS)"),
    estimate  = signif(c(ols_beta, iv_beta), 4),
    std_error = signif(c(ols_se, iv_se), 4),
    ci_low    = signif(c(ols_ci[1], iv_ci[1]), 4),
    ci_high   = signif(c(ols_ci[2], iv_ci[2]), 4),
    p_value   = c(fmt_p(ols_p), fmt_p(iv_p)),
    interpretation = c(
      paste0("Confounded: mixes the effect of ", treatment_name,
             " with anything unmeasured that moves both it and ", outcome_name, "."),
      if (usable)
        paste0("Uses only the variation in ", treatment_name,
               " driven by the instrument; valid only if the exclusion restriction holds.")
      else
        paste0("NOT USABLE — the first-stage F of ", fmt_val(f_stat),
               " is below 10, so this estimate is biased toward the ordinary one and its interval is too narrow.")
    ),
    stringsAsFactors = FALSE
  )

Step 15: Diagnostic table — the checks that decide whether the

instrumented estimate may be believed at all.

inst_list <- paste(inst_names[model_insts], collapse = ", ")
  diag_rows <- list(
    data.frame(
      check = "Instrument strength(first-stage F on the excluded instruments)",
      result = paste0("F = ", fmt_val(f_stat), " on ", f_df1, " and ",
                      format(f_df2, big.mark = ","), " degrees of freedom, ",
                      p_phrase(f_p)),
      interpretation = if (usable)
        paste0("Above the conventional threshold of 10, so ", inst_list,
               " moves ", treatment_name, " strongly enough for the instrumented estimate to be read.")
      else
        paste0("Below the conventional threshold of 10. With an instrument this weak, two-stage least squares is biased back toward the ordinary least squares estimate it was meant to correct, and its confidence interval covers the truth less often than 95 percent of the time. The instrumented estimate is not usable."),
      stringsAsFactors = FALSE
    ),
    data.frame(
      check = "Instrument relevance(extra variation explained)",
      result = paste0(fmt_val(100 * partial_r2), " percent of the variation in ",
                      treatment_name, " left by the controls"),
      interpretation = paste0("The share of ", treatment_name,
                              " that the instrument explains beyond the control columns. This is the only part of ",
                              treatment_name, " the instrumented estimate uses."),
      stringsAsFactors = FALSE
    ),
    data.frame(
      check = "Over-identification(Sargan test of the exclusion restriction)",
      result = if (overid_testable)
        paste0("Statistic = ", fmt_val(sargan_stat), " on ", sargan_df,
               " degrees of freedom, ", p_phrase(sargan_p))
      else
        "Not testable with one instrument",
      interpretation = if (!overid_testable)
        paste0("With exactly one instrument and one confounded regressor the model is exactly identified, and the exclusion restriction is UNTESTABLE — no statistic in this or any other analysis can check it. It is assumed, not established.")
      else if (is.finite(sargan_p) && sargan_p < 0.05)
        paste0("The test REJECTS: the instruments(", inst_list,
               ") disagree with each other about the effect of ", treatment_name,
               " by more than sampling noise, so at least one of them fails the exclusion restriction. The instrumented estimate should not be trusted until the offending instrument is identified and removed.")
      else
        paste0("The test does not reject, which is consistent with the instruments agreeing. This is a weak check, not a clearance: it can only detect DISAGREEMENT between instruments, and instruments that are all invalid in the same direction pass it."),
      stringsAsFactors = FALSE
    ),
    data.frame(
      check = "Endogeneity(Wu-Hausman test of ordinary vs instrumented)",
      result = if (is.finite(hausman_p))
        paste0("t = ", fmt_val(hausman_t), ", ", p_phrase(hausman_p))
      else "Not available",
      interpretation = if (is.finite(hausman_p) && hausman_p < 0.05)
        paste0("The two estimates differ by more than sampling noise, which is evidence that ",
               treatment_name, " is confounded and the ordinary least squares estimate is biased.")
      else if (is.finite(hausman_p))
        paste0("The two estimates are not distinguishable from each other, so this data gives no evidence that ",
               treatment_name, " is confounded. When that is so, the ordinary least squares estimate is the more precise of the two.")
      else "The endogeneity test could not be computed on this data.",
      stringsAsFactors = FALSE
    ),
    data.frame(
      check = "Exclusion restriction",
      result = "Assumed, never tested",
      interpretation = paste0(
        "The analysis assumes ", inst_list, " affects ", outcome_name,
        " ONLY by changing ", treatment_name, ", and is unrelated to whatever else moves ",
        outcome_name, ". No amount of data can verify that. It is a claim about the world that you make, and the whole estimate rests on it."),
      stringsAsFactors = FALSE
    ),
    data.frame(
      check = "What the estimate applies to",
      result = "A local average treatment effect",
      interpretation = paste0(
        "Two-stage least squares recovers the effect of ", treatment_name, " on ",
        outcome_name, " for the subgroup whose ", treatment_name,
        " actually responded to ", inst_list,
        " — not the average effect across everyone in the data. If the effect differs between people who respond to the instrument and people who do not, this number does not describe the second group. It also assumes the instrument pushes every unit in the same direction."),
      stringsAsFactors = FALSE
    )
  )
  diagnostic_df <- do.call(rbind, diag_rows)
  rownames(diagnostic_df) <- NULL

Step 16: KPI metrics (user-facing keys)

metrics <- list(
    `Observations`          = final_rows,
    `OLS Estimate`          = signif(ols_beta, 4),
    `IV Estimate`           = signif(iv_beta, 4),
    `Confounding Gap`       = signif(gap, 4),
    `First-Stage F`         = signif(f_stat, 4),
    `Instrument Strength`   = if (usable) "strong enough" else "too weak"
  )

Step 17: json_output machine channel. When the instrument is weak the

refusal comes FIRST — a number a reader should not use must not be the first thing they read.

gap_direction <- if (!is.finite(gap)) "differs from"
    else if (abs(gap) < 1e-12) "matches"
    else if (gap > 0) "is higher than" else "is lower than"

  answer_head <- if (!usable) {
    paste0(
      "NOT USABLE: the instrument", if (length(model_insts) > 1) "s" else "", " ",
      inst_list, " ", if (length(model_insts) > 1) "are" else "is",
      " too weak — the first-stage F on the excluded instrument", if (f_df1 > 1) "s" else "",
      " is ", fmt_val(f_stat), ", below the conventional threshold of 10. ",
      "At that strength two-stage least squares is biased back toward the ordinary least squares estimate it was meant to correct and its confidence interval is too narrow, so the instrumented number below should be read as a diagnostic, not as an answer. ")
  } else ""

  json_output <- list(
    answer = paste0(
      answer_head,
      "Two-stage least squares of ", outcome_name, " on ", treatment_name,
      ", instrumented by ", inst_list,
      if (length(model_covs) > 0)
        paste0(" and controlling for ", paste(cov_names[model_covs], collapse = ", "))
      else "",
      ", across ", format(final_rows, big.mark = ","), " rows: the ordinary least squares estimate is ",
      fmt_signed(ols_beta), " (95% CI ", fmt_val(ols_ci[1]), " to ", fmt_val(ols_ci[2]),
      ") and the instrumented estimate is ", fmt_signed(iv_beta),
      " (95% CI ", fmt_val(iv_ci[1]), " to ", fmt_val(iv_ci[2]), ", ", p_phrase(iv_p),
      "), a difference of ", fmt_signed(gap),
      " — the instrumented estimate ", gap_direction, " the ordinary one, and that distance is what the confounding was doing to it. ",
      "The first-stage F on the excluded instrument", if (f_df1 > 1) "s" else "", " is ",
      fmt_val(f_stat), " (", if (usable) "at or above" else "below",
      " the conventional threshold of 10). ",
      if (overid_testable) {
        if (is.finite(sargan_p) && sargan_p < 0.05)
          paste0("The Sargan over-identification test rejects(", p_phrase(sargan_p),
                 "), so the instruments disagree with each other and at least one of them fails the exclusion restriction. ")
        else
          paste0("The Sargan over-identification test does not reject(", p_phrase(sargan_p),
                 "), which is consistent with the instruments agreeing but does not clear them. ")
      } else {
        "With exactly one instrument the exclusion restriction is untestable — no statistic here can check it. "
      },
      "The exclusion restriction — that ", inst_list, " affects ", outcome_name,
      " only through ", treatment_name,
      " — is assumed, not established, and the estimate applies to the subgroup whose ",
      treatment_name, " responded to the instrument rather than to everyone in the data."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "instrument_strength",
        "first_stage_plot", "estimate_comparison", "estimate_table", "diagnostics"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows, rows_removed = rows_removed,
    outcome_name = outcome_name, treatment_name = treatment_name,
    inst_names = inst_names, cov_names = cov_names,
    model_insts = model_insts, model_covs = model_covs,
    dropped_insts = dropped_insts, dropped_covs = dropped_covs,
    inst_list = inst_list,
    ols_beta = ols_beta, ols_se = ols_se, ols_ci = ols_ci, ols_p = ols_p,
    ols_sig = ols_sig,
    iv_beta = iv_beta, iv_se = iv_se, iv_ci = iv_ci, iv_p = iv_p, iv_sig = iv_sig,
    iv_se_naive = iv_se_naive, sigma2_2sls = sigma2_2sls,
    gap = gap, gap_direction = gap_direction,
    f_stat = f_stat, f_p = f_p, f_df1 = f_df1, f_df2 = f_df2,
    partial_r2 = partial_r2, strength_verdict = strength_verdict, usable = usable,
    n_inst_params = n_inst_params, overid_testable = overid_testable,
    sargan_stat = sargan_stat, sargan_df = sargan_df, sargan_p = sargan_p,
    hausman_t = hausman_t, hausman_p = hausman_p,
    first_stage_df = first_stage_df, first_stage_points = first_stage_points,
    fs_plot_kind = fs_plot_kind, fs_plot_label = fs_plot_label,
    estimate_compare_df = estimate_compare_df,
    estimate_results_df = estimate_results_df,
    diagnostic_df = diagnostic_df,
    metrics = metrics, json_output = json_output
  )
}

This is the exact R code executed to produce the analysis above. Run it with the same dataset to reproduce these results.

Your data has more stories to tell. Run any analysis on your own data — validated R modules, interactive reports, AI insights, and PDF export. 500 free credits on signup.
Try Free — No Signup Sign Up Free
Something wrong? Report it and we'll fix it — or get your credits back
Download PDF · Print · Share ·Refresh this analysis ·Edit this report

Report an Issue

Tell us what's wrong. You'll get a free re-run of this analysis so you can try again with different parameters. If the re-run still doesn't meet your expectations, we'll refund your credits.

Want to run this analysis on your own data? Upload CSV — Free Analysis See Pricing