Econometric Models

The Linear Model, Binary Choice, and a Tour of the Toolbox

What Is an Econometric Model?

// ---------------------------------------------------------------------------
// ldk — a tiny, dependency-free plotting + statistics library for these decks.
//
// Everything runs in the browser as plain JavaScript building SVG nodes, so a
// slider move costs a fraction of a millisecond instead of a round trip
// through a WebAssembly R runtime.
// ---------------------------------------------------------------------------
ldk = {
  const NS = "http://www.w3.org/2000/svg";
  let uid = 0;

  // Deck palette, matching ldkslides.css
  const C = {
    green:  "#74b83d",
    blue:   "#4072c2",
    pink:   "#ED017D",
    orange: "#d18b2a",
    grey:   "#6E6E6E",
    ink:    "#222222",
    axis:   "#9aa0a6",
    faint:  "#efefef"
  };

  function el(tag, attrs, text) {
    const e = document.createElementNS(NS, tag);
    if (attrs) for (const k in attrs) if (attrs[k] != null) e.setAttribute(k, attrs[k]);
    if (text != null)
      // SVG collapses runs of whitespace, which would eat the wide gaps that
      // separate the readout columns; non-breaking spaces survive. (xml:space
      // is a namespaced attribute and does not survive setAttribute here.)
      e.textContent = tag === "text"
        ? String(text).replace(/ {2,}/g, m => " ".repeat(m.length))
        : text;
    return e;
  }

  function alpha(hex, a) {
    const n = parseInt(hex.slice(1), 16);
    return `rgba(${(n >> 16) & 255},${(n >> 8) & 255},${n & 255},${a})`;
  }

  // ------------------------------------------------------------------ numerics

  const LGC = [676.5203681218851, -1259.1392167224028, 771.32342877765313,
               -176.61502916214059, 12.507343278686905, -0.13857109526572012,
               9.9843695780195716e-6, 1.5056327351493116e-7];

  function lgamma(z) {
    if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);
    z -= 1;
    let x = 0.99999999999980993;
    for (let i = 0; i < 8; i++) x += LGC[i] / (z + i + 1);
    const t = z + 7.5;
    return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
  }

  const lchoose = (n, k) => lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1);

  // Chebyshev erfc (Numerical Recipes); plenty of accuracy for plotting.
  const ERFC_COF = [
    -1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2,
    -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4,
     4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6,
     1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
     6.529054439e-9, 5.059343495e-9, -9.91364156e-10, -2.27365122e-10,
     9.6467911e-11, 2.394038e-12, -6.886027e-12, 8.94487e-13,
     3.13092e-13, -1.12708e-13, 3.81e-16, 7.106e-15];

  function erfc(x) {
    const z = Math.abs(x), t = 2 / (2 + z), ty = 4 * t - 2;
    let d = 0, dd = 0;
    for (let j = ERFC_COF.length - 1; j > 0; j--) {
      const tmp = d;
      d = ty * d - dd + ERFC_COF[j];
      dd = tmp;
    }
    const ans = t * Math.exp(-z * z + 0.5 * (ERFC_COF[0] + ty * d) - dd);
    return x >= 0 ? ans : 2 - ans;
  }

  const dnorm = (x, m = 0, s = 1) => {
    const z = (x - m) / s;
    return Math.exp(-0.5 * z * z) / (s * Math.sqrt(2 * Math.PI));
  };
  const pnorm = (x, m = 0, s = 1) => 0.5 * erfc(-(x - m) / (s * Math.SQRT2));

  // Acklam's inverse normal CDF (relative error < 1.2e-9)
  const QN_A = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
                 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00];
  const QN_B = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
                 6.680131188771972e+01, -1.328068155288572e+01];
  const QN_C = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
                -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00];
  const QN_D = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
                3.754408661907416e+00];

  function qnorm(p, m = 0, s = 1) {
    if (p <= 0) return -Infinity;
    if (p >= 1) return Infinity;
    const lo = 0.02425, hi = 1 - lo;
    let x, q, r;
    if (p < lo) {
      q = Math.sqrt(-2 * Math.log(p));
      x = (((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
          ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    } else if (p <= hi) {
      q = p - 0.5; r = q * q;
      x = (((((QN_A[0] * r + QN_A[1]) * r + QN_A[2]) * r + QN_A[3]) * r + QN_A[4]) * r + QN_A[5]) * q /
          (((((QN_B[0] * r + QN_B[1]) * r + QN_B[2]) * r + QN_B[3]) * r + QN_B[4]) * r + 1);
    } else {
      q = Math.sqrt(-2 * Math.log(1 - p));
      x = -(((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
           ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    }
    return m + s * x;
  }

  function dbinom(k, n, p) {
    if (k < 0 || k > n) return 0;
    if (p <= 0) return k === 0 ? 1 : 0;
    if (p >= 1) return k === n ? 1 : 0;
    return Math.exp(lchoose(n, k) + k * Math.log(p) + (n - k) * Math.log1p(-p));
  }
  const dpois = (k, lam) => k < 0 ? 0 : Math.exp(-lam + k * Math.log(lam) - lgamma(k + 1));
  const dexp  = (x, rate) => x < 0 ? 0 : rate * Math.exp(-rate * x);
  const pexp  = (x, rate) => x <= 0 ? 0 : 1 - Math.exp(-rate * x);
  const dunif = (x, a, b) => (x >= a && x <= b) ? 1 / (b - a) : 0;
  const punif = (x, a, b) => x <= a ? 0 : x >= b ? 1 : (x - a) / (b - a);
  const dlnorm = (x, ml = 0, sl = 1) => x <= 0 ? 0 : dnorm(Math.log(x), ml, sl) / x;
  const dchisq1 = (x) => x <= 0 ? 0 : dnorm(Math.sqrt(x)) / Math.sqrt(x);

  // ---- regularized lower incomplete gamma P(a,x): series + continued fraction
  function gammaP(a, x) {
    if (!(x >= 0) || !(a > 0)) return NaN;
    if (x === 0) return 0;
    const pref = Math.exp(-x + a * Math.log(x) - lgamma(a));
    if (x < a + 1) {                       // series
      let ap = a, sum = 1 / a, del = sum;
      for (let i = 0; i < 500; i++) {
        ap++; del *= x / ap; sum += del;
        if (Math.abs(del) < Math.abs(sum) * 1e-16) break;
      }
      return sum * pref;
    }
    const TINY = 1e-300;                   // continued fraction for Q = 1 - P
    let b = x + 1 - a, c = 1 / TINY, d = 1 / b, h = d;
    for (let i = 1; i <= 500; i++) {
      const an = -i * (i - a);
      b += 2;
      d = an * d + b; if (Math.abs(d) < TINY) d = TINY;
      c = b + an / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 1e-16) break;
    }
    return 1 - pref * h;
  }

  // ---- regularized incomplete beta I_x(a,b)
  function betacf(a, b, x) {
    const TINY = 1e-300, qab = a + b, qap = a + 1, qam = a - 1;
    let c = 1, d = 1 - qab * x / qap;
    if (Math.abs(d) < TINY) d = TINY;
    d = 1 / d;
    let h = d;
    for (let m = 1; m <= 400; m++) {
      const m2 = 2 * m;
      let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d; h *= d * c;
      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 3e-16) break;
    }
    return h;
  }
  function betai(a, b, x) {
    if (!(x > 0)) return 0;
    if (x >= 1) return 1;
    const bt = Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
                        a * Math.log(x) + b * Math.log1p(-x));
    return x < (a + 1) / (a + b + 2)
      ? bt * betacf(a, b, x) / a
      : 1 - bt * betacf(b, a, 1 - x) / b;
  }

  // ---- chi-squared and Student t
  const dchisq = (x, k) => x <= 0 ? 0 :
    Math.exp((k / 2 - 1) * Math.log(x) - x / 2 - lgamma(k / 2) - (k / 2) * Math.LN2);
  const pchisq = (x, k) => x <= 0 ? 0 : gammaP(k / 2, x / 2);

  const dt = (x, nu) => Math.exp(lgamma((nu + 1) / 2) - lgamma(nu / 2)
    - 0.5 * Math.log(nu * Math.PI) - (nu + 1) / 2 * Math.log1p(x * x / nu));
  const pt = (x, nu) => {
    const half = 0.5 * betai(nu / 2, 0.5, nu / (nu + x * x));
    return x > 0 ? 1 - half : half;
  };

  // quantiles by bisection — a handful of calls per redraw, so speed is moot
  function invert(cdf, p, lo, hi) {
    for (let i = 0; i < 90; i++) {
      const mid = (lo + hi) / 2;
      if (cdf(mid) < p) lo = mid; else hi = mid;
    }
    return (lo + hi) / 2;
  }
  const qchisq = (p, k) => invert(x => pchisq(x, k), p, 0, Math.max(60, k * 25));
  const qt = (p, nu) => invert(x => pt(x, nu), p, -1e4, 1e4);

  // Cumulative sums of a pmf enumerated over 0..n — cheaper and more accurate
  // here than a general incomplete-beta / incomplete-gamma implementation.
  function cdfFromPmf(pmf) {
    const out = new Array(pmf.length);
    let s = 0;
    for (let i = 0; i < pmf.length; i++) { s += pmf[i]; out[i] = Math.min(s, 1); }
    return out;
  }

  // ----------------------------------------------------------------- sequences

  function seq(a, b, n) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = a + (b - a) * i / (n - 1);
    return out;
  }
  const cumsum = (v) => { let s = 0; return v.map(x => (s += x)); };
  const sum = (v) => v.reduce((a, b) => a + b, 0);
  const mean = (v) => sum(v) / v.length;

  function sd(v) {
    const m = mean(v);
    return Math.sqrt(v.reduce((a, x) => a + (x - m) * (x - m), 0) / (v.length - 1));
  }

  function quantileSorted(s, p) {
    const h = (s.length - 1) * p, i = Math.floor(h);
    const j = Math.min(i + 1, s.length - 1);
    return s[i] + (h - i) * (s[j] - s[i]);
  }

  function hist(v, lo, hi, nb) {
    const cnt = new Float64Array(nb), w = (hi - lo) / nb;
    for (let i = 0; i < v.length; i++) {
      const k = Math.floor((v[i] - lo) / w);
      if (k >= 0 && k < nb) cnt[k]++;
      else if (v[i] === hi) cnt[nb - 1]++;
    }
    const dens = new Array(nb);
    for (let i = 0; i < nb; i++) dens[i] = cnt[i] / (v.length * w);
    return { lo, hi, w, nb, dens, count: cnt };
  }

  // --------------------------------------------------------------------- rngs

  // mulberry32: small, fast, seedable — the demos stay reproducible.
  function rng(seed) {
    let a = seed >>> 0;
    const u = () => {
      a |= 0; a = (a + 0x6D2B79F5) | 0;
      let t = Math.imul(a ^ (a >>> 15), 1 | a);
      t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
      return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
    };
    let spare = null;
    u.normal = () => {                       // Marsaglia polar
      if (spare !== null) { const v = spare; spare = null; return v; }
      let x, y, s;
      do { x = 2 * u() - 1; y = 2 * u() - 1; s = x * x + y * y; } while (s >= 1 || s === 0);
      const f = Math.sqrt(-2 * Math.log(s) / s);
      spare = y * f;
      return x * f;
    };
    u.exp = () => -Math.log(1 - u());
    u.cauchy = () => Math.tan(Math.PI * (u() - 0.5));
    return u;
  }

  // Weighted sampling without replacement (Efraimidis–Spirakis exponential race)
  function sampleWeighted(n, w, u) {
    const N = w.length;
    const keys = new Array(N);
    for (let i = 0; i < N; i++) keys[i] = w[i] > 0 ? -Math.log(1 - u()) / w[i] : Infinity;
    const idx = Array.from({ length: N }, (_, i) => i);
    idx.sort((a, b) => keys[a] - keys[b]);
    return idx.slice(0, n);
  }

  // Sampling *with* replacement — what the bootstrap needs, and what
  // sampleWeighted deliberately is not.
  function sampleReplace(n, v, u) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = v[Math.floor(u() * v.length)];
    return out;
  }

  // Beta density — the conjugate posterior for a proportion.
  const dbeta = (x, a, b) => (x <= 0 || x >= 1) ? 0 :
    Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
             (a - 1) * Math.log(x) + (b - 1) * Math.log(1 - x));

  // Beta quantile, by inverting the regularized incomplete beta above. Used for
  // the Clopper-Pearson interval, whose endpoints are exactly beta quantiles.
  const qbeta = (p, a, b) => invert(x => betai(a, b, x), p, 0, 1);

  // Gini coefficient of a sorted, non-negative sample: twice the area between
  // the Lorenz curve and the diagonal. No closed-form standard error exists,
  // which is exactly why the bootstrap earns its place.
  function gini(sorted) {
    const n = sorted.length;
    let num = 0, tot = 0;
    for (let i = 0; i < n; i++) { num += (i + 1) * sorted[i]; tot += sorted[i]; }
    return tot > 0 ? (2 * num) / (n * tot) - (n + 1) / n : 0;
  }

  // ------------------------------------------------------- regression + series

  // Simple OLS of y on x, with both the classical and the heteroskedasticity-
  // robust (HC1) standard error. Deck 08 leans on the pair: the estimate and its
  // standard error fail independently, and the standard error fails first.
  function ols(x, y) {
    const n = x.length;
    const mx = mean(x), my = mean(y);
    let sxx = 0, sxy = 0;
    for (let i = 0; i < n; i++) { const dx = x[i] - mx; sxx += dx * dx; sxy += dx * (y[i] - my); }
    const b = sxx > 0 ? sxy / sxx : 0;
    const a = my - b * mx;
    let rss = 0, meat = 0, tss = 0;
    for (let i = 0; i < n; i++) {
      const e = y[i] - a - b * x[i], dx = x[i] - mx;
      rss += e * e; meat += dx * dx * e * e; tss += (y[i] - my) * (y[i] - my);
    }
    const se = sxx > 0 && n > 2 ? Math.sqrt(rss / (n - 2) / sxx) : NaN;
    const seR = sxx > 0 && n > 2
      ? Math.sqrt(meat * n / (n - 2)) / sxx                    // HC1
      : NaN;
    return { a, b, se, seR, t: b / se, tR: b / seR,
             r2: tss > 0 ? 1 - rss / tss : 0, n };
  }

  const corr = (x, y) => {
    const mx = mean(x), my = mean(y);
    let sxy = 0, sxx = 0, syy = 0;
    for (let i = 0; i < x.length; i++) {
      const dx = x[i] - mx, dy = y[i] - my;
      sxy += dx * dy; sxx += dx * dx; syy += dy * dy;
    }
    return sxx > 0 && syy > 0 ? sxy / Math.sqrt(sxx * syy) : 0;
  };

  // Stationary AR(1), started from its stationary distribution so the first
  // observations are not systematically closer to zero than the rest.
  function ar1(n, rho, u, sd = 1) {
    const out = new Float64Array(n);
    let v = u.normal() * sd / Math.sqrt(Math.max(1 - rho * rho, 1e-9));
    for (let i = 0; i < n; i++) { v = rho * v + sd * u.normal(); out[i] = v; }
    return out;
  }

  const randomWalk = (n, u, sd = 1) => {
    const out = new Float64Array(n);
    let v = 0;
    for (let i = 0; i < n; i++) { v += sd * u.normal(); out[i] = v; }
    return out;
  };

  // G clusters of m units with intra-cluster correlation rho: a shared group
  // draw plus an idiosyncratic one, weighted so the total variance stays 1.
  function clustered(G, m, rho, u) {
    const out = new Float64Array(G * m);
    const sg = Math.sqrt(Math.max(rho, 0)), si = Math.sqrt(Math.max(1 - rho, 0));
    for (let g = 0; g < G; g++) {
      const shared = sg * u.normal();
      for (let i = 0; i < m; i++) out[g * m + i] = shared + si * u.normal();
    }
    return out;
  }

  // ---------------------------------------------------------------- formatting

  const fmt = (v, d = 2) => (Number.isFinite(v) ? v.toFixed(d) : "—");
  const signed = (v, d = 2) => (Number.isFinite(v) ? (v >= 0 ? "+" : "") + v.toFixed(d) : "—");

  function smart(v) {
    if (!Number.isFinite(v)) return "—";
    const a = Math.abs(v);
    if (a >= 1e5 || (a > 0 && a < 1e-3)) return v.toExponential(1);
    return v.toFixed(a >= 100 ? 0 : a >= 10 ? 1 : 2);
  }

  function nice(lo, hi, n) {
    const span = (hi - lo) || Math.abs(hi) || 1;
    const raw = span / n;
    const mag = Math.pow(10, Math.floor(Math.log10(raw)));
    const norm = raw / mag;
    const step = mag * (norm < 1.5 ? 1 : norm < 3 ? 2 : norm < 7 ? 5 : 10);
    const out = [];
    for (let v = Math.ceil(lo / step - 1e-9) * step; v <= hi + 1e-9 * span; v += step)
      out.push(Math.abs(v) < step * 1e-9 ? 0 : v);
    return { ticks: out, step };
  }

  function logTicks(lo, hi) {
    const decades = Math.log10(hi / lo);
    const mults = decades > 2.2 ? [1] : decades > 1.1 ? [1, 3] : [1, 2, 5];
    const out = [];
    for (let e = Math.floor(Math.log10(lo)); e <= Math.ceil(Math.log10(hi)); e++)
      for (const k of mults) {
        const v = k * Math.pow(10, e);
        if (v >= lo * 0.9999 && v <= hi * 1.0001) out.push(v);
      }
    return out;
  }

  const decimalsFor = (step) =>
    Math.max(0, Math.min(6, -Math.floor(Math.log10(Math.abs(step)) + 1e-9)));

  // -------------------------------------------------------------------- charts

  function chart(W, H, extra) {
    const svg = el("svg", Object.assign({
      viewBox: `0 0 ${W} ${H}`,
      style: `width:100%;height:auto;display:block;font-family:inherit;color:${C.ink}`
    }, extra || {}));
    svg.__W = W;
    svg.__H = H;
    return svg;
  }

  function panel(svg, o) {
    const L = o.left ?? 0, T = o.top ?? 0;
    const W = o.w ?? svg.__W, H = o.h ?? svg.__H;
    const m = Object.assign({ l: 74, r: 22, t: 52, b: 54 }, o.margin || {});
    const fs = o.fontSize ?? 15;
    const iw = W - m.l - m.r, ih = H - m.t - m.b;
    const px0 = L + m.l, py0 = T + m.t;

    let [xa, xb] = o.xdom, [ya, yb] = o.ydom;
    if (!(xb > xa)) xb = xa + 1;
    if (!(yb > ya)) yb = ya + 1;

    const lg = o.xlog ? Math.log : (v => v);
    const xA = lg(xa), xB = lg(xb);
    const x = v => px0 + (lg(v) - xA) / (xB - xA) * iw;
    const y = v => py0 + ih - (v - ya) / (yb - ya) * ih;

    const gAxis = el("g");
    const id = `ldkclip${++uid}`;
    const defs = el("defs");
    const cp = el("clipPath", { id });
    cp.append(el("rect", { x: px0 - 3, y: py0 - 8, width: iw + 6, height: ih + 11 }));
    defs.append(cp);
    const gc = el("g", { "clip-path": `url(#${id})` });
    const gt = el("g");
    svg.append(defs, gAxis, gc, gt);

    // ticks
    let xt, xfmt, yt, yfmt;
    if (o.xticks) {
      xt = o.xticks;
      xfmt = o.xfmt || (v => String(v));
    } else if (o.xlog) {
      xt = logTicks(xa, xb);
      xfmt = o.xfmt || (v => v >= 1 ? String(Math.round(v)) : String(v));
    } else {
      const n = nice(xa, xb, o.nx ?? 6);
      xt = n.ticks;
      xfmt = o.xfmt || (v => v.toFixed(decimalsFor(n.step)));
    }
    if (o.yticks) {
      yt = o.yticks;
      yfmt = o.yfmt || (v => String(v));
    } else {
      const n = nice(ya, yb, o.ny ?? 5);
      yt = n.ticks;
      yfmt = o.yfmt || (v => v.toFixed(decimalsFor(n.step)));
    }

    if (o.axes === false) { xt = []; yt = []; }   // bare drawing surface

    if (o.axes !== false) {
    gAxis.append(el("line", { x1: px0, y1: py0 + ih, x2: px0 + iw, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    gAxis.append(el("line", { x1: px0, y1: py0, x2: px0, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    }
    for (const t of xt) {
      const px = x(t);
      if (px < px0 - 0.5 || px > px0 + iw + 0.5) continue;
      gAxis.append(el("line", { x1: px, y1: py0 + ih, x2: px, y2: py0 + ih + 5,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px, y: py0 + ih + 6 + fs, "text-anchor": "middle",
                                "font-size": fs, fill: C.grey }, xfmt(t)));
    }
    for (const t of yt) {
      const py = y(t);
      if (py < py0 - 0.5 || py > py0 + ih + 0.5) continue;
      gAxis.append(el("line", { x1: px0 - 5, y1: py, x2: px0, y2: py,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px0 - 10, y: py + fs * 0.35, "text-anchor": "end",
                                "font-size": fs, fill: C.grey }, yfmt(t)));
    }
    if (o.xlab)
      gAxis.append(el("text", { x: px0 + iw / 2, y: T + H - 8, "text-anchor": "middle",
                                "font-size": fs + 2, fill: C.ink }, o.xlab));
    if (o.ylab) {
      const ty = py0 + ih / 2, tx = L + 20;
      gAxis.append(el("text", { x: tx, y: ty, "text-anchor": "middle", "font-size": fs + 2,
                                fill: C.ink, transform: `rotate(-90 ${tx} ${ty})` }, o.ylab));
    }

    const P = {
      x, y, svg, gAxis, gc, gt, fs,
      xdom: [xa, xb], ydom: [ya, yb],
      area: { x0: px0, y0: py0, w: iw, h: ih },

      add(node, top) { (top ? gt : gc).append(node); return node; },

      path(pts, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        return P.add(el("path", Object.assign(
          { d, fill: "none", "stroke-linejoin": "round", "stroke-linecap": "round" }, a)));
      },

      fill(pts, base, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        d += "L" + x(pts[pts.length - 1][0]).toFixed(2) + "," + y(base).toFixed(2);
        d += "L" + x(pts[0][0]).toFixed(2) + "," + y(base).toFixed(2) + "Z";
        return P.add(el("path", Object.assign({ d, stroke: "none" }, a)));
      },

      seg(x1, y1, x2, y2, a) {
        return P.add(el("line", Object.assign(
          { x1: x(x1), y1: y(y1), x2: x(x2), y2: y(y2) }, a)));
      },
      vline(v, a) { return P.seg(v, ya, v, yb, a); },
      hline(v, a) { return P.seg(xa, v, xb, v, a); },

      bar(x1, x2, y1, y2, a) {
        const X1 = x(x1), X2 = x(x2), Y1 = y(y1), Y2 = y(y2);
        return P.add(el("rect", Object.assign({
          x: Math.min(X1, X2), y: Math.min(Y1, Y2),
          width: Math.max(Math.abs(X2 - X1), 0.5),
          height: Math.max(Math.abs(Y2 - Y1), 0)
        }, a)));
      },

      dot(vx, vy, r, a) {
        return P.add(el("circle", Object.assign({ cx: x(vx), cy: y(vy), r }, a)));
      },

      // mtext-style annotation above the plot region
      header(s, a = {}) {
        const adj = a.adj ?? 0;
        const py = py0 - 14 - (a.line ?? 0) * (fs + 7);
        return P.add(el("text", {
          x: adj === 1 ? px0 + iw : adj === 0.5 ? px0 + iw / 2 : px0,
          y: py,
          "text-anchor": adj === 1 ? "end" : adj === 0.5 ? "middle" : "start",
          "font-size": a.size ?? fs + 2,
          "font-weight": a.weight ?? "bold",
          fill: a.fill ?? C.ink
        }, s), true);
      },

      legend(items, opt = {}) {
        const size = opt.size ?? fs;
        const pad = opt.pad ?? 8;
        const rowH = size + 8, swW = 30;
        const g = el("g");
        let maxw = 0;
        items.forEach((it, i) => {
          const ty = i * rowH + size;
          if (it.swatch)
            g.append(el("rect", { x: 0, y: ty - size * 0.82, width: swW,
                                  height: size * 0.85, fill: it.swatch, stroke: "none" }));
          else
            g.append(el("line", { x1: 0, y1: ty - size * 0.35, x2: swW, y2: ty - size * 0.35,
                                  stroke: it.color, "stroke-width": it.lw ?? 3.5,
                                  "stroke-dasharray": it.dash ?? null,
                                  "stroke-linecap": "round" }));
          g.append(el("text", { x: swW + 9, y: ty, "font-size": size,
                                fill: opt.fill ?? C.ink }, it.label));
          maxw = Math.max(maxw, swW + 9 + it.label.length * size * 0.54);
        });
        // opt-in backing panel, for legends that must sit over the data
        if (opt.bg) {
          const q = 8;
          g.insertBefore(el("rect", {
            x: -q, y: -q + 3, width: maxw + 2 * q,
            height: items.length * rowH + 2 * q - 6,
            fill: "white", "fill-opacity": opt.bg === true ? 0.85 : opt.bg, rx: 5
          }), g.firstChild);
        }
        const corner = opt.corner ?? "topright";
        const gx = corner.includes("right") ? px0 + iw - maxw - pad : px0 + pad;
        const gy = corner.includes("top") ? py0 + pad
                                          : py0 + ih - items.length * rowH - pad;
        g.setAttribute("transform", `translate(${gx},${gy})`);
        gt.append(g);
        return g;
      }
    };
    return P;
  }

  // Convenience: one chart with exactly one panel filling it.
  function single(W, H, o) {
    const svg = chart(W, H);
    return { svg, p: panel(svg, o) };
  }

  return {
    C, el, alpha, chart, panel, single,
    lgamma, lchoose, erfc, dnorm, pnorm, qnorm, dbinom, dpois, dexp, pexp,
    dunif, punif, dlnorm, dchisq1, cdfFromPmf,
    gammaP, betai, dbeta, qbeta, dchisq, pchisq, qchisq, dt, pt, qt, invert,
    gini,
    seq, cumsum, sum, mean, sd, quantileSorted, hist, rng, sampleWeighted,
    sampleReplace, ols, corr, ar1, randomWalk, clustered,
    fmt, signed, smart, nice, logTicks
  };
}

An econometric model links an outcome $Y$ to explanatory variables $X$ and to everything else we do not observe, $u$:

\[Y = f(X) + u\]

Most of the time we model the conditional mean $E[Y \mid X]$ — the average outcome among people who share the same $X$.

Describe

How do wages differ across education levels?

Predict

How likely is this applicant to default?

Explain

What does one more year of school do to a wage?

The same regression serves all three — but the third needs the most assumptions.

Three Outcomes, Three Shapes

The type of outcome decides which model fits. This deck walks through them in order.

Linear model for the first, logit / probit for the second, and a family of specialised models for the rest.

1 · The Linear Model

Does Education Pay?

Data: CPS1985 from the AER package — 534 US workers.

Each grey dot is one worker
Each pink dot is the average wage at that education level: $E[\text{wage} \mid \text{educ}]$
Dot size = number of workers

The idea of regression

Summarise the pink dots with a simple formula — a line.

Why Log Wages?

Statistical reason

Wages are skewed; the log pulls in the tail, so a few extreme earners do not dominate the fit.

Economic reason

Education plausibly raises wages by a percentage, not by a fixed dollar amount. In logs, coefficients read as % changes.

The Linear Model

\[\log(\text{wage}_i) = \underbrace{\beta_0 + \beta_1\,\text{educ}_i}_{E[Y \mid X]:\ \text{systematic part}} \; + \underbrace{u_i}_{\text{everything else}}\]

The pink line is $E[Y \mid X]$: where the average sits at each education level
The blue curves are the spread of $u$ around it: many people above, many below

The key assumption

$E[u \mid \text{educ}] = 0$: whatever is left in $u$ averages to zero at every education level. Then the line really is the conditional mean.

Least Squares, Live

Choose the line yourself. Least squares picks the one with the smallest sum of squared residuals (SSR).

viewof ls_b1 = Inputs.range([-0.05, 0.2], {step: 0.001, value: 0.01, label: "slope b₁"})
viewof ls_shift = Inputs.range([-0.6, 0.6], {step: 0.01, value: 0.35, label: "shift up / down"})
viewof ls_show = Inputs.toggle({label: "show the OLS line", value: false})
ls_jit = { const u = ldk.rng(3); return ls_x.map(() => (u() - 0.5) * 0.4); }
ls_ymax = {
  // fixed y-range: the worst SSR the sliders can reach, so the axis never jumps
  const mx = ldk.mean(ls_x), my = ldk.mean(ls_y);
  let w = 0;
  for (const b of [-0.05, 0.2]) for (const sh of [-0.6, 0.6]) {
    let s = 0;
    for (let i = 0; i < ls_x.length; i++) { const e = ls_y[i] - (my - b * mx + sh) - b * ls_x[i]; s += e * e; }
    w = Math.max(w, s);
  }
  return Math.ceil(w / 10) * 10;
}

60 workers drawn from CPS1985. Each vertical segment is a residual $y_i - \hat y_i$: pink above the line, blue below.

The orange dot is the point of means $(\bar x, \bar y)$. The least-squares line always passes through it, so shift = 0 is always part of the answer.

{
  const C = ldk.C;
  const x = ls_x, y = ls_y, n = x.length;
  const mx = ldk.mean(x), my = ldk.mean(y);
  const fit = ldk.ols(x, y);
  const ssr = (a, b) => {
    let s = 0;
    for (let i = 0; i < n; i++) { const e = y[i] - a - b * x[i]; s += e * e; }
    return s;
  };
  const b1 = ls_b1, b0 = my - b1 * mx + ls_shift;
  const S = ssr(b0, b1), Smin = ssr(fit.a, fit.b);

  const svg = ldk.chart(900, 480);

  // ---- left: the data, the line, the residuals
  const pl = ldk.panel(svg, {
    left: 0, w: 520, xdom: [5, 19], ydom: [0.3, 3.9],
    xlab: "years of education", ylab: "log(wage)",
    margin: { l: 62, r: 14, t: 74, b: 50 }, nx: 7
  });
  for (let i = 0; i < n; i++) {
    const xi = x[i] + ls_jit[i], yh = b0 + b1 * x[i];
    pl.seg(xi, y[i], xi, yh, { stroke: ldk.alpha(y[i] > yh ? C.pink : C.blue, 0.75),
                               "stroke-width": 2.2 });
  }
  for (let i = 0; i < n; i++)
    pl.dot(x[i] + ls_jit[i], y[i], 4.2, { fill: ldk.alpha(C.ink, 0.75) });
  pl.path([[5, b0 + b1 * 5], [19, b0 + b1 * 19]], { stroke: C.ink, "stroke-width": 3.5 });
  if (ls_show)
    pl.path([[5, fit.a + fit.b * 5], [19, fit.a + fit.b * 19]],
            { stroke: C.green, "stroke-width": 3.5, "stroke-dasharray": "10,6" });
  pl.dot(mx, my, 8, { fill: C.orange, stroke: "white", "stroke-width": 2 });
  pl.header(`your line:  ŷ = ${ldk.fmt(b0, 2)} + ${ldk.fmt(b1, 3)} · educ`, { line: 1, size: 17 });
  if (ls_show)
    pl.header(`OLS line:  ŷ = ${ldk.fmt(fit.a, 2)} + ${ldk.fmt(fit.b, 3)} · educ`,
              { size: 15, fill: C.green, weight: "normal" });

  // ---- right: SSR as a function of the slope
  const grid = ldk.seq(-0.05, 0.2, 101);
  const curve = grid.map(b => [b, ssr(my - b * mx + ls_shift, b)]);
  const best  = grid.map(b => [b, ssr(my - b * mx, b)]);
  const pr = ldk.panel(svg, {
    left: 520, w: 380, xdom: [-0.05, 0.2], ydom: [0, ls_ymax],
    xlab: "slope b₁", ylab: "SSR",
    margin: { l: 58, r: 16, t: 74, b: 50 }, nx: 5
  });
  pr.path(best, { stroke: ldk.alpha(C.grey, 0.5), "stroke-width": 2, "stroke-dasharray": "5,5" });
  pr.path(curve, { stroke: C.ink, "stroke-width": 3 });
  if (ls_show) pr.vline(fit.b, { stroke: C.green, "stroke-width": 2.5, "stroke-dasharray": "8,5" });
  pr.dot(b1, S, 8, { fill: C.pink, stroke: "white", "stroke-width": 2 });
  pr.header(`SSR = ${ldk.fmt(S, 2)}`, { line: 1, size: 18, fill: S - Smin < 0.05 ? C.green : C.pink });
  pr.header(`smallest possible: ${ldk.fmt(Smin, 2)}`, { size: 15, fill: C.grey, weight: "normal" });
  return svg;
}

Running It in R

m1 <- lm(log(wage) ~ education, data = CPS1985)
summary(m1)

Call:
lm(formula = log(wage) ~ education, data = CPS1985)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.98099 -0.37155  0.03391  0.34975  1.66098 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.059890   0.107432   9.866   <2e-16 ***
education   0.076759   0.008091   9.487   <2e-16 *** 123
---

Residual standard error: 0.4885 on 532 degrees of freedom 4
Multiple R-squared:  0.1447,    Adjusted R-squared:  0.1431  5
F-statistic: 90.01 on 1 and 532 DF,  p-value: < 2.2e-16

1 Estimate. One more year of education goes with $0.077$ higher log wage — about +8% ($e^{0.077}-1$).

2 Std. Error. How much the estimate would move from sample to sample.

3 t value $=$ estimate / SE, and its p-value for $H_0: \beta_1 = 0$. Here: no doubt that $\beta_1 \neq 0$.

4 Residual SE. The typical size of $u$: about $0.49$ log points.

5 R². Education accounts for 14% of the variation in log wages. The F-statistic tests all slopes $=0$ at once.

Standard Errors Are Sampling Variability

confint(m1)

                 2.5 %     97.5 %
(Intercept) 0.84884618 1.27093355
education   0.06086475 0.09265238

The Std. Error is the width of that histogram — no resampling needed, R computes it from one sample. The 95% interval is roughly estimate ± 2 SE: a return of 6% to 9% per year.

Reading Coefficients: Units Matter

Model	Formula in R	Interpretation of $\beta_1$
level–level	`lm(wage ~ educ)`	+1 year $\Rightarrow$ wage changes by $\beta_1$ dollars
log–level	`lm(log(wage) ~ educ)`	+1 year $\Rightarrow$ wage changes by about $100\cdot\beta_1$ %
level–log	`lm(wage ~ log(sales))`	+1% sales $\Rightarrow$ wage changes by $\beta_1/100$ dollars
log–log	`lm(log(q) ~ log(p))`	+1% price $\Rightarrow$ quantity changes by $\beta_1$ % — an elasticity

Exact, not approximate

In log–level models the exact change is $e^{\beta_1} - 1$. For $\beta_1 = 0.077$: $7.98\%$. For small $\beta$ the two agree; for $\beta = -0.26$ they do not ($-22.7\%$).

Always state the units

“The coefficient is 0.077” means nothing on its own. “One more year of schooling goes with about 8% higher hourly wages” is an answer.

More Than One Regressor

m2 <- lm(log(wage) ~ education + experience + I(experience^2) +
           gender, data = CPS1985)
summary(m2)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      0.6007445  0.1194927   5.027 6.81e-07 ***
education        0.0912936  0.0080049  11.405  < 2e-16 *** 1
experience       0.0360522  0.0054352   6.633 8.14e-11 *** 2
I(experience^2) -0.0005412  0.0001197  -4.520 7.64e-06 *** 3
genderfemale    -0.2570355  0.0387066  -6.641 7.77e-11 *** 4

Residual standard error: 0.4442 on 529 degrees of freedom
Multiple R-squared:  0.2968,    Adjusted R-squared:  0.2915  5
F-statistic: 55.81 on 4 and 529 DF,  p-value: < 2.2e-16

Each coefficient is now an effect holding the other regressors fixed — ceteris paribus.

1 Return to education rises from 0.077 to 0.091. Why? Next slide.

2 3 Experience raises wages, but by less and less — the square lets the profile bend.

4 Women earn about 23% less ($e^{-0.257}-1$) at the same education and experience.

5 Adjusted R² penalises extra regressors; it doubled.

Why Did the Education Coefficient Move?

Leaving experience out of the model puts it into $u$. It is not harmless there:

experience raises wages ($+$)
experience is lower for the more educated ($-$)

Omitted variable bias

bias $=$ (effect of the omitted variable on $Y$) $\times$ (how it moves with $X$)

Here $(+)\times(-)$: the short regression understates the return to education.

Curves Are Still “Linear” Models

\[\ldots + \beta_2\,\text{exper} + \beta_3\,\text{exper}^2\]

“Linear” means linear in the coefficients $\beta$, not in $X$. Squares, logs and interactions all work in lm().

One more year of experience now adds $\beta_2 + 2\beta_3\,\text{exper}$: a lot early in a career, zero at the peak $-\hat\beta_2 / (2\hat\beta_3) \approx 33$ years.

In R

Write I(experience^2) — inside a formula, a bare ^ means something else.

Dummy Variables Shift the Line

gender is a factor; R turns it into a 0/1 dummy genderfemale automatically. The first level (male) is the reference group.

The dummy moves the intercept:

men: $\beta_0 + \ldots$
women: $(\beta_0 - 0.257) + \ldots$

Parallel by construction

The gap is the same at every education level — because we assumed it, not because the data said so. To let it differ, interact.

Interactions Let the Slopes Differ

m3 <- lm(log(wage) ~ education * gender + experience +
           I(experience^2), data = CPS1985)

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.7896897  0.1402019   5.633 2.89e-08 ***
education               0.0762191  0.0099290   7.676 7.97e-14 *** 1
genderfemale           -0.7539256  0.1992145  -3.784 0.000172 *** 2
experience              0.0369215  0.0054182   6.814 2.59e-11 ***
I(experience^2)        -0.0005585  0.0001193  -4.681 3.64e-06 ***
education:genderfemale  0.0381424  0.0150037   2.542 0.011300 *   3

1 Return to education for men: 7.6%.

2 Gender gap at zero education — an extrapolation, do not read it alone.

3 Extra return for women: $0.076 + 0.038 = 0.114$ per year.

The gender gap narrows with education.

Testing Several Coefficients at Once

The $t$-test in summary() checks one coefficient. To ask “do experience and gender matter at all?” compare the two models with an F-test:

anova(m1, m2)

Analysis of Variance Table

Model 1: log(wage) ~ education
Model 2: log(wage) ~ education + experience + I(experience^2) + gender
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    532 126.97                                   1
2    529 104.39  3    22.576 38.134 < 2.2e-16 *** 23

1 SSR of the restricted model — experience and gender set to zero.

2 SSR of the full model: 3 extra coefficients, 22.6 less squared error.

3 Is that drop larger than chance would produce? $F = 38$, $p \approx 0$: yes.

The logic from the testing deck

Same recipe as before: a test statistic, its distribution under $H_0$, a p-value. car::linearHypothesis() tests any other linear restriction, e.g. $\beta_\text{exper} = \beta_\text{educ}$.

The Assumptions, in Plain Words

Assumption	In words	What fails without it
Linear in parameters	$Y$ is a sum of $\beta \times$ (something built from $X$)	The fitted line is the wrong shape
Random sample	Observations are independent draws from the population	Standard errors — see clustering
No perfect collinearity	No regressor is an exact combination of others	R cannot separate them: `NA` coefficient
Exogeneity $E[u \mid X] = 0$	The omitted stuff is unrelated to $X$	The estimates themselves are biased
Homoskedasticity	The spread of $u$ is the same for every $X$	Only the usual standard errors
Normal errors	$u$ is normally distributed	Nothing, once $n$ is moderate (CLT)

Two very different kinds of failure

The first four give you the right number. The last two only give you the right uncertainty — and are easy to fix. Exogeneity is the one that keeps econometricians busy.

Heteroskedasticity: Look at the Residuals

m0 <- lm(wage ~ education + experience +
           I(experience^2) + gender, data = CPS1985)
round(cbind(usual  = sqrt(diag(vcov(m0))),
            robust = sqrt(diag(vcovHC(m0, type = "HC1")))), 4)

                 usual robust
(Intercept)     1.1888 1.2493
education       0.0796 0.0880
experience      0.0541 0.0631
I(experience^2) 0.0012 0.0013
genderfemale    0.3851 0.3936

The estimates do not change — only the standard errors
Robust (“sandwich”, HC1) SEs are valid whether or not the spread is constant
In a table of results: coeftest(m0, vcov = vcovHC(m0, type = "HC1"))

Rule of thumb

Use robust standard errors by default.

Regression Is Not Yet Causation

Ability raises both schooling and wages. It sits in $u$, so $E[u \mid \text{educ}] \neq 0$ and $\hat\beta_1$ mixes the effect of school with the effect of ability.

What $\hat\beta_1 = 0.077$ is

A well-estimated conditional mean difference: how wages differ between people with different schooling.

What it is not (yet)

The effect of sending someone to school for one more year. That needs a design — panel data, IV, DiD in part 3.

2 · Binary Choice

When the Outcome Is Yes or No

Labour-force participation, loan default, buying a car, migrating, voting: many outcomes are 0 or 1.

Data: HMDA — 2,380 Boston mortgage applications, 12% denied. The key regressor is the ratio of monthly debt payments to income.

The conditional mean is a probability

For $Y \in \{0,1\}$: \[E[Y \mid X] = P(Y = 1 \mid X)\] So a model for $E[Y\mid X]$ is a model for the probability of a yes.

The Linear Probability Model

Just run OLS on the 0/1 outcome — a linear probability model (LPM):

lpm <- lm(I(deny == "yes") ~ pirat, data = HMDA)
coeftest(lpm, vcov = vcovHC(lpm, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -0.079910   0.031967 -2.4998   0.01249 *   3
pirat        0.603535   0.098483  6.1283 1.036e-09 *** 12

1 Coefficients are changes in probability: $+0.1$ in P/I raises the denial probability by $0.06$ — 6 percentage points.

2 Robust SEs are mandatory here: with a 0/1 outcome, $\text{Var}(u \mid X) = p(1-p)$ changes with $X$ by construction.

3 The intercept is a negative probability. That is a warning sign.

Why people like it

Easy to fit, easy to read, and its slopes are usually close to the average effects of the fancier models. A sensible first look.

Where the Straight Line Breaks

Impossible predictions. At P/I = 3 (it happens) the LPM predicts a denial probability of 1.73
Constant effects. +0.1 in P/I adds 6 points whether you start at 1% or at 95% — but near 0 and 1 a probability cannot keep moving at the same pace
Built-in heteroskedasticity — fixed by robust SEs

What we want

A curve that stays inside $[0, 1]$ and flattens at both ends.

Bend the Line: The S-Curve

\[P(Y = 1 \mid X) = G(\beta_0 + \beta_1 X), \qquad G:\ \mathbb{R} \to [0, 1] \text{ increasing}\]

Two choices of G

Logit uses the logistic CDF, probit the normal CDF. Same shape; the logistic has fatter tails.

Different scales

The logistic is more spread out, so logit coefficients are about 1.6× probit coefficients. Only the probabilities are comparable.

Where the S-Curve Comes From, Live

Each applicant has a latent “denial score” $y^* = \beta_0 + \beta_1 x + \varepsilon$. The bank denies if $y^* > 0$. So $P(\text{deny}) = P(\varepsilon > -\beta_0 - \beta_1 x) = G(\beta_0 + \beta_1 x)$.

viewof lv_x = Inputs.range([0, 1.2], {step: 0.01, value: 0.33, label: "P/I ratio x"})
viewof lv_link = Inputs.radio(["logit", "probit"], {label: "error ε", value: "logit"})

Coefficients are the ones fitted to HMDA on the next slides.

Left: the distribution of $y^*$ for applicants with this $x$. The shaded area beyond 0 is the probability of denial.

Right: that area, traced out over $x$. The dashed tangent is the marginal effect at this $x$: density $\times\ \beta_1$.

{
  const C = ldk.C;
  const logit = lv_link === "logit";
  const [b0, b1] = logit ? lv_logit : lv_probit;
  const dens = logit ? (z => { const e = Math.exp(-Math.abs(z)); return e / ((1 + e) * (1 + e)); })
                     : (z => ldk.dnorm(z, 0, 1));
  const cdf  = logit ? (z => 1 / (1 + Math.exp(-z))) : (z => ldk.pnorm(z, 0, 1));
  const idx = b0 + b1 * lv_x;
  const P = cdf(idx), me = dens(idx) * b1;

  const svg = ldk.chart(940, 450);

  // ---- left: latent density
  const lo = logit ? -9 : -5, hi = logit ? 6 : 4;
  const pl = ldk.panel(svg, {
    left: 0, w: 470, xdom: [lo, hi], ydom: [0, logit ? 0.3 : 0.45],
    xlab: "latent score y*", ylab: "density",
    margin: { l: 62, r: 16, t: 70, b: 50 }, nx: 6
  });
  const g = ldk.seq(lo, hi, 300).map(v => [v, dens(v - idx)]);
  pl.fill(g.filter(p => p[0] >= 0), 0, { fill: ldk.alpha(C.pink, 0.45) });
  pl.fill(g.filter(p => p[0] <= 0), 0, { fill: ldk.alpha(C.grey, 0.15) });
  pl.path(g, { stroke: C.ink, "stroke-width": 2.6 });
  pl.vline(0, { stroke: C.pink, "stroke-width": 2.5 });
  pl.vline(idx, { stroke: C.blue, "stroke-width": 1.8, "stroke-dasharray": "6,5" });
  pl.header(`index β₀ + β₁x = ${ldk.fmt(idx, 2)}`, { line: 1, size: 17, fill: C.blue });
  pl.header("shaded: y* > 0  →  denied", { size: 15, fill: C.pink, weight: "normal" });

  // ---- right: the S-curve with data and tangent
  const pr = ldk.panel(svg, {
    left: 470, w: 470, xdom: [0, 1.2], ydom: [0, 1],
    xlab: "P/I ratio x", ylab: "P(deny)",
    margin: { l: 62, r: 16, t: 70, b: 50 }, nx: 6
  });
  for (let i = 0; i < lv_bin_x.length; i++)
    pr.dot(lv_bin_x[i], lv_bin_y[i], 5.5, { fill: ldk.alpha(C.pink, 0.55) });
  pr.path(ldk.seq(0, 1.2, 200).map(v => [v, cdf(b0 + b1 * v)]),
          { stroke: C.ink, "stroke-width": 3 });
  pr.path([[lv_x - 0.25, P - 0.25 * me], [lv_x + 0.25, P + 0.25 * me]],
          { stroke: C.orange, "stroke-width": 3, "stroke-dasharray": "8,5" });
  pr.seg(lv_x, 0, lv_x, P, { stroke: ldk.alpha(C.blue, 0.6), "stroke-width": 1.5,
                            "stroke-dasharray": "4,4" });
  pr.dot(lv_x, P, 8, { fill: C.pink, stroke: "white", "stroke-width": 2 });
  pr.header(`P(deny) = ${ldk.fmt(P, 3)}`, { line: 1, size: 17, fill: C.pink });
  pr.header(`marginal effect = ${ldk.fmt(me, 3)} per unit of x`,
            { size: 15, fill: C.orange, weight: "normal" });
  return svg;
}

Estimation: Maximum Likelihood

There is no “residual” to square. Instead ask: which $\beta$ makes the observed denials and approvals most probable?

a denied applicant contributes $p_i = G(\beta_0 + \beta_1 x_i)$
an approved one contributes $1 - p_i$

\[\ell(\beta) = \sum_{i} \Big[ y_i \log p_i + (1 - y_i) \log(1 - p_i) \Big]\]

R climbs this function numerically — no formula like OLS.

From the estimation deck

This is the MLE recipe exactly. Every likelihood is a loss; OLS was the special case with normal errors.

Logit in R

lg <- glm(deny ~ pirat, family = binomial(link = "logit"),
          data = HMDA)
summary(lg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -4.0284     0.2686 -14.999  < 2e-16 ***
pirat         5.8845     0.7336   8.021 1.05e-15 *** 12

    Null deviance: 1744.2  on 2379  degrees of freedom
Residual deviance: 1660.2  on 2378  degrees of freedom 3
AIC: 1664.2 4

Number of Fisher Scoring iterations: 5 5

1 The coefficient is on the index $\beta_0 + \beta_1 x$, i.e. on the log-odds. Its sign and significance are readable; its size is not a probability.

2 z value, not $t$: MLE standard errors rely on large samples.

3 Deviance $= -2 \times$ log-likelihood. The drop from the null deviance (1744) is a likelihood-ratio test of the slope.

4 AIC: deviance + a penalty for parameters; smaller is better when comparing models.

5 Steps of the numerical optimiser.

Making the Coefficient Speak

Odds ratios — logit only:

exp(coef(lg))

 (Intercept)        pirat 
  0.01780222 359.42220584

exp(0.1 * coef(lg)["pirat"])   # per 0.1 of P/I

   pirat 
1.801194

Raising P/I by 0.1 multiplies the odds of denial, $p/(1-p)$, by 1.8. Common in medicine; many find odds hard to think in.

Predicted probabilities — always work:

predict(lg, newdata = data.frame(pirat = c(0.3, 0.5)),
        type = "response")

         1          2 
0.09422692 0.25233766

Moving from P/I 0.3 to 0.5 raises the denial probability from 9% to 25%. type = "response" returns $G(\hat\beta_0 + \hat\beta_1 x)$; the default "link" returns the index.

What you can read straight from the logit table

Sign (direction), significance (whether), and ratios of two coefficients (relative importance). Not: “how many percentage points”.

Logit, Probit, LPM: Same Story

rbind(logit  = coef(lg),
      probit = coef(pb),
      LPM    = coef(lpm))

       (Intercept)     pirat
logit  -4.02843220 5.8844978
probit -2.19414501 2.9678666
LPM    -0.07990965 0.6035349

Coefficients differ a lot: 5.9 vs 3.0 vs 0.6 — different scales, not different findings
Predicted probabilities are nearly identical where the data are
Logit vs probit is rarely a substantive choice. Pick one and report probabilities or marginal effects

Marginal Effects Vary Along the Curve

In logit/probit the effect of $x$ on the probability is

\[\frac{\partial P}{\partial x} = g(\beta_0 + \beta_1 x)\cdot\beta_1\]

— it depends on where you stand. Two ways to boil it down to one number:

Two one-number summaries

AME (average marginal effect): the effect for every person, then averaged — the standard choice.

MEM (effect at the mean): plug in the average person, who may not exist.

Average Marginal Effects in R

eta <- predict(lg, type = "link")      # index, per person
mean(dlogis(eta)) * coef(lg)["pirat"]  # logit AME

    pirat 
0.5907417

eta_p <- predict(pb, type = "link")
mean(dnorm(eta_p)) * coef(pb)["pirat"] # probit AME

    pirat 
0.5664755

coef(lpm)["pirat"]                     # LPM slope

    pirat 
0.6035349

All three models now say the same: +0.1 in P/I ≈ +6 percentage points in the denial probability, on average.

In practice

The package marginaleffects does this for any model, with standard errors:

avg_slopes(lg)

For a dummy it averages the change in predicted probability from switching 0 → 1.

A Richer Model

lg3 <- glm(deny ~ pirat + lvrat + phist + insurance + afam,
           family = binomial(link = "logit"), data = HMDA)
summary(lg3)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -5.7309     0.4705 -12.181  < 2e-16 ***
pirat          4.8496     0.7811   6.209 5.34e-10 *** 1
lvrat          1.9120     0.4815   3.971 7.16e-05 ***
phistyes       1.6702     0.1898   8.801  < 2e-16 *** 2
insuranceyes   4.3560     0.5450   7.993 1.32e-15 *** 3
afamyes        0.9622     0.1663   5.787 7.17e-09 *** 4

Controls: loan-to-value ratio lvrat, public bad-credit record phist, denied mortgage insurance insurance, and whether the applicant is African American afam.

1 P/I still matters, holding the rest fixed.

2 3 Bad credit and denied insurance are strong predictors — insurance denial almost settles it.

4 The afam coefficient is positive and clearly significant: a race gap after these controls — the question the original study asked.

From Log-Odds to Percentage Points

# helper in this deck's source; the same numbers
# as marginaleffects::avg_slopes(lg3)
ame(lg3)

                AME Std. Error z value Pr(>|z|)
pirat        0.3870     0.0623    6.21 5.28e-10
lvrat        0.1526     0.0387    3.95 7.89e-05
phistyes     0.2129     0.0333    6.40 1.59e-10
insuranceyes 0.7345     0.0679   10.82 2.75e-27
afamyes      0.0959     0.0201    4.78 1.71e-06 1

1 At the same P/I, loan-to-value, credit history and insurance status, African American applicants were on average 9.6 percentage points more likely to be denied (base rate 12%).

Predicted Probabilities Tell the Story

Fix the controls at a typical applicant and vary one thing at a time. predict(..., type = "response") on a grid of values; build bands on the link scale, then transform.

The gap widens as P/I rises: in logit and probit the effect of a dummy in percentage points depends on the other regressors, even without an interaction term.

Good practice

A plot of predicted probabilities for a few meaningful profiles usually communicates better than any table.

How Well Does It Fit?

# McFadden pseudo-R²
1 - lg3$deviance / lg3$null.deviance

[1] 0.2246673

table(predicted = fitted(lg3) > 0.5, actual = HMDA$deny)

         actual
predicted   no  yes
    FALSE 2073  212
    TRUE    22   73

Pseudo-R² compares log-likelihoods; values of 0.2–0.4 are already good. Not comparable to OLS R²
Accuracy is 90% — but saying “approve” to everyone scores 88%. With rare outcomes, look at the distributions, not the hit rate

Binary Choice: Practical Advice

Choosing a model

LPM with robust SEs: fine for average effects, easy with fixed effects and IV
Logit / probit: when predictions must stay in $[0,1]$ or you care about effects away from the average
Report AMEs or predicted probabilities, not raw coefficients

Traps

Perfect separation: a regressor predicts $Y$ perfectly; R warns fitted probabilities numerically 0 or 1 and the coefficient runs off to infinity
Do not compare raw coefficients across models or samples — the scale is not fixed
Interactions in logit/probit: the coefficient on the product is not the interaction effect on probabilities

3 · A Tour of Other Models

Match the Model to the Outcome

Outcome	Example	Model	R
Continuous	log wage	linear regression	`lm()`
Binary	mortgage denied	logit / probit	`glm(family = binomial)`
Ordered categories	health: poor … excellent	ordered logit / probit	`MASS::polr()`
Unordered categories	car, bus, bike	multinomial logit	`nnet::multinom()`
Count	doctor visits	Poisson / negative binomial	`glm(family = poisson)`, `MASS::glm.nb()`
Many zeros, then continuous	hours worked	Tobit	`AER::tobit()`
Observed only for a selected group	wages of workers	Heckman selection	`sampleSelection::selection()`
Time until an event	unemployment spell	duration / survival	`survival::coxph()`

And by data structure or design

The same unit observed many times $\to$ panel / fixed effects. A policy that hits some units at some date $\to$ difference-in-differences. An endogenous regressor with an outside source of variation $\to$ instrumental variables. Effects that differ across the distribution $\to$ quantile regression.

Ordered and Unordered Choices

Ordered: one latent scale, several cut-points

The binary latent model with more thresholds. $X$ shifts the distribution; the areas between cut-points are the probabilities. MASS::polr()

Unordered: one equation per alternative

Car, bus or bike: no natural order. Multinomial logit models each option relative to a base option. nnet::multinom()

Count Data: Poisson and Negative Binomial

Counts are $0, 1, 2, \ldots$ and skewed. Model the log of the mean: $E[Y\mid X] = e^{X\beta}$, so $e^{\beta}-1$ is a % change.

Poisson forces variance $=$ mean. Here the mean is 5.8 and the variance 45.7: overdispersion. The negative binomial adds a parameter for it — and matches the zeros.

nb <- glm.nb(visits ~ health + chronic + age +
               gender + insurance, data = NMES1988)
round(exp(coef(nb)[c("chronic", "insuranceyes")]), 3)

     chronic insuranceyes 
       1.216        1.388

Each chronic condition: +22% visits; insurance: +39%.

Corner Solutions: Tobit

A pile-up at zero: many women choose not to work. A straight line through this cloud is pulled toward the zeros.

Tobit: latent desired hours $y^* = X\beta + u$; we see $y = \max(0, y^*)$. Estimated by maximum likelihood, like probit for the zeros and OLS for the rest.

PSID1976$nwinc <- with(PSID1976,
                       (fincome - hours * wage) / 1000)
f  <- hours ~ nwinc + education + experience +
              age + youngkids
tb <- tobit(f, data = PSID1976)
round(cbind(OLS   = coef(lm(f, data = PSID1976)),
            Tobit = coef(tb)), 1)

               OLS  Tobit
(Intercept) 1306.9 1287.8
nwinc         -3.5   -8.5
education     32.7   86.1
experience    46.8   77.4
age          -30.3  -58.3
youngkids   -444.6 -919.9

Tobit coefficients refer to desired hours $y^*$ and are larger than OLS. For effects on actual hours, compute marginal effects — as with probit.

Sample Selection

We see wages only for people who work. The low-educated work mainly if their offer is unusually high, so the slope comes out too flat.

Heckman’s two steps

a probit for who works
the wage regression plus a correction term from step 1

Needs a variable that shifts working but not wages (e.g. young children).

Simulated data. In R: sampleSelection::selection().

Panel Data: Fixed Effects

Across states, higher beer taxes go with more deaths — because rural, southern states have both. That is omitted variable bias again, from things that differ between states but barely change over time.

Fixed effects

Give every state its own intercept. The slope is then estimated only from changes within a state over time — everything constant about a state drops out, observed or not.

The within-state slopes point down.

Fixed Effects in R

library(fixest)
pooled <- feols(frate ~ beertax, Fatalities, cluster = ~state)
fe_s   <- feols(frate ~ beertax | state, Fatalities, cluster = ~state)
fe_st  <- feols(frate ~ beertax | state + year, Fatalities,
                cluster = ~state)
etable(pooled, fe_s, fe_st)

                           pooled              fe_s             fe_st
Dependent Var.:             frate             frate             frate
                                                                     
Constant        1.853*** (0.1185)                                    
beertax         0.3646** (0.1197) -0.6559* (0.2919) -0.6400. (0.3571) 12
Fixed-Effects:  ----------------- ----------------- -----------------
state                          No               Yes               Yes
year                           No                No               Yes 3
_______________ _________________ _________________ _________________
S.E.: Clustered         by: state         by: state         by: state 4
Observations                  336               336               336
R2                        0.09336           0.90501           0.90893
Within R2                      --           0.04075           0.03606

1 Pooled OLS: taxes raise deaths?

2 State fixed effects: the sign flips. A $1 tax increase goes with 0.66 fewer deaths per 10,000.

3 Year effects also remove nationwide shocks (recessions, safety laws) — the two-way FE model.

4 Clustered SEs: the 7 years of one state are not independent draws.

Difference-in-Differences

A policy hits some units from some date on. Compare the change in the treated group with the change in the control group:

\[\hat\delta = (\bar y^{T}_\text{after} - \bar y^{T}_\text{before}) - (\bar y^{C}_\text{after} - \bar y^{C}_\text{before})\]

Parallel trends

Without the policy, both groups would have moved in parallel. Check the pre-period; that is where the credibility comes from.

In R: feols(y ~ treated:post | unit + year, data) — a two-way fixed-effects regression.

Instrumental Variables

Prices and quantities are set together by supply and demand, so price is endogenous in a demand equation.

An instrument $Z$ moves $X$ but has no other route to $Y$. Use only the part of $X$ that $Z$ moves:

regress $X$ on $Z$ (first stage) and keep the fitted $\hat X$
regress $Y$ on $\hat X$ (second stage)

This is two-stage least squares (2SLS).

Two conditions

Relevance — $Z$ really moves $X$: testable.

Exogeneity — $Z$ affects $Y$ only through $X$: an argument, not a test.

IV in R: Cigarette Demand

c95 <- subset(CigarettesSW, year == "1995")
c95 <- transform(c95, rprice = price / cpi, salestax = (taxs - tax) / cpi,
                      rincome = income / population / cpi)
iv <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + salestax,
            data = c95)
summary(iv, diagnostics = TRUE)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.4307     1.3584   6.943 1.24e-08 ***
log(rprice)   -1.1434     0.3595  -3.181  0.00266 **  1
log(rincome)   0.2145     0.2686   0.799  0.42867    

Diagnostic tests:
                 df1 df2 statistic  p-value    
Weak instruments   1  45    45.158 2.65e-08 *** 2
Wu-Hausman         1  44     1.102      0.3     3

Residual standard error: 0.1896 on 45 degrees of freedom
Multiple R-Squared: 0.4189, Adjusted R-squared: 0.3931 
Wald test: 6.534 on 2 and 45 DF,  p-value: 0.003227

After the |: the instruments — exogenous regressors (income) plus the excluded instrument (sales tax).

1 Price elasticity of demand: a 1% price rise cuts packs sold by 1.14%. OLS gives -1.41.

2 Weak-instrument test = first-stage F. Far above the rule-of-thumb 10: the tax is relevant.

3 Wu–Hausman: OLS and IV do not differ significantly here. With one instrument, exogeneity cannot be tested.

Quantile Regression

OLS models the mean. Quantile regression models a percentile of $Y$ given $X$ — e.g. the wage at the bottom 10% of workers with 12 years of school.

qr <- rq(log(wage) ~ education,
         tau = c(0.1, 0.5, 0.9), data = CPS1985)
round(coef(qr), 3)

            tau= 0.1 tau= 0.5 tau= 0.9
(Intercept)    0.717    1.019    1.453
education      0.056    0.082    0.093

Education pays more at the top: 5.6% per year at the 10th percentile, 9.3% at the 90th. It also spreads wages out. The asymmetric loss from the estimation deck, put to work.

Also Worth Knowing

Time series

ARIMA for one series, VAR for several that move together. Beware trends: two unrelated trending series look correlated.

Duration models

Time until an event — finding a job, a firm exiting. Handle spells still running at the end of the data (censoring). survival::coxph().

Regression discontinuity

Treatment switches on at a cut-off in a running variable (a test score, an age). Compare units just below and just above. rdrobust.

Machine learning meets econometrics

Lasso and random forests for prediction and for choosing controls; double machine learning for causal effects with many controls.

The common thread

Each is a conditional mean with a new shape (fitted by least squares or ML), or a design that makes $E[u \mid X] = 0$ believable.

Summary

Key Takeaways

A regression models the conditional mean $E[Y\mid X]$. Read every coefficient in units
summary() output: estimate, standard error, $t$ / $z$, p-value, fit. Use robust standard errors by default
“Holding other things fixed” only holds for what is in the model — omitted variables bias the rest
For 0/1 outcomes: LPM, logit or probit. Report marginal effects or predicted probabilities, not raw logit coefficients
Let the outcome choose the model — and the design decide whether the answer is causal

The one thing to remember

Every model in this deck answers: how does the distribution of $Y$ change with $X$?

The models differ in the shape they allow. Whether the answer is causal depends on $E[u\mid X]=0$ — and no choice of model can promise that.

Model	Formula in R	Interpretation of \(\beta_1\)
level–level	`lm(wage ~ educ)`	+1 year \(\Rightarrow\) wage changes by \(\beta_1\) dollars
log–level	`lm(log(wage) ~ educ)`	+1 year \(\Rightarrow\) wage changes by about \(100\cdot\beta_1\) %
level–log	`lm(wage ~ log(sales))`	+1% sales \(\Rightarrow\) wage changes by \(\beta_1/100\) dollars
log–log	`lm(log(q) ~ log(p))`	+1% price \(\Rightarrow\) quantity changes by \(\beta_1\) % — an elasticity

Assumption	In words	What fails without it
Linear in parameters	\(Y\) is a sum of \(\beta \times\) (something built from \(X\))	The fitted line is the wrong shape
Random sample	Observations are independent draws from the population	Standard errors — see clustering
No perfect collinearity	No regressor is an exact combination of others	R cannot separate them: `NA` coefficient
Exogeneity \(E[u \mid X] = 0\)	The omitted stuff is unrelated to \(X\)	The estimates themselves are biased
Homoskedasticity	The spread of \(u\) is the same for every \(X\)	Only the usual standard errors
Normal errors	\(u\) is normally distributed	Nothing, once \(n\) is moderate (CLT)