BLUE: a guarantee with conditions

Best Linear Unbiased Estimator — OLS is BLUE when its assumptions hold

• Theorem about properties of the OLS sampling distribution
• Unbiased: on average, the estimate equals the true value — \(E[\hat{\beta}] = \beta\)
• Efficient (Best): among all linear unbiased estimators, OLS has the smallest variance
• These are different properties, protected by different assumptions

Unbiasedness and efficiency: visualised

Blue: correct answer, precise. Gold: correct direction, imprecise. Red: wrong answer entirely.

The assumption checklist

Why \(E[u \mid X] = 0\) fails: five routes

• Omitted variable bias — a relevant determinant of \(y\) is left in \(u\) and correlated with \(X\) (e.g. ability in wage regressions)
• Measurement error in regressors — if \(X\) is observed with error correlated with the true regressor, estimates are attenuated or biased
• Simultaneity / reverse causality — \(X\) affects \(y\) and \(y\) affects \(X\); the regressor is jointly determined with the outcome (price–quantity, income–health)
• Dynamic models with serially correlated errors — including \(y_{t-1}\) as a regressor while \(u_t\) is autocorrelated makes the lag endogenous
• Functional form misspecification — wrong transformation or omitted interactions push systematic structure into \(u\), which then correlates with included regressors

Our diagnostic dashboard

One function covers five of the six checks:

performance::check_model(model_base)

Linearity in parameters ≠ linearity in variables

OLS requires the model to be linear in the parameters — not in the raw variables.

Residuals vs. fitted: the primary diagnostic

The red smoother should be flat. A systematic curve means something is missing.

Hotels-europe: linear vs. log-log

The RESET test

Formal test for functional form misspecification (lmtest::resettest()).

resettest(m_linear)

    RESET test

data:  m_linear
RESET = 7.1487, df1 = 2, df2 = 9644, p-value = 0.00079

Fix 1 — Log transformation for a multiplicative relationship

lm(log(price) ~ log(distance + 0.1), data = hotels)

Fix 2 — Quadratic term for a curved relationship

lm(price ~ distance + I(distance^2), data = hotels)

Non-constant error variance

OLS assumes \(\text{Var}(\varepsilon_i) = \sigma^2\) for all \(i\) — the same uncertainty everywhere.

What goes wrong when this fails:

• OLS estimates are still unbiased — the direction and magnitude are correct on average
• But standard errors are inconsistent — too narrow in some places, too wide in others
• t-tests and confidence intervals mislead you about precision

Scale-location plot: the visual test

The red smoother should be flat and horizontal.

Hotels-europe: heteroskedasticity in practice

The Breusch-Pagan test

bptest(model_base)

    studentized Breusch-Pagan test

data:  model_base
BP = 30.437, df = 2, p-value = 2.459e-07

The fix: Heteroskedasticity-robust (HC) standard errors

• Replace standard formula with a “sandwich” estimator
• Uses each observation’s actual squared residual to estimate uncertainty.
• Result: larger SEs where the model struggles, smaller where it fits well.
  • Estimates do not change.

modelsummary(list("Standard" = model_base, "Robust (HC1)" = model_base),
             vcov = list("classical", "HC1"), stars = TRUE)

When predictors move together

Perfectly correlated predictors: OLS breaks down

VIF: the Variance Inflation Factor

VIF: by how much is variance of \(\hat{\beta}_j\) inflated relative to a model with no correlation among predictors.

\[\text{VIF}_j = \frac{1}{1 - R^2_j}\]

\(R^2_j\): the \(R^2\) from regressing predictor \(j\) on all other predictors.

Hotels-europe: low vs. problematic VIF

Base model — distance and stars measure different things: VIF is fine.

vif(model_base)   # price ~ distance + stars

distance    stars 
1.003385 1.003385 

model_collinear <- lm(price ~ distance + stars + log(stars), data = hotels)
vif(model_collinear)

  distance      stars log(stars) 
  1.004235  20.543854  20.526349 

What to do about multicollinearity

• Do not drop variables mechanically based on VIF — think about what each variable represents
• Dropping an important control to reduce VIF can introduce omitted variable bias (a far worse problem)
• If two variables genuinely measure the same construct, keep only one — or combine them into a single index (e.g. a simple average, or the first principal component)
• Sometimes collinearity is just a feature of reality — acknowledge it and note that individual coefficients should be interpreted with caution

Three types of unusual data points

Cook’s distance

\[D_i = \frac{(\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(-i)})^\top (X^\top X) (\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}_{(-i)})}{p \cdot \hat{\sigma}^2}\]

Example: one influential point

The dashed line is the original fit. One point pulls the slope down substantially.

Hotels-europe: who are the influential hotels?

augment(model_base) |>
  bind_cols(hotels |> select(hotel_id, city)) |>
  arrange(desc(.cooksd)) |>
  select(hotel_id, city, stars, price, .cooksd) |>
  head(8)

# A tibble: 8 × 5
  hotel_id city      stars price .cooksd
     <dbl> <chr>     <dbl> <dbl>   <dbl>
1    14797 Prague        4  3714  0.0477
2    14797 Prague        4  3714  0.0477
3     1392 Barcelona     1  1213  0.0313
4     1392 Barcelona     1  1213  0.0313
5     8321 London        5  1848  0.0295
6    13394 Paris         5  1364  0.0139
7    14307 Paris         5  1265  0.0112
8    12644 Paris         5  1218  0.0101

Investigate — don’t delete

• An influential observation is not automatically an error. It may be the most interesting case in the dataset.
• Ask: why is this point influential? Data entry error, genuine extreme case, or a different sub-population?
• Workflow: re-estimate without the flagged observations. Report both results and note the sensitivity.
• Silently removing influential points without disclosure is a serious methodological error.

The Q-Q plot: good vs. violated

A widespread misconception

What normality actually protects: Exact finite-sample inference — t- and F-statistics follow their named distributions exactly only when errors are normal.

The Shapiro-Wilk trap

shapiro.test(resid(model_base)[1:5000])  # max n for shapiro.test is 5000

    Shapiro-Wilk normality test

data:  resid(model_base)[1:5000]
W = 0.80384, p-value < 2.2e-16

What the PPC shows

The PPC simulates new response values from the fitted model many times and compares their distribution to the observed data.

PPC: level-level vs. log-log

• Blue line: density of the observed outcome \(y\)
• Green lines: densities of outcomes simulated from the fitted model
• If green ≈ blue → the model generates plausible data → distributional assumptions are reasonable
• If they diverge → the model produces values in the wrong range → mis-specified DGP

Dashboard panels 1–2: functional form and variance

• Linearity — residuals vs. fitted; the red smoother should be flat. Systematic curve → functional form problem
• Homogeneity of variance — scale-location; the smoother should be horizontal. Remaining fan shape → use robust SEs even with the log model

Dashboard panels 3–4: collinearity and normality

• Collinearity (VIF) — bars should be short (< 5). Both predictors are fine; recall the stars + log(stars) example where VIF ≈ 20
• Normality of residuals (Q-Q) — much closer to the diagonal than model_base (skewness −0.1 vs. 7.8); with \(n > 9{,}000\) the remaining departure does not affect inference

Dashboard panels 5–6: influential observations and overall fit

• Outliers / influential observations — Cook’s D bar and leverage–residuals scatter; points outside the contour lines warrant investigation (report with and without, do not silently delete)
• Posterior predictive check — green lines now closely track the blue observed density. The log transformation resolved the distributional mismatch we saw with model_base.

What residuals cannot tell you

All diagnostic panels show properties of the residuals — the part of \(y\) the model did not explain.

But every route to endogeneity (\(E[u \mid X] \neq 0\)) is invisible in the residuals:

• Omitted variable bias — a confounder is absorbed into \(\hat{\beta}\); the fit looks fine
• Measurement error in regressors — attenuation bias hides behind well-behaved residuals
• Simultaneity / reverse causality — the endogenous regressor leaves no obvious residual pattern

The mechanism

Suppose the true model is:

\[y = \beta_0 + \beta_1 x + \beta_2 z + \varepsilon\]

but we estimate \(y = \beta_0 + \beta_1 x + \varepsilon\) (omitting \(z\)).

OVB in practice: what changes when we add city?

m_no_city   <- lm(log(price) ~ log(distance + 0.1), data = hotels)
m_with_city <- lm(log(price) ~ log(distance + 0.1) + city, data = hotels)

Why the residuals looked fine

Run the full diagnostic suite on the biased model (no city controls):

The way forward

Panel data and fixed effects (panel data session):

• If the omitted variable is time-invariant (e.g. which city), fixed effects eliminate it
• One of the most powerful practical solutions to OVB

The diagnostic workflow

For every regression, ask:

1. Residuals vs. fitted — is the functional form correct? Is there a pattern?
2. Scale-location — is the variance constant? Do I need robust SEs?
3. Cook’s distance / leverage — are a few observations driving my results?
4. VIF — are my predictors too correlated to separate their effects?
5. Q-Q plot — are there serious departures from normality? (Large \(n\): rarely a concern)
6. Thinking — what important variable might I have left out?

performance::check_model(your_model)  # panels 1–5 in one call