The big ideas of data analysis and modeling, with the fine print

Statement

For $y=X\beta+\varepsilon$ with $X\in\mathbb R^{n\times p}$, the fitted values are $\hat y=Hy$ with $H=X(X^\top X)^{-1}X^\top$. The leverage $h_{ii}=H_{ii}$ satisfies $1/n\le h_{ii}\le1$ and $\sum_i h_{ii}=p$. For simple regression $h_{ii}=\frac1n+\frac{(x_i-\bar x)^2}{S_{xx}}$. The residual variance is $\operatorname{Var}(e_i)=\sigma^2(1-h_{ii})$, and the influence of row $i$ on the coefficients is summarized by Cook's distance:

Why it's true

The Sherman–Morrison formula gives $(X_{(i)}^\top X_{(i)})^{-1}=(X^\top X)^{-1}+\frac{(X^\top X)^{-1}x_ix_i^\top(X^\top X)^{-1}}{1-h_{ii}}$, and substituting it into the normal equations gives the leave-one-out identity above without refitting. Since $H$ is an orthogonal projection, its trace equals its rank, $p$, so the average leverage is $p/n$. A point with $h_{ii}$ near 1 forces $\hat y_i\approx y_i$ whatever its residual would otherwise be, which is why its residual looks small.

Statement

In the linear model with homoskedastic errors, the variance of the $j$-th OLS coefficient is

where $S_{jj}=\sum_i(x_{ij}-\bar x_j)^2$ and $R_j^2$ is the $R^2$ from regressing column $j$ on the other columns. Equivalently, the eigenvalues of $X^\top X$ set the variance in each direction: $\operatorname{Var}(v^\top\hat\beta)=\sigma^2\sum_k (v^\top u_k)^2/\lambda_k$.

Why it's true

By Frisch–Waugh–Lovell, $\hat\beta_j$ is the slope of $y$ on the part of $x_j$ not explained by the other columns. That residualized column has sum of squares $S_{jj}(1-R_j^2)$, and the slope's variance is $\sigma^2$ divided by that. A thin direction of the cloud has a small eigenvalue $\lambda_k$ and a large variance. Because the data barely vary along it, nothing there constrains the fit.

Statement

Let $X$ be an exposure, $Y$ an outcome and $Z$ a third variable. A set $Z$ identifies the total effect of $X$ on $Y$ by adjustment if it satisfies the backdoor criterion: it blocks every path from $X$ to $Y$ that starts with an arrow into $X$, and contains no descendant of $X$. Three canonical roles:

Why it's true

Conditioning on $Z$ blocks a path through a non-collider and opens a path through a collider. If $X$ and $Y$ are independent causes of $Z$, then among units with the same $Z$, a high $X$ implies a low $Y$: this is Berkson's bias, and in a regression it appears as a nonzero coefficient on $X$ where there's no effect. Adjusting for a mediator removes the part of the effect that travels through it, leaving only the direct effect.

Statement

If $\operatorname{Var}(\varepsilon_i)=\sigma_i^2$, the exact covariance of the OLS estimator is $(X^\top X)^{-1}X^\top\Omega X(X^\top X)^{-1}$ with $\Omega=\operatorname{diag}(\sigma_i^2)$, not $\sigma^2(X^\top X)^{-1}$. White's estimator replaces $\sigma_i^2$ with squared residuals:

Why it's true

The middle matrix $\sum_i\sigma_i^2x_ix_i^\top$ is estimated by $\sum_i e_i^2x_ix_i^\top$. Each $e_i^2$ is a poor estimate of its own $\sigma_i^2$, but the sum averages the errors out, so $\frac1n\sum e_i^2x_ix_i^\top\to\frac1n\sum\sigma_i^2x_ix_i^\top$ and the sandwich is consistent whatever the form of $\sigma_i^2$. Residuals are shrunk by $\sqrt{1-h_{ii}}$ (idea 1), so $e_i^2$ underestimates $\sigma_i^2$ most where leverage is high. HC3 divides that out and is approximately the jackknife variance.

Statement

Create $m$ completed data sets by drawing the missing values from their posterior predictive distribution, estimate $\hat Q_k$ and its variance $U_k$ in each, and combine:

Intervals use a $t$ distribution with $\nu=(m-1)\big(1+\bar U/((1+1/m)B)\big)^2$ degrees of freedom. The fraction of missing information is $\lambda\approx(1+1/m)B/T$.

Why it's true

Total variance splits into within- and between-imputation parts: $\operatorname{Var}(\hat Q)=E[\operatorname{Var}(\hat Q\mid Y_{\rm obs},Y_{\rm mis})]+\operatorname{Var}(E[\hat Q\mid Y_{\rm obs},Y_{\rm mis}])$. The first term is $\bar U$, and the second is the spread of the estimate across plausible fills, $B$. A single imputation has $B=0$ by construction. Mean or regression imputation also distorts the relationship: the filled values sit exactly on the prediction line, so correlations are inflated or deflated and the variance of the variable is understated.

Statement

Suppose $K$ candidates have validation scores $\hat s_k=s_k+\sigma\xi_k$ with $\xi_k$ i.i.d. standard normal. The selected score $\max_k\hat s_k$ is biased upward relative to the true score $s_{k^*}$ of the model it picked. When all $s_k$ are equal,

Why it's true

The maximum of $K$ standard normals concentrates near $\sqrt{2\ln K}$, because the tail $\Pr[\xi>t]\approx e^{-t^2/2}$ hits $1/K$ there. The selected model's true score hasn't changed, only its observed score. If the candidates differ, selection still picks a mixture of real quality and luck, and the luck component is positive on average.

Statement

For a model with $k$ parameters and maximized log-likelihood $\hat\ell$ on $n$ observations,

For Gaussian regression, $-2\hat\ell=n\ln(\mathrm{RSS}/n)+\text{const}$. AIC is asymptotically equivalent to leave-one-out cross-validation (Stone, 1977) and is asymptotically efficient: it picks the candidate with the lowest prediction error. BIC is consistent: if the true model is among the candidates it is selected with probability tending to 1. AIC isn't consistent, because with probability bounded away from 0 it selects an over-large nested model.

Why it's true

The optimism of a maximized log-likelihood as an estimate of expected out-of-sample log-likelihood is $k$ in regular models, which gives the $2k$. BIC arises from a Laplace approximation to the log marginal likelihood, $\ln p(y\mid M)\approx\hat\ell-\frac k2\ln n$. Adding a spurious parameter raises $-2\hat\ell$ by a $\chi^2_1$ variable, and AIC's penalty of 2 is below the $\chi^2_1$ 95th percentile of 3.84, so it admits a useless parameter about 16% of the time, however large $n$ is. BIC's $\ln n$ eventually exceeds any fixed quantile.

Statement

A scoring rule $S(p,y)$ is proper if the expected score is optimized by reporting the true probability, and strictly proper if only by it. The Brier score $(p-y)^2$ and the log score $-\log p_y$ are strictly proper. A predictor is calibrated if $\Pr[Y=1\mid \hat p=p]=p$ for all $p$. Murphy's decomposition of the Brier score is

Why it's true

Condition on $\hat p$: $\mathbb E[(\hat p-Y)^2\mid\hat p]=(\hat p-\pi)^2+\pi(1-\pi)$, minimized at $\hat p=\pi$, which is propriety and the reliability term. Any strictly increasing transform of $\hat p$ leaves the ranking, and so AUC and the best-threshold accuracy, unchanged, while moving reliability. A predictor can therefore have an AUC of 0.95 and say "99%" when it's right 80% of the time.

Statement

Let training inputs follow $p(x)$ and deployment inputs $q(x)$, with $p(y\mid x)$ unchanged. For any loss $\ell$,

so the weighted training loss is an unbiased estimate of deployment loss. The price is variance: the effective sample size is $n_{\rm eff}=(\sum w_i)^2/\sum w_i^2$.

Why it's true

This is the change-of-measure identity, valid when $q\ll p$. If the model class contains the truth, both the weighted and unweighted fits are consistent for it and weighting only costs variance. If the model is misspecified, the unweighted fit minimizes error where the training data are dense, not where the deployment data are, and the two best-in-class fits differ (Shimodaira, 2000). Weighting trades that bias for variance.

Statement

For $y_{ij}=\theta_j+\varepsilon_{ij}$ with $\theta_j\sim N(\mu,\tau^2)$ and $\varepsilon_{ij}\sim N(0,\sigma^2)$, the posterior mean (BLUP) of group $j$ with $n_j$ observations is

The intraclass correlation is $\rho=\tau^2/(\tau^2+\sigma^2)$, and a cluster of size $m$ carries the information of $m/(1+(m-1)\rho)$ independent observations (the design effect is $1+(m-1)\rho$).

Why it's true

$\bar y_j$ is an unbiased but noisy estimate with variance $\sigma^2/n_j$, and the prior says $\theta_j$ is near $\mu$ with variance $\tau^2$. The precision-weighted combination has weights proportional to the reciprocal variances, which gives $B_j$. Groups with little data (large $\sigma^2/n_j$) are shrunk heavily and groups with lots of data are barely moved. For $J\ge4$ groups, this has lower total squared error than the separate means whatever the true $\theta_j$ (Stein; Efron & Morris, 1975).

Statement

For a stationary series with autocorrelations $\rho_k$, the variance of the sample mean is

Why it's true

Covariance between observations $k$ apart adds $2\rho_k\sigma^2/n$ for each lag to the variance. For AR(1), $\rho_k=\varphi^k$ and the geometric sum gives the factor above. With $\varphi=0.8$, $n_{\rm eff}=n/9$: a thousand daily observations are worth a hundred independent ones. Negative autocorrelation does the reverse, so alternating series are more informative about their mean.

Statement

A parameter is structurally identifiable if different values give different distributions of the data, and practically identifiable if the data at hand can distinguish them. Practical identifiability is read from the profile likelihood:

For a model with Jacobian $J$, the Fisher information $J^\top J/\sigma^2$ having eigenvalues spread across many orders of magnitude is the signature of a sloppy model (Brown & Sethna, 2003).

Why it's true

In $y=Ae^{-k_1t}+Be^{-k_2t}$ the two columns $\partial y/\partial k_1$ and $\partial y/\partial k_2$ become nearly parallel as $k_2\to k_1$, so $J^\top J$ is nearly singular and the curvature of the likelihood along the corresponding direction is tiny. The likelihood is flat along a ridge: raising $k_2$ while lowering $A$ and $B$'s balance leaves the curve almost unchanged. Moving along the ridge costs little in $\ell$, so the profile interval is wide.