Session 3: Adequacy, Validity, and Robustness of Regression Models

Published course presentation · repository edition

Thierry Warin, PhD

HEC Montréal

Outline · The hour ahead

  1. From data to models
  2. Why evaluate a model at all?
  3. Adequacy · how much does the line explain?
  4. Validity · do the assumptions hold?
  5. Robustness · which observations move the answer?
  6. Parsimony · comparing candidate models
  7. Best practice, and what to take away

What session 3 is for

The mathematics

  • Derive the least-squares slope, and read “controlling for” as a residual regression
  • Read residuals versus fitted and say what structure it reveals
  • Compute leverage from the hat matrix, and show it sums to \(p\)
  • Separate an outlier from an influential point with Cook’s distance
  • Compare models on AIC, and say why that is not model selection

The Python

  • Fit a simple and a multiple regression with smf.ols, and read .params
  • Pull diagnostics from a fitted model with get_influence()
  • Build the four diagnostic plots yourself rather than calling one function
  • Compare candidate specifications side by side on AIC and BIC
  • Report the comparison as a publication table with qmib.regtable

The international business

  • Quote a slope in units a manager can act on
  • Say why a slope is an association, and what more is needed to call it a cause
  • Say whether a published coefficient survives its own diagnostics
  • Decide what to do about the one country that moves the whole result
  • Refuse a model whose residuals show it has missed the structure

Your Python so far

Everything you need to type is on one sheet — how to start Python in VS Codium on your platform, how to make a file and run it, the DataFrame basics, and every name from sessions 1 to 3, each with a line showing it in use.

Python Cheatsheet · Session 03 (PDF) · in the repository at 03-…/PYTHON-CHEATSHEET.pdf

One sheet per package, functions explained with pictures: pandas · NumPy · matplotlib · statsmodels

Print it, or keep it open beside the editor. Nothing on it is examinable; all of it is assumed.

From data to models

The question

Do more productive European countries have higher income per head?

Thirty countries, fifteen years, both numbers already in the repository.

No download, no API key: data/spine/core.parquet.

Always look first

import pandas as pd
import matplotlib.pyplot as plt

core = pd.read_parquet("data/spine/core.parquet")

# Drop rows missing either variable. subset= matters: without it, a row missing
# any column at all would go, and you would silently lose a third of the panel.
d = core.dropna(subset=["gdp_pc_eur", "productivity_idx"])

# Look before you model. s= is the marker size and alpha= the transparency —
# with 428 points, overlapping markers hide density unless you can see through
# them.
plt.scatter(d["productivity_idx"], d["gdp_pc_eur"], s=12, alpha=0.6)
plt.xlabel("Labour productivity index (EU average = 100)")
plt.ylabel("GDP per capita (EUR)")
plt.show()

Draw the picture before you fit anything. A straight line through a curved cloud is still a straight line, and the software will not warn you.

The model, and what each piece is

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

Symbol Name Here
\(y_i\) outcome GDP per capita
\(x_i\) predictor productivity index
\(\beta_0\) intercept predicted \(y\) when \(x = 0\)
\(\beta_1\) slope change in \(y\) per unit of \(x\)
\(\varepsilon_i\) error what \(x_i\) does not explain

Choosing the line

\[\min_{\beta_0,\,\beta_1} \sum_{i=1}^{n}\big(y_i - \beta_0 - \beta_1 x_i\big)^2\]

Squaring makes every miss positive, so misses above and below cannot cancel — and it punishes one large miss more than several small ones.

This is the mean, done conditionally. The mean minimises \(\sum(x_i - c)^2\); the regression line minimises the same quantity over lines.

Three lines of Python

import statsmodels.formula.api as smf

# The formula reads "gdp_pc_eur explained by productivity_idx". statsmodels
# adds the intercept for you; write ~ 0 + x if you ever want it suppressed,
# and know that R-squared stops being interpretable when you do.
model = smf.ols("gdp_pc_eur ~ productivity_idx", data=d).fit()

# .params holds the two coefficients. The slope is in euros of GDP per capita
# per one index point — always say the units.
print(model.params)
Intercept          -86118.1
productivity_idx     1111.3

Read the ~ as “explained by”. .fit() does the arithmetic; model.params holds the estimates.

The slope is the sentence you will write

Comparing two country-years that differ by one point of productivity, the more productive one has on average €1,111 more GDP per capita.

Always say the units aloud — euros per index point. A slope without units is not a finding.

The intercept says a country at productivity 0 would have −€86,118. The index runs from 72 to 129 in our data. Do not interpret an intercept you have no data near.

Fitted values and residuals

Germany, 2022 — productivity 121.2, actual GDP per capita €54,239:

\[\hat y = -86{,}118 + 1{,}111.3 \times 121.2 = 48{,}606\]

\[e = 54{,}239 - 48{,}606 = \mathbf{+5{,}633}\]

Germany sits €5,633 above the line. That residual is everything about German income that productivity alone does not account for — and reading residuals is the rest of this hour.

How much did the line explain?

# R-squared is the share of variance the line accounts for — and nothing about
# whether the line is the right shape. The rest of this hour is that difference.
print(f"R-squared = {model.rsquared:.3f}")
R-squared = 0.627

\(R^2\) is the share of the variation in \(y\) the model accounts for: 62.7%, on a scale from 0 to 1.

It does not say the model is correct, that the relationship is causal, or that the line will predict a new country well.

Two predictors: what “controlling for” means

# The third column has its own holes, so the sample is rebuilt: d2 is not d.
d2 = core.dropna(subset=["gdp_pc_eur", "productivity_idx", "gfcf_meur"])

# Adding a second predictor changes what the first coefficient means: it is now
# the association holding investment fixed, not the raw association. The two
# answer different questions, and the number will move.
model2 = smf.ols("gdp_pc_eur ~ productivity_idx + gfcf_meur", data=d2).fit()
Intercept          -71211.8
productivity_idx      928.3
gfcf_meur              68.7
R-squared = 0.692

The productivity slope falls 1,111 → 928, about 16%. It now answers a different question: comparing country-years with the same investment, one point more productivity goes with €928 more income.

Correlation is not causation

The slope says the two move together. It does not say raising productivity by a point would raise income by €928.

  1. Productivity raises income
  2. Higher income buys the capital and training that raise productivity — the arrow points the other way
  3. Something else — institutions, education, market size — drives both

Nothing in the output distinguishes these. Until Sessions 10–11, write “is associated with” — and mean it.

Why evaluate a model at all?

What a regression table asks you to believe

Every number in it is an estimate, and every star is a claim about a population you never saw. Today is about which of those claims survive.

“All models are wrong, but some are useful”

— George Box

  • Regression results must be tested for adequacy, not just reported.
  • Ignoring diagnostics can lead to:
    • False significance
    • Overconfidence in predictions
    • Poor generalization

The four questions of the hour

  • Model Accuracy (Adequacy): \(R^2\), Adjusted \(R^2\), RSE.
  • Model Validity (Assumptions): Linearity, homoscedasticity, normality of residuals.
  • Model Robustness: Outliers, leverage, Cook’s distance (sensitivity to unusual points).
  • Model Parsimony/Selection: AIC, BIC — and, in session 05, what to do instead of stepwise.

Adequacy · how much does the line explain?

Adequacy: two things to assess

We can assess the adequacy of

  • the coefficient estimates [individually]

  • the model [the whole]

First, assessing the adequacy of the coefficient estimates:

We continue the analogy with the estimation of the population mean \(\mu\) of a random variable \(Y\).

A natural question is as follows:

  • how accurate is the sample mean \(\hat{\mu}\) as an estimate of \(\mu\)?

We have established that the average of \(\hat{\mu}\)’s over many data sets will be very close to \(\mu\), but that a single estimate \(\hat{\mu}\) may be a substantial underestimate or overestimate of \(\mu\).

Residual Standard Error

Residuals are the leftovers from the model fit: Data = Fit + Residual

A residual is the difference between the observed \(y_i\) and predicted \(\hat{y}_i\).

\[\epsilon_i = y_i - \hat{y}_i\]

  • living in poverty in DC is 5.44% more than predicted.

  • living in poverty in RI is 4.16% less than predicted.

How far off will that single estimate of \(\hat{\mu}\) be?

  • In general, we answer this question by computing the standard error of \(\hat{\mu}\), written as \(SE(\hat{\mu})\).

We have the well-known formula:

\[Var(\hat{\mu}) = SE(\hat{\mu})^2 = \frac{\sigma^2}{n}\] where \(\sigma\) is the standard deviation of each of the realizations \(y_i\) of \(Y\).

Roughly speaking, the standard error tells us the average amount that this estimate \(\hat{\mu}\) differs from the actual value of \(\mu\).

The previous equation also tells us how this deviation shrinks with \(n\) — the more observations we have, the smaller the standard error of \(\hat{\mu}\).

Both appear in the formula above, and they are not the same quantity.

measures the spread of with more data
standard deviation \(\sigma\) the observations, around their mean does not shrink — you estimate a fixed \(\sigma\) more precisely
standard error \(\mathrm{SE}(\hat\mu) = \sigma/\sqrt{n}\) an estimate, around the truth, across repeated samples shrinks as \(1/\sqrt{n}\)

They carry the same units, which is why they get swapped. But \(\sigma\) is a fact about the population your sample came from, and \(\mathrm{SE}\) is a fact about your estimate of it: the standard error is the standard deviation of \(\hat\mu\) over the samples you might have drawn but did not.

So: an error bar that stays the same width as \(n\) grows is a standard deviation; one that narrows is a standard error. Collecting more data makes you more certain where the mean is; it does not make the observations less spread out.

In a similar vein, we can wonder how close \(\hat{\beta}_0\) and \(\hat{\beta}_1\) are to the true values \(\beta_0\) and \(\beta_1\)?

  • To compute the standard errors associated with \(\hat{\beta}_1\), it comes down to the following formula:

\[ SE(\hat{\beta_1}) \]

Models: Assessing the adequacy of the coefficient estimates

Standard errors can be used to compute our confidence intervals, i.e. that the true value of \(\beta_1\), will within this interval with a 95% chance:

\[ \hat{\beta}_1 - 1.96 \times SE(\hat{\beta}_1),~\hat{\beta}_1 + 1.96 \times SE(\hat{\beta}_1) \]

Standard errors can also be used for hypothesis tests:

\[ H_0:~\beta_1=0 \]

\[ H_1:~\beta_1 \neq 0 \]

As already seen, we compute the the t-statistic, which measures the number of standard deviations that \(\hat{\beta}_1\) is away from 0 (see p-values from last session):

\[ t=\frac{\hat{\beta}_1-0}{SE(\hat{\beta}_1)} \]

From residuals to \(R^2\)

Assessing the adequacy of the model now:

Once we have rejected the null hypothesis in favor of the alternative hypothesis, we want to quantify the extent to which the model fits the data.

For that, we use:

  • the residual standard error, and

  • the \(R^2\) statistic.

The RSE is an estimate of the standard deviation of \(\epsilon\).

Definition: It is the average amount that the response will deviate from the true regression line.

\[ RSE = \sqrt{\frac{RSS}{n-2}} = \sqrt{\frac{\sum_{i=1}^{n}(y_i-\hat{y}_i)^2}{n-2}} \]

recall:

\[ RSS= \sum_{i=1}^{n}(y_i-\hat{y}_i)^2 \]

Let’s define TSS as the total sum of squares.

\[ TSS = \sum_{i=1}^{n}(y_i - \bar{y})^2 \]

  • TSS measures the total variance in the response Y.

Definition: The \(R^2\) statistic is the proportion of variance explained.

To calculate \(R^2\), we use the formula:

\[R^2 = \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS}\]

  • It tells us what percent of variability in the dependent variable is explained by the model.

  • The remainder of the variability is explained by variables not included in the model or by inherent randomness in the data.

Adjusted \(R^2\) and RSE, in Python

\(R^2\) can never fall when a predictor is added, so on its own it cannot say whether that predictor earned its place. Two answers:

  • Adjusted \(R^2\): corrects for number of predictors.
  • RSE (Residual Standard Error): average deviation from regression line.
import qmib
import statsmodels.formula.api as smf

core = qmib.load("core")
d = core.dropna(subset=["gdp_pc_eur", "productivity_idx"])

# the model we fitted at the top of the hour — now: can it be trusted?
fit = smf.ols("gdp_pc_eur ~ productivity_idx", data=d).fit()

# Adjusted R^2 penalises each extra predictor, so it can fall when you add one.
print(f"R2 {fit.rsquared:.4f}   adjusted {fit.rsquared_adj:.4f}")

# RSE is the typical miss, in euros — the same units as the outcome, which is
# what makes it the more useful of the two to quote to a non-specialist.
import numpy as np
print(f"RSE {np.sqrt(fit.mse_resid):,.0f} EUR per head")
R2 0.6273   adjusted 0.6264
RSE 9,376 EUR per head

Interpretation:

  • Adjusted \(R^2\) can decrease if irrelevant variables are added.
  • RSE gives an absolute measure, in units of Y.

Quiz: reading an \(R^2\)

Which of the below is the correct interpretation of \(R^2 = 0.38\)?

  • 38% of the variability in the % of HG graduates among the 51 states is explained by the model.

  • 38% of the variability in the % of residents living in poverty among the 51 states is explained by the model.

  • 38% of the time % HS graduates predict % living in poverty correctly.

  • 62% of the variability in the % of residents living in poverty among the 51 states is explained by the model.

Which of the below is the correct interpretation of \(R^2 = 0.38\)?

  • 38% of the variability in the % of HG graduates among the 51 states is explained by the model.

  • 38% of the variability in the % of residents living in poverty among the 51 states is explained by the model.

  • 38% of the time % HS graduates predict % living in poverty correctly.

  • 62% of the variability in the % of residents living in poverty among the 51 states is explained by the model.

Validity · do the assumptions hold?

The residual plot is the instrument

Residual = observed \(y_i\) – predicted \(\hat{y}_i\).

  • Residuals contain all the information the model did not explain.
  • Checking them reveals whether assumptions of linear regression hold.
# Residuals vs Fitted: is there structure the straight line missed?
plt.scatter(fit.fittedvalues, fit.resid, s=12, alpha=0.6)
plt.axhline(0, color="black", lw=1)
plt.xlabel("Fitted values"); plt.ylabel("Residuals")
plt.show()

What you are looking for

Scattered at random around zero is the model behaving. A curve, a funnel or a drift is the model telling you what it missed.

Three conditions to check

    1. Linearity
    1. Nearly normal residuals
    1. Constant variability

Validity: linearity

  • The relationship between the explanatory and the dependent variable should be linear.

Check using a scatterplot of the data, or a residuals plot.

Anatomy of a residuals plot

\[\%~HSgrad = 81\%~\&~\%~in~poverty = 10.3\]

\[\widehat{\%~in~poverty} = 64.68 - 0.62 \times 81 = 14.46\]

\[\epsilon = {\%~in~poverty} - \widehat{\%~in~poverty}\]

\[\epsilon = 10.3 - 14.46 = -4.16\]

\[\%~HSgrad = 86\%~\&~\%~in~poverty = 16.8\]

\[\widehat{\%~in~poverty} = 64.68 - 0.62 \times 86 = 11.36\]

\[\epsilon = {\%~in~poverty} - \widehat{\%~in~poverty}\]

\[\epsilon = 16.8 - 11.36 = 5.44\]

  • The relationship between predictors and outcome should be linear.
  • Violations show up as curved patterns in residuals.

Remedies: add polynomial terms, interactions, or reconsider functional form.

Validity: nearly normal residuals

The residuals should be nearly normal — and nearly is the operative word. Least squares does not need normal errors to give you unbiased coefficients; it needs them for the t-tests and confidence intervals to mean what they claim. In a large sample even that requirement softens, because the central limit theorem does the work for you.

So the question is not “are these residuals normal” — they never are — but how badly does the departure distort the inference. Heavy tails on both sides usually signal a handful of unusual observations rather than a broken model; a single long tail on one side usually means the dependent variable wanted a log.

Check using a histogram.

Validity: constant variance

  • Residuals should have constant variance across fitted values.
  • Violation = heteroscedasticity (funnel shape).
  • Remedies: transform Y, add missing predictors, or use robust errors.
import matplotlib.pyplot as plt

# Scale-Location: is the spread of the residuals constant?
plt.scatter(fit.fittedvalues, np.sqrt(np.abs(fit.get_influence().resid_studentized_internal)),
            s=12, alpha=0.6)
plt.xlabel("Fitted values"); plt.ylabel("sqrt(|standardised residual|)")
plt.show()

  • The variability of points around the least squares line should be roughly constant.

  • This implies that the variability of residuals around the 0 line should be roughly constant as well.

  • Also called homoskedasticity.

Check using a residuals plot.

Quiz: reading a residual plot

What condition is this linear model obviously violating?

  • Constant variability

  • Relationship

  • Normal residuals

  • No extreme outliers

What condition is this linear model obviously violating?

  • Constant variability

  • Relationship

  • Normal residuals

  • No extreme outliers

What condition is this linear model obviously violating?

  • Constant variability

  • Linear relationship

  • Normal residuals

  • No extreme outliers

What condition is this linear model obviously violating?

  • Constant variability = heteroskedastic = non-homoskedastic

  • Linear relationship

  • Normal residuals

  • No extreme outliers

Practise the eye

Guess the correlation, then check. Ten minutes here is worth a chapter.

Correlation and Regression Game

Robustness · which observations move the answer?

Outlier, leverage point, influential point

Three words that get used interchangeably in conversation and mean quite different things here. Keeping them apart is most of the work.

Outlier unusual Y given X the model predicts it badly
Leverage point unusual X it sits far from the other predictor values
Influential point removing it moves the coefficients the answer depends on it

A point can be any one of these without being the others. An outlier sitting in the middle of the X range is badly predicted but has little power to tilt the line: there are neighbours on both sides holding it in place. A leverage point that happens to fall exactly where the line already goes is unusual in X and perfectly consistent with everything else — it does no damage at all, and in fact it stabilises the slope.

The damage is done by points that are both. Far out in X, so the line pivots around them, and badly predicted, so the pivot is large. That combination is what influence means, and it is why neither residual size nor leverage on its own is the diagnostic you want.

Which is the argument for a measure that multiplies the two together. That measure is Cook’s distance, next.

Cook’s Distance

One question, asked of every row in turn: delete this observation, refit, and how far do all the fitted values move? Cook’s distance is the answer, and it is a single number per row.

\[D_i \;=\; \underbrace{\frac{e_i^2}{p\,s^2}}_{\text{is it badly predicted?}} \;\times\; \underbrace{\frac{h_i}{(1-h_i)^2}}_{\text{does it have leverage?}}\]

Those are the two ingredients from the previous slide, multiplied rather than added — which is the point. A big residual in the middle of the X range gives a large first factor and a tiny second one, so the product stays small. It takes both to move the line, and only the product knows that.

influence = fit.get_influence()
cooks_d = influence.cooks_distance[0]

print(f"largest Cook's D {cooks_d.max():.4f}")
print(f"above 4/n ({4/len(d):.4f}): {(cooks_d > 4/len(d)).sum()} observations")
largest Cook's D 0.0425
above 4/n (0.0093): 25 observations

\(4/n\) is a screen, not a test. It has no distributional justification — it simply flags the most extreme few per cent of a well-behaved sample. Crossing it is a reason to go and look at the row. It is never, on its own, a reason to delete it.

One spike per observation

One spike per row, its height the distance the fitted values would travel if that row were deleted. Read it for shape, not for threshold-crossing: a forest of similar spikes means no single country-year is carrying the result, while one tower above the rest means your finding is really a statement about one observation.

# The dashed line is the conventional 4/n screen — a prompt to go and look at
# the row, not a rule for dropping it.
plt.stem(cooks_d)
plt.axhline(4 / len(d), color="red", ls="--")
plt.xlabel("Observation"); plt.ylabel("Cook's distance")
plt.show()

When a spike stands out, the next move is always the same: find out which row it is and go and look at it.

d.iloc[cooks_d.argmax()]      # the observation, not just its index

Then ask what it is. A data-entry error you fix. A small economy that genuinely behaves differently is not an error — it is a finding, and dropping it to tidy the plot is how a result gets manufactured. If you do exclude a row, the exclusion and the reason for it belong in the write-up, alongside the fit with the row left in.

Parsimony · comparing candidate models

Information criteria — AIC and BIC

Adding a variable can never make the fit worse, so “which model fits best” always answers “the biggest one”. Information criteria fix that by charging rent for each parameter:

\[\text{AIC} = -2\log(L) + 2k \qquad\qquad \text{BIC} = -2\log(L) + k\ln(n)\]

The first term is the same in both — it is the fit, and it falls as the model grows. The second term is the rent, and it is the only thing that differs:

Rent per parameter Which means
AIC \(2\) a flat charge, whatever the sample size
BIC \(\ln(n)\) the charge rises with \(n\) — at \(n = 410\) it is \(6.02\), three times AIC’s

So BIC always prefers a model at least as small as AIC’s, and the gap widens with the sample. Neither is “correct”: AIC is aiming at the model that predicts best, BIC at the model that is true. Those are different targets, and they are allowed to disagree.

Watching them disagree

Four nested models on the same 410 rows — because two fits estimated on different samples are not a comparison.

import numpy as np

cols = ["gdp_pc_eur", "productivity_idx", "gfcf_meur", "employment_ths", "falling_behind"]
dic  = core.dropna(subset=cols)

specs = {
    "productivity":     "productivity_idx",
    "+ investment":     "productivity_idx + gfcf_meur",
    "+ employment":     "productivity_idx + gfcf_meur + employment_ths",
    "+ falling_behind": "productivity_idx + gfcf_meur + employment_ths + falling_behind",
}
fits = {k: smf.ols(f"gdp_pc_eur ~ {v}", data=dic).fit() for k, v in specs.items()}

ic = pd.DataFrame({"k":   [int(m.df_model) + 1 for m in fits.values()],
                   "AIC": [m.aic for m in fits.values()],
                   "BIC": [m.bic for m in fits.values()]}, index=list(fits))
ic["dAIC"] = ic.AIC - ic.AIC.min()      # distance from the best, which is what you read
ic["dBIC"] = ic.BIC - ic.BIC.min()

print(f"n = {len(dic)},  AIC's rent = 2,  BIC's rent = ln(n) = {np.log(len(dic)):.2f}\n")
print(ic.round(2).to_string())
print(f"\nAIC prefers: {ic.AIC.idxmin()}\nBIC prefers: {ic.BIC.idxmin()}")
n = 410,  AIC's rent = 2,  BIC's rent = ln(n) = 6.02

                  k      AIC      BIC    dAIC    dBIC
productivity      2  8666.05  8674.08  223.18  214.41
+ investment      3  8588.88  8600.93  146.02  141.26
+ employment      4  8443.60  8459.67    0.74    0.00
+ falling_behind  5  8442.87  8462.95    0.00    3.28

AIC prefers: + falling_behind
BIC prefers: + employment

Reading that table

The first two steps are unanimous: investment and employment are worth their rent several times over, and both criteria drop by more than a hundred points.

The last step is where they part. Adding falling_behind improves AIC by 0.74 and worsens BIC by 3.28 — and the difference between those two movements is not a coincidence:

\[3.28 - (-0.74) \;=\; 4.02 \;=\; \ln(410) - 2\]

That is the extra rent BIC charges for one more parameter. The entire disagreement is that one number. The data has not changed; the price of complexity has.

Read differences, not levels. An AIC of 8442 means nothing on its own — it depends on the units of \(y\) and the size of the sample. A \(\Delta\)AIC of 0.74 between two models means they are indistinguishable; a \(\Delta\) of 145 means one is out of the running. A rough convention: under 2 is a tie, over 10 is decisive.

And the standing condition: same rows, same \(y\). Two AICs computed on different samples, or one on \(y\) and one on \(\log y\), are not comparable at all — and dropna() silently changes the sample every time you change the column list.

The table you would actually publish

qmib.regtable is R’s stargazer, here: one column per model, stars on the coefficients, standard errors beneath, fit statistics at the foot.

import qmib

d2 = core.dropna(subset=["gdp_pc_eur", "productivity_idx", "gfcf_meur"])

# Both models on the SAME rows. dropna gives a different sample for each set of
# columns, and two fits estimated on different samples are not a comparison —
# their AIC, their R-squared and their coefficients are not on one scale.
m1 = smf.ols("gdp_pc_eur ~ productivity_idx", data=d2).fit()
m2 = smf.ols("gdp_pc_eur ~ productivity_idx + gfcf_meur", data=d2).fit()

print(qmib.regtable([m1, m2], names=["simple", "+ investment"],
                    order=["productivity_idx", "gfcf_meur", "Intercept"]))

===============================================
                    (1) simple (2) + investment
-----------------------------------------------
productivity_idx   1119.8***   928.3***        
                   (42.9)      (43.7)          
gfcf_meur                      68.7***         
                               (7.2)           
Intercept          -86949.7*** -71211.8***     
                   (4526.6)    (4421.4)        
Adjusted R-squared 0.623       0.690           
Observations       414         414             
R-squared          0.624       0.692           
Residual SE        9,434       8,548           
===============================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01

.as_latex() for a paper, .as_html() for a slide. Pass your own stats= dict to put different rows at the foot.

Choosing a model is not a diagnostic

AIC and BIC compare models you have already written down. They do not write them down for you.

The automated alternative — stepwise selection, which adds and drops predictors by criterion until nothing improves — is a selection procedure, not a check on a fit. Teaching it beside Cook’s distance would imply it belongs to model checking.

It also cannot be judged on a case like this one. With five candidates and 367 rows forward selection simply keeps all five, and nothing is learned. The procedure only looks decisive when there are many candidates and few rows — which is the case it is actually used in, and the case where it is least trustworthy.

So it opens session 05 instead, on 29 candidates and 37 rows, where it does real work: watch its choices fail to survive a resample and you have the argument for shrinking coefficients continuously rather than including or excluding them. That argument is ridge and lasso.

Best practice, and what to take away

Best Practices

  • Check residuals first — assumptions matter.
  • Use multiple fit criteria, not just \(R^2\).
  • Always investigate influential points.
  • Use AIC/BIC to compare, but remember they are relative.
  • If you selected on a criterion, say so, and treat the inference as exploratory.

Final message: Balance predictive adequacy, interpretability, and parsimony.

Conclusion

  • We fitted the smallest possible model, and read its slope in units
  • Adequacy asked how much it explains; validity, whether its assumptions hold; robustness, which observations are doing the work; parsimony, whether a simpler model would do as well
  • Not one of the four is answered by the coefficient table alone
  • Session 04 keeps the machinery and changes the outcome: the thing you are predicting becomes a category

One more thing…

After Wainer’s “most dangerous equation” — the standard error of the mean, and what forgetting it costs.

References

  • James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An Introduction to Statistical Learning with Applications in R (2nd ed.). Springer.
  • Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005). Applied Linear Statistical Models (5th ed.). McGraw-Hill Irwin.
  • Wooldridge, J. M. (2019). Introductory Econometrics: A Modern Approach (7th ed.). Cengage.

Before the next session

Required: the next session’s article, its Dataverse package, and the self-check. That is the pre-session deck.

Optional: the QMIB Lab for this session — a knowledge check on the hour you have just sat through: warin.ca/qmib-labs

Never a prerequisite, never marked for correctness. A completed report counts toward participation because attempting it is the engagement being measured.