Executive Summary
Ordinary versus instrumented effect of years education on log wage
The short answer
College proximity as an instrument yields a causal return to education of +0.1323 log points per year, compared to +0.074 from ordinary regression. The 0.0583 difference reveals confounding: unmeasured traits that raise both education and wages were inflating the ordinary estimate downward.
The detail
Across 3,010 observations, ordinary least squares estimates the education effect at 0.074 per year (95% CI 0.0671 to 0.0809). The instrumented estimate is 0.1323 (95% CI 0.0358 to 0.229, p = 0.00725). The first-stage F of 16.7 clears the threshold of 10, indicating the instrument is strong enough. The confounding gap of 0.0583 is the magnitude of bias in the ordinary estimate. The exclusion restriction—that near 4yr college affects log wage only through years education—cannot be tested with one instrument; it is an untestable assumption. The instrumented effect applies to the local-average-treatment-effect population: those whose education actually responded to college proximity.
What this can't tell you
Whether college proximity has any effect on wages independent of education. A single instrument cannot test the exclusion restriction. Causation is only as valid as this assumption holds.
Analysis Overview
Two-stage least squares of log wage on years education, instrumented by near 4yr college, across 3,010 observations.
The short answer
Ordinary regression of wages on education is confounded: it mixes the true effect of education with unmeasured traits that raise both. This analysis isolates the causal effect using college proximity as an instrument, which explains only the part of education variation tied to proximity, leaving confounding traits behind.
The detail
Two-stage least squares across 3,010 observations with near 4yr college as the single instrument and experience, experience squared, race, region, and urbanization as controls. The first stage regresses years education on the instrument; the second regresses log wage on the predicted education from stage one. The first-stage F-statistic is 16.72, meeting the conventional threshold of 10 for instrument strength. Standard errors are rebuilt from residuals against actual years education, not the second stage's predicted values, because the latter would misstate precision in two-stage least squares.
What this can't tell you
The exclusion restriction—that college proximity affects wages only through education—is untestable with a single instrument. No statistic here can verify that proximity does not independently influence wages through networks, local labor-market conditions, or other unmeasured pathways. The causal claim rests entirely on this assumption.
Data Quality
Row accounting, column typing, and which mapped columns were excluded.
The short answer
All 3,010 rows were complete on log wage, years education, and the instrument. No rows were dropped. The analysis uses 1 instrument (near 4yr college) and 5 controls (experience, experience squared, black, south, urban smsa), all of which entered the models without issue.
The detail
The dataset loaded with 3,010 rows and remained at 3,010 rows after completeness checking: all observations had non-missing values for the outcome (log wage), treatment (years education), and instrument (near 4yr college). Numeric columns were median-filled for missing values before the completeness check; text columns were treated as categories with blanks coded as 'Missing' and rare levels beyond twelve grouped into 'Other'. All mapped columns were usable; none were excluded as constant or identifier-like.
What this can't tell you
The preprocessing report does not address whether the median imputation for numeric columns biased the estimates, or whether any rows should have been excluded on substantive grounds (e.g., implausible wage or education values). A finer-grained account of which columns received median imputation and how many values were filled would clarify the scope of that step.
First Stage
First-stage coefficients and the F statistic on the excluded instruments.
| Term | Role | Estimate | Std Error | T Stat | P Value |
|---|---|---|---|---|---|
| near 4yr college | Instrument | 0.3373 | 0.0825 | 4.089 | < 0.001 |
| experience | Control | -0.41 | 0.0337 | -12.17 | < 0.001 |
| experience sq | Control | 0.0007 | 0.0016 | 0.4438 | 0.657 |
| black | Control | -1.006 | 0.0896 | -11.22 | < 0.001 |
| south | Control | -0.2915 | 0.0792 | -3.679 | < 0.001 |
| urban smsa | Control | 0.4039 | 0.0849 | 4.758 | < 0.001 |
The short answer
College proximity moves years of education with a first-stage F of 16.7, well above the threshold of 10. The relationship is statistically significant (t = 4.089, p < 0.001) and explains 0.554 percent of the remaining variation in education after controls. The instrument is strong enough for the analysis to proceed.
The detail
The first stage regresses years education on near 4yr college, experience, experience sq, black, south, and urban smsa. The F statistic on the excluded instrument is 16.7 on 1 and 3,003 degrees of freedom (p < 0.001), clearing the conventional threshold of 10. The near 4yr college coefficient is 0.3373 (SE 0.0825, t = 4.089). This instrument explains 0.554 percent of the variation in years education beyond the control columns. A large first-stage F indicates relevance—the instrument moves the treatment—but says nothing about validity, which rests on the untestable exclusion restriction.
What this can't tell you
Whether college proximity is valid—that is, whether it reaches log wage only through education. A strong first stage and an invalid instrument are indistinguishable in the data. The analysis cannot test the exclusion restriction with one instrument.
Instrument versus Treatment
The first-stage relationship the instrumented estimate runs through.
The short answer
The scatter shows a visible upward slope: people near a 4-year college report more years of education than those far away. The cloud is diffuse but directional, consistent with the first-stage F of 16.7. The pattern is the mechanism through which the instrumented estimate runs.
The detail
The plot samples 1,000 of the 3,010 rows, with near 4yr college on the horizontal axis (binary: 0 or 1) and years education on the vertical. The slope is the first-stage relationship; the F-statistic of 16.7 quantifies its strength. A tighter band would indicate the instrument explains more of education variation; this diffuse cloud means college proximity explains only a portion, so the instrumented estimate is correspondingly less precise.
What this can't tell you
Whether the instrument is valid. A strong first stage and a weak but invalid instrument are visually identical. The exclusion restriction cannot be assessed from this plot alone.
Ordinary versus Instrumented
The confounded estimate next to the instrumented one, with 95 percent intervals.
The short answer
The instrumented effect (+0.1323) is 0.0583 larger in magnitude than the ordinary estimate (+0.074), indicating the ordinary regression was downward-biased by confounding. The instrumented interval is wider because the estimate uses only the education variation college proximity explains.
The detail
Ordinary least squares: 0.074 per year (95% CI 0.0671 to 0.0809). Instrumental variables: 0.1323 per year (95% CI 0.0357 to 0.2288). The 0.0583 gap is the confounding the instrument removed. The instrumented interval is wider because the estimate is built from only the instrumented share of years education variation. Both intervals are bounded away from zero on the lower end, though the instrumented interval's upper bound is 0.2288.
What this can't tell you
Whether the confounding direction or magnitude would hold under alternative instruments or in other populations. The validity rests on the untestable exclusion restriction.
Effect Estimates
Ordinary and instrumented estimates with intervals and p-values.
| Estimate Type | Estimate | Std Error | CI Low | CI High | P Value | Interpretation |
|---|---|---|---|---|---|---|
| Ordinary least squares | 0.074 | 0.0035 | 0.0671 | 0.0809 | < 0.001 | Confounded: mixes the effect of years education with anything unmeasured that moves both it and log wage. |
| Instrumental variables (2SLS) | 0.1323 | 0.0492 | 0.0357 | 0.2288 | 0.00725 | Uses only the variation in years education driven by the instrument; valid only if the exclusion restriction holds. |
The short answer
The instrumented estimate is 0.1323 per year of education (SE 0.0492, p = 0.00725), compared to 0.074 in ordinary regression (SE 0.0035, p < 0.001). The instrumented standard error is rebuilt from residuals against actual education, not the second stage's own regression residuals, which would have yielded 0.0504—a difference of 0.976 times.
The detail
Ordinary least squares: estimate 0.074, SE 0.0035, 95% CI 0.0671 to 0.0809, p < 0.001. This uses all 3,010 rows and all variation in years education, including the confounded part. Instrumental variables (2SLS): estimate 0.1323, SE 0.0492, 95% CI 0.0357 to 0.2288, p = 0.00725. This uses only the 0.554 percent of education variation attributable to near 4yr college. The reported standard error is the two-stage least squares standard error, measured against actual years education; the second stage's own figure of 0.0504 would understate the true variance.
What this can't tell you
Whether confounding is present—the endogeneity test did not detect it (p = 0.215). Whether the exclusion restriction holds. Whether the effect applies to all workers or only those whose education was influenced by college proximity.
Assumptions & Diagnostics
What was tested, what was assumed, and what the estimate actually applies to.
| Check | Result | Interpretation |
|---|---|---|
| Instrument strength (first-stage F on the excluded instruments) | F = 16.7 on 1 and 3,003 degrees of freedom, p < 0.001 | Above the conventional threshold of 10, so near 4yr college moves years education strongly enough for the instrumented estimate to be read. |
| Instrument relevance (extra variation explained) | 0.554 percent of the variation in years education left by the controls | The share of years education that the instrument explains beyond the control columns. This is the only part of years education the instrumented estimate uses. |
| Over-identification (Sargan test of the exclusion restriction) | Not testable with one instrument | 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. |
| Endogeneity (Wu-Hausman test of ordinary vs instrumented) | t = -1.24, p = 0.215 | The two estimates are not distinguishable from each other, so this data gives no evidence that years education is confounded. When that is so, the ordinary least squares estimate is the more precise of the two. |
| Exclusion restriction | Assumed, never tested | The analysis assumes near 4yr college affects log wage ONLY by changing years education, and is unrelated to whatever else moves log wage. 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 to | A local average treatment effect | Two-stage least squares recovers the effect of years education on log wage for the subgroup whose years education actually responded to near 4yr 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. |
The short answer
Instrument strength passes: the first-stage F of 16.7 clears the threshold of 10. The exclusion restriction—that college proximity affects wages only through education—is untestable with one instrument and cannot be verified. The Wu-Hausman endogeneity test (p = 0.215) found no evidence of confounding, so ordinary regression would be more precise if confounding is absent.
The detail
Relevance (testable): first-stage F = 16.7 on 1 and 3,003 degrees of freedom, p < 0.001. The instrument is strong enough. Over-identification (testable): not applicable; the model is exactly identified with one instrument and one confounded regressor, so the exclusion restriction is untestable. Endogeneity (testable): Wu-Hausman t = −1.24, p = 0.215. The two estimates are not distinguishable, so no evidence of confounding emerges from the data. Exclusion restriction (assumed, untestable): the analysis assumes near 4yr college affects log wage only by changing years education. Monotonicity (assumed, untestable): the analysis assumes the instrument pushes every unit in the same direction. Local average treatment effect: the estimate applies to the subgroup whose education responded to proximity, not the average effect across all workers.
What this can't tell you
Whether college proximity has unmeasured pathways to wages. Whether the effect is the same for compliers (those whose education responds to the instrument) and non-compliers. The endogeneity test's non-significance suggests confounding may not be present; if so, the ordinary estimate is more reliable.
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 'outcome' column (the numeric result you want to explain)")
}
if (!"treatment" %in% names(df)) {
stop("column_mapping must map a '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('%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_rowsStep 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('%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('%s') is constant — there is nothing to explain.",
outcome_name))
}
cov_rhs <- if (length(model_covs) > 0) paste(model_covs, collapse = " + ") else NULLStep 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_FStep 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_seStep 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_betaStep 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) <- NULLStep 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) <- NULLStep 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
)
}