Prof D’s Regression Sessions - Vol 4

Good Skool / Bad Skool Breaks - GLMs and logistic regression

Adam Dennett - a.dennett@ucl.ac.uk

Week GLM

CASA0007 Quantitative Methods

Recap

  • Over three weeks we built a regression model up from a single scatter plot to a linear mixed effects model with random intercepts and random slopes
  • We learned a recipe - understand your data, start simple, add complexity only when it earns its place, and check your assumptions at every step
  • Last week we relaxed one particular assumption: that our observations were independent. Schools sit inside local authorities, and that grouping matters

But one thing never changed

  • Every single model we have fitted has had a continuous outcome variable
  • Attainment 8. Progress 8. Numbers that can take any value, that are roughly bell-shaped, and that happily sit on either side of a straight line
  • This choice let us use Ordinary Least Squares throughout - but we know that our data don’t always look like this
  • What happens when the thing you want to explain is a yes/no, a category, or a count?

This week’s session aims

  • Understand why OLS fails when the outcome is not continuous, and see it fail
  • Meet the Generalised Linear Model (GLM) - one framework, three moving parts
  • Logistic regression for binary outcomes - in depth, because if you understand this one the rest follow
  • Where counts and ordered categories fit in, and where to go for them
  • How to judge one when there is no R² - deviance, AIC, the pseudo-R²s, and whether the predictions are actually any use
  • Throughout: can we predict a school’s Ofsted rating from everything else we know about it?

Back to Week 1

Remember these?

In week 1 we spent a long time on types of data, and on the distributions that go with them.

Data type Example from our data Distribution
Real-valued Attainment 8 Normal
Binary Is this school single-sex? Bernoulli / Binomial
Count How many pupils on free school meals? Poisson
Count, overdispersed … when the variance exceeds the mean Negative binomial
Ordinal Ofsted rating (no single answer)

We then spent three weeks modelling only the first row!

The whole idea of a GLM, in one table

If your outcome is… …its distribution is …the model is
A continuous number Normal OLS - weeks 5, 6, 7
A yes/no Binomial Logistic regression
A proportion out of a known total Binomial Logistic regression
An ordered category - Ordinal logistic
An unordered category - Multinomial logistic
A count Poisson Poisson regression
A count that is too variable Negative binomial Negative binomial regression

This is not seven new methods. It is one method, with a slot you fill differently depending on what you are trying to explain.

What actually changes

Every GLM has exactly three pieces. You already know the middle one.

  1. The random component - what distribution does the outcome follow? (Normal, binomial, Poisson…)
  2. The linear predictor - \(\beta_0 + \beta_1x_1 + \beta_2x_2 + ...\) . This is identical to everything we have done for three weeks
  3. The link function - the bit of maths that connects the two, so that a straight line on one scale becomes a sensible prediction on the other

OLS is just the GLM where the distribution is Normal and the link function does nothing at all.

The same model, three times

Written out in full, with the right-hand side identical every time. Only the left changes.

Ordinary regression — the outcome goes in as it is:

\[y_i = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \epsilon_i\]

Logistic regression — the probability goes in, through the logit:

\[\log\!\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots\]

Poisson regression — the average count goes in, through a log:

\[\log\big(E[y_i]\big) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots\]

\(E[y_i]\) is standard notation for the expected value of \(y\) — the average count for a school like this one.

Notice what is missing from the bottom two: the \(\epsilon\). OLS adds an error term on the end; a GLM does not, because the randomness is already described by the distribution you chose. That is the swap — an error term for a distribution.

Can we predict School Ofsted ratings?

Our question for today

  • We have used Ofsted ratings as a predictor and as a grouping variable already
  • Today we turn it around and make it the outcome
  • Given everything else we know about a school - its results, its intake, its absence rates, its size - can we predict what Ofsted said about it?
  • And if we can, what does that tell us about what Ofsted is actually measuring?
  • Note, Ofsted has changed how it rates schools now, so single-word ratings are no longer given

The outcome variable

3,170 schools with complete data. Note immediately how lopsided this is - we will come back to that, because it causes more trouble than anything else today.

First attempt: just use OLS?

  • Let’s re-code our ratings: “below Good” as 1 and “Good or better” as 0
  • We can regress it on Progress 8 exactly as we have been doing for the last three weeks.
lpm <- lm(below_good ~ P8MEA, data = model_data)
summary(lpm)

Call:
lm(formula = below_good ~ P8MEA, data = model_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.05109 -0.19520 -0.10193  0.03523  1.13124 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.145824   0.005812   25.09   <2e-16 ***
P8MEA       -0.274324   0.011145  -24.61   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3271 on 3168 degrees of freedom
Multiple R-squared:  0.1606,    Adjusted R-squared:  0.1603 
F-statistic: 605.9 on 1 and 3168 DF,  p-value: < 2.2e-16

The model runs. It produces coefficients. It produces an R². Nothing errors.

And its fitted values run from -0.504 to 1.051.

A probability of minus fifty percent?

  • The straight line gives 457 schools a negative probability of being rated below good.
  • At a Progress 8 of +1 it predicts a probability of being rated below good of about -13% (impossible).
  • The logistic curve cannot do this - it is bounded at 0 and 1 by construction.

Where the curve comes from

We need something that takes a straight line, which runs from \(-\infty\) to \(+\infty\), and squashes it into \([0,1]\).

  • Start with the probability \(p\) - bounded at 0 and 1, e.g. 0.4 (40%), awkward
  • Turn it into odds, \(\frac{p}{1-p} = \frac{probability of happening}{probability of NOT happening}\) or \(\frac{0.4}{0.6} = 0.67\) - now bounded at 0 and \(+\infty\), better
  • Take the log of that (\(\log\left(\frac{0.4}{0.6}\right) = -0.40\)) - now unbounded in both directions 🎉

\[\log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots\]

  • That left-hand side is the logit, and it is the link function. The right-hand side is exactly the model we have been fitting since week 5.
  • A logit simply is a mathematical function that converts a probability (between 0 and 1) into a log-odds value that can range from \(-\infty\) to \(+\infty\)
  • The link-function converts the expected value of the outcome (the mean) into something a straight line can cope with

Enough theory — let’s fit one

The simplest possible version: is a school rated below Good, predicted from Progress 8 and nothing else.

glm1 <- glm(below_good ~ P8MEA, data = model_data, family = binomial)
coef(summary(glm1))
             Estimate Std. Error   z value      Pr(>|z|)
(Intercept) -2.323537 0.07385577 -31.46047 3.019107e-217
P8MEA       -2.698673 0.13672741 -19.73761  1.025074e-86

The Std. Error, z value and Pr(>|z|) columns do exactly what their equivalents did in OLS — how precise the estimate is, and whether it is distinguishable from zero. It is the Estimate column that needs new thinking.

So what is that coefficient?

Put the estimates back into the equation from the last slide:

\[\log\left(\frac{p}{1-p}\right) = -2.32 + (-2.70 \times \text{Progress 8})\]

The left-hand side is log-odds. So whatever is multiplying Progress 8 must be in log-odds too:

-2.70 is the amount the log-odds change for each extra Progress 8 point.

One more point of Progress 8, and the log-odds of being rated below Good drop by 2.70.

And the intercept is just the same thing at zero: a school with a Progress 8 of exactly 0 has log-odds of -2.32.

Both of those sentences are true, and both are useless. Nobody has an intuition for log-odds — you cannot tell a head teacher their log-odds fell by 2.70. The next few slides are entirely about turning this number into something you can say out loud.

So we undo the log

We only took a logarithm in the first place to give the straight line somewhere to live. It was a convenience for the maths — nothing obliges us to report our answer on that scale.

So undo it. The opposite of a log is \(\exp\), and once you exponentiate a coefficient it stops adding and starts multiplying.

That one move — exp() — turns a number nobody can read into a sentence you could say to a head teacher.

Here is the recipe.

A recipe for reading your model

You have run glm(..., family = binomial) and you are staring at a table of coefficients. Four steps, every time:

  1. Ignore the raw coefficients. They are log-odds. Nobody can read them, including me
  2. Exponentiate them — exp(coef(model)). Each one is now an odds ratio
  3. Read each as a multiplier, for a one unit increase in that predictor:
    • above 1 → the odds go up
    • below 1 → the odds go down
    • exactly 1 → that predictor is doing nothing
  4. Say it in a sentence. “One more unit of x multiplies the odds of the outcome by [OR].” If the OR is below 1, a cut is easier to hear: \(1 - \text{OR}\)

That is the whole job. Everything else on the next few slides is just this recipe, done on our data.

Step 2 and 3, on our model

exp(coef(glm1))
(Intercept)       P8MEA 
 0.09792662  0.06729478

The one we care about is Progress 8: 0.067. It is below 1, so more progress means lower odds of being rated below Good.

Step 4 — the sentence:

One more Progress 8 point multiplies the odds of being rated below Good by 0.067.

Or, since that is under 1, say it as a cut: \(1 - 0.067 = 0.93\) —

One more Progress 8 point cuts those odds by about 93%.

Why it is called a ratio

Because that is literally what it is. Work out the odds for two schools, and divide. Our model, from a few slides back:

\[\log(\text{odds}) = -2.32 + (-2.70 \times x)\]

where \(x\) is the school’s Progress 8 score — the only thing that changes from school to school. −2.32 is the intercept and −2.70 is the coefficient on Progress 8; both stay put.

So put three values of \(x\) through it, one step at a time:

\(x\) (Progress 8) coefficient \(\times\ x\) log-odds \(=\) intercept \(+\) (coefficient \(\times\ x\)) odds \(= e^{\text{log-odds}}\)
−1 −2.70 × (−1) = +2.70 −2.32 + (+2.70) = 0.38 1.4552
0 −2.70 × (0) = 0.00 −2.32 + (0.00) = −2.32 0.0979
+1 −2.70 × (+1) = −2.70 −2.32 + (−2.70) = −5.02 0.0066

And now divide

Take the odds column from that table and divide each row by the one above it:

\[\frac{0.0979}{1.4552} = 0.067 \qquad \frac{0.0066}{0.0979} = 0.067\]

The same number both times. Going from −1 to 0 multiplies the odds by 0.067, and so does going from 0 to +1.

And 0.067 is exactly \(e^{-2.70}\) — the exponentiated coefficient. You never need the table. The coefficient already is the ratio, and it is the same one wherever you start.

A warning about that intercept

With one predictor, \(e^{\beta_0}\) was genuinely useful — the odds for a school on average progress. That worked because Progress 8 is already centred: zero means expected progress, not no progress.

Add more predictors and the intercept becomes the odds when every single predictor is zero at once. In our seven-predictor model that is a school with:

  • no disadvantaged pupils, no absence at all, no persistent absence
  • no high prior attainers
  • and log(cohort) = 0, meaning a cohort of one pupil

No such school exists. In our data, 0 schools have a cohort of one and 0 have zero absence. The intercept there — \(e^{-5.09}\) = 0.0062 — is a mathematical anchor, not a description of anything.

If you want an interpretable intercept, centre your predictors — subtract each one’s mean, so zero means average rather than none. Then \(e^{\beta_0}\) is the odds for a thoroughly average school, which is a sentence worth saying. Progress 8 arrives centred already, which is why it worked here for free.

And if you want a probability

Odds ratios are for comparing. If you want a plain answer for one school — what are the chances? — do not do it by hand. Ask the model:

predict(glm1,
        newdata = data.frame(P8MEA = c(-1, 0, 1)),
        type    = "response")
          1           2           3 
0.592699378 0.089192317 0.006546807

type = "response" means “answer on the scale of the thing I measured”. Here that is a probability between 0 and 1 that the school is rated below Good. Leave it out and R gives you log-odds — which is where we came in.

The three numbers are the three schools you asked about, in the order you asked for them:

P8MEA Which school? log-odds odds probability
−1 one point below the national average 0.375 1.4552 59.3%
0 exactly the national average -2.324 0.0979 8.9%
+1 one point above the national average -5.022 0.0066 0.7%

Read the last column only. The two middle ones are the scales we passed through on the way — and type = "response" is what saves you from making that journey by hand.

The reason we report ratios, not probabilities

Look down the two right-hand columns of that table again.

Going from The odds The probability
−1 → 0 ×0.067 59.3% → 8.9%
0 → +1 ×0.067 8.9% → 0.7%

The odds fall by the same multiplier every time — that multiplier is our odds ratio, and it is the same whatever school you start from.

The probabilities do not. The first step costs 50 percentage points, the second costs 8. Same one-point change in Progress 8, wildly different effect on the probability.

That is why we report an odds ratio and never “the probability changes by X%”. There is no single such number — it depends entirely on where the school started.

More than one predictor

Same recipe, more predictors

One predictor answered “what does Progress 8 do?”. A real question needs controls — what does Progress 8 do, for schools otherwise alike?

Same formula interface you have used for three weeks. Two differences: glm() not lm(), and a family argument.

logit_below <- glm(
  below_good ~ P8MEA + PTFSM6CLA1A + PERCTOT + PPERSABS10 +
               log(TPUP) + PTPRIORHI + MINORGROUP,
  data   = train,
  family = binomial
)
printCoefmat(coef(summary(logit_below)))
                              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)                 -5.0873362  1.0287222 -4.9453 7.603e-07 ***
P8MEA                       -2.2525783  0.2254148 -9.9930 < 2.2e-16 ***
PTFSM6CLA1A                 -0.0190684  0.0062864 -3.0333  0.002419 ** 
PERCTOT                      0.0717345  0.0798043  0.8989  0.368716    
PPERSABS10                   0.0477336  0.0249613  1.9123  0.055838 .  
log(TPUP)                    0.2654518  0.1759237  1.5089  0.131324    
PTPRIORHI                   -0.0055148  0.0097352 -0.5665  0.571066    
MINORGROUPMaintained school -0.0405095  0.1704089 -0.2377  0.812099    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(That is the coefficients block of summary(logit_below) — which also prints the call, the deviance and the AIC. We come back to those later.)

Eight numbers, all in log-odds, all unreadable. Step 1 of the recipe: ignore them.

Step 2: exponentiate the lot

exp(coef(logit_below))
                (Intercept)                       P8MEA 
                0.006174446                 0.105127821 
                PTFSM6CLA1A                     PERCTOT 
                0.981112244                 1.074370086 
                 PPERSABS10                   log(TPUP) 
                1.048891173                 1.304020028 
                  PTPRIORHI MINORGROUPMaintained school 
                0.994500370                 0.960300013

Step 3 — sort them by which side of 1 they fall:

Odds ratio Which means
Progress 8 0.105 well below 1 — odds drop sharply
% disadvantaged 0.981 just below 1
% persistently absent 1.049 just above 1 — odds rise
the rest ≈ 1 nothing detectable

Step 4: say them out loud

Progress 8. One more point multiplies the odds of being rated below Good by 0.105 — a cut of about 89%.

Persistent absence. One more percentage point multiplies those odds by 1.049 — a rise of about 5%.

Small ratios compound. 1.049 sounds trivial, but it applies per percentage point. Ten points worse is 1.04910 = 1.61 — the odds over half as big again.

And note Progress 8 has changed. On its own it was 0.067; here it is 0.105. Same arithmetic, different claim — this one holds absence, disadvantage, size and school type constant. Always say which one you mean.

Reading odds ratios

Anything whose interval crosses 1 is doing nothing detectable. Note how much work Progress 8 is doing compared with everything else.

How do you score a model like this?

For three weeks we judged models with R² — the share of the variation in the outcome the model managed to explain.

That worked because the outcome was a number, with variation in it to explain. Ours is a yes or no. There is no spread to account for, and R² as you know it simply does not exist here.

So we ask two questions instead, and they are genuinely different:

  1. Does the model fit the data? — likelihood, deviance, AIC, the pseudo-R²s, Hosmer–Lemeshow
  2. Does it make good decisions? — accuracy, the confusion matrix, ROC and AUC

A model can do respectably on one and badly on the other. We take them in turn.

1. Does the model fit?

Likelihood: the probability of what happened

Every school in the data either was or was not rated below Good — no ambiguity. And the model gives every school a probability. So we can ask one question of each school in turn:

What probability did the model put on the thing that actually happened?

School What actually happened Model’s P(below Good) Probability it put on what happened
A Below Good 0.80 0.80
B Good or better 0.05 0.95
C Below Good 0.03 0.03

The last column is the whole idea. For A and B the model was confident and right, so it put a high probability on the truth. For C it was confident and wrong — it gave that school only 0.03 of being below Good, and below Good is exactly what it was.

From three schools to one number

Do that for every school, then multiply the last column together. That product is the likelihood of the whole dataset.

And now the definition of fitting falls out of it:

Fitting the model is choosing the coefficients that make that number as large as possible.

That is all maximum likelihood means. There is no formula to solve as there was in OLS — the computer searches until it cannot make the number any bigger.

We report its log, for a blunt practical reason: multiplying 2,219 numbers that are all smaller than 1 gives a result far too small for a computer to hold. Taking logs turns that product into a sum.

logLik(logit_below)
'log Lik.' -717.1662 (df=8)

Always negative, and closer to zero is better. You will also hear this described as how surprised the model is: school C above was a nasty surprise, and surprises are what drag the number down.

Deviance: distance from a perfect fit

A log-likelihood on its own is hard to read. Deviance rescales it into something more like an error score — how far this model sits from one that predicted every school perfectly.

\[\text{deviance} = -2 \times \text{log-likelihood}\]

What it is
Null deviance 1883.4 with no predictors at all — just the overall rate
Residual deviance 1434.3 with our seven predictors in 
logit_below$null.deviance    # no predictors
[1] 1883.425
logit_below$deviance         # our model
[1] 1434.332

Lower is better, and the drop — 449 — is what the predictors bought us. This is the direct equivalent of the sum-of-squares drop you used to read off an OLS table.

AIC: fit, minus a fine for complexity

There is a problem with judging a model on deviance alone: adding any predictor, however useless, will always lower it slightly. Left to itself, deviance always says “more variables”.

So we charge for them:

\[\text{AIC} = \text{deviance} + 2k \quad\text{where } k = \text{number of coefficients}\]

For our model: 1434.3 + 2 × 8 = 1450.3. A predictor only improves the AIC if it reduces the deviance by more than the 2 points it costs.

AIC(logit_below)
[1] 1450.332

AIC has no meaning on its own. 1450 is neither good nor bad — it is only useful for comparing models fitted to the same data. Lower wins. That is exactly how we used it earlier to choose between a raw and a logged predictor.

The pseudo-R²s

People badly want an R², so several substitutes exist. The two you will meet in software output:

Cox & Snell — built from the likelihood ratio between your model and the null model. For us: 0.183.

That looks dismal — until you learn its catch. Cox & Snell cannot reach 1. With a binary outcome its ceiling here is 0.572, so 0.183 is really 0.183 out of a possible 0.572.

Nagelkerke fixes precisely that, by dividing by the ceiling:

\[R^2_{\text{Nagelkerke}} = \frac{0.183}{0.572} = \mathbf{0.320}\]

Now it can reach 1, and it reads much more like the R² you are used to.

Getting them in R

Unlike logLik(), AIC() and the deviances, the pseudo-R²s are not in base R — summary() will not give you one. The performance package will:

performance::r2_coxsnell(logit_below)
Cox & Snell's R2 
       0.1832199 
performance::r2_nagelkerke(logit_below)
Nagelkerke's R2 
      0.3202809

So the whole fit toolkit is four lines:

What How
Log-likelihood logLik(model)
Deviances model$null.deviance, model$deviance
AIC AIC(model)
Pseudo-R² performance::r2_nagelkerke(model)

summary(model) prints the deviances and the AIC for you at the bottom of its output — you have been looking at them all along without knowing what they were.

…and the same four in Python

In statsmodels the first three are attributes of the fitted model rather than functions — no brackets:

What R Python (statsmodels)
Log-likelihood logLik(model) model.llf
Deviances model$null.deviance / model$deviance model.null_deviance / model.deviance
AIC AIC(model) model.aic
Pseudo-R² performance::r2_nagelkerke(model) nothing built in

For the pseudo-R²s you write them out. Two lines, and they agree with performance to six decimal places:

n   = model.nobs
cs  = 1 - np.exp((model.llnull - model.llf) * 2 / n)   # Cox & Snell
nag = cs / (1 - np.exp(2 * model.llnull / n))          # Nagelkerke

One trap worth knowing about

Fit with smf.logit() instead of smf.glm() and Python does hand you a .prsquared. It is tempting. Do not assume it is the one we just calculated.

.prsquared is McFadden’s pseudo-R² — a different measure, built a different way. On our model:

Nagelkerke 0.320
McFadden (.prsquared) 0.238

Same data, same model, same fit — different number, because they are different quantities with the same nickname.

Whichever you report, name it. “Pseudo-R² = 0.3” tells your reader almost nothing on its own.

…and how much to trust them

Report them if the journal or the marker expects them — but know what you are saying.

  • They are not “variance explained”. Nothing here is variance. The name is borrowed, and the interpretation does not come with it
  • They are not comparable across datasets, and barely across models. A Nagelkerke of 0.32 is not “half as good” as one of 0.64
  • There is no accepted threshold for a good one. Values that would be embarrassing in OLS are routine here

Treat a pseudo-R² as a rough indication that the predictors are doing something, and nothing more. If you want to know whether the model is actually any use, ask the second question — the one we come to next.

2. Does it make good decisions?

Pause — this part is not just for logistic regression

Everything in this section applies to any model you will ever fit — the OLS models from Vol 1 and 2, the mixed effects models from Vol 3, and every machine learning method you go on to meet.

What transfers, unchanged:

  • Hold data back. Fit on some, judge on the rest. A model marked on its own training data flatters itself, whatever kind of model it is
  • Compare against a baseline. Never ask “is this good?” — ask “is this better than the dumbest thing I could have done?”
  • Fit and prediction are different questions. A model can describe this data well and still be useless on the next school

…but the scoreboard changes

What does not transfer is how you keep score. Accuracy, the confusion matrix and AUC all need a categorical outcome — none of them means anything if you are predicting a number.

Swap the outcome type and you swap the metric:

If the outcome is… Baseline to beat Typical metrics
A category always predict the commonest class accuracy, AUC, precision/recall
A number always predict the mean RMSE, MAE, R² on held-out data
A count always predict the mean count RMSE, or deviance on held-out data

So learn the shape of what follows, not just the logistic version of it. We come back to what this means for the earlier lectures at the end.

A recipe for judging the decisions

Same shape as the fit recipe. Six steps, and the order matters:

  1. Hold data back — fit on some schools, judge on the rest
  2. Predict on the ones it never saw
  3. Compare against a baseline — never quote accuracy on its own
  4. Read the confusion matrix — which kind of mistake is it making?
  5. Check the ranking with AUC — this needs no cut-off at all
  6. Only then pick a cut-off — and that is a policy decision, not a statistical one

Steps 5 and 6 are the ones people skip, and they are where the interesting part lives.

First, the thing we are guarding against

Here are fifteen points, and two models fitted to them. The red one is better at its job — it passes through every single point and makes no error at all.

On the left the red model is flawless and the teal one is not. On the right — the same two models, fifteen points neither has seen — the red one is now the worse of the two, by a distance.

Overfitting

The red model did not learn the pattern. It learned the noise — the particular wobbles of the particular fifteen points it was handed, which the next fifteen do not share.

And here is the uncomfortable part. You can always make the fit better. Add a predictor, add an interaction, add a squared term — the fit to the data in front of you never gets worse. Push it far enough and you reach error 0.00, like the red line.

So “it fits well” is not evidence of very much on its own — and a perfect fit is not a triumph, it is a warning sign.

The only honest test is data the model has never seen. Which is exactly what we do next.

Steps 1 and 2: don’t mark your own homework

So we split the schools before fitting anything: learn from one part, judge on the other.

set.seed(42)
train_idx <- sample(nrow(model_data), floor(0.7 * nrow(model_data)))
train <- model_data[train_idx, ]   # fit the model on these
test  <- model_data[-train_idx, ]  # judge it on these

2,219 schools to learn from. 951 kept back to test on.

Test-set accuracy is then one question: “of the 951 schools the model never met, what share did it guess were below good, correctly?

Step 3: from a probability to a yes or a no

The model does not hand back a yes or a no. For each school it hands back a probability — a number between 0 and 1.

To count how often it was right, we first have to turn each of those into an answer. That needs a cut-off.

We will use 0.5, which is just this rule:

If this school’s probability is above 0.5, call it below Good. Otherwise call it Good or better.

0.5 is not “half the time”. It is a line drawn on each school’s own probability, one school at a time — and most schools are nowhere near it.

Applied to all 951 test schools, it flags 55 of them. Not 476.

So, one school at a time

Four of our 951 test schools, taken through the rule:

School Model’s probability Above 0.5? So we predict… It actually was… Right?
A 0.03 no Good or better Good or better ✔
B 0.70 yes Below Good Below Good ✔
C 0.30 no Good or better Below Good ✘
D 0.59 yes Below Good Good or better ✘

C was missed and D was a false alarm. Each school is judged on its own probability — the cut-off simply sorts them either side of the line, one at a time.

Now do that for every school and count the ticks.

accuracy <- mean((p_below > 0.5) == test$below_good)
baseline <- max(mean(test$below_good), 1 - mean(test$below_good))

820 ticks out of 951 — test-set accuracy 86.2%. Nearly nine schools in ten.

Now watch.

The single most important slide today

A model that makes the same prediction every single time scores 85.3%. Ours beats it by 0.9 percentage points.

And that red bar is a model you have already met. glm(below_good ~ 1) fits an intercept and nothing else, so it predicts the base rate — 0.151 — for every school. Under 0.5, so it answers “not below Good” to all of them. Its deviance is the null deviance (1883) from the fit section: R has been fitting this model for you every time you call glm().

Step 4: why accuracy lied to us

Teal = right, red = wrong. Both models are judged on the same 951 test schools.

  • The baseline never flags a single school — its entire left-hand column is zero — and still scores 85.3%
  • Ours flags 55 schools and gets 32 of them right. That is the whole of the 0.9-point gain
  • It still misses 108 of the 140 really below Good — close to useless for deciding where to send inspectors at this threshold

Step 5: judge the ranking, not the labels

Everything so far hung on 0.5 — and nothing made 0.5 the right cut-off. Move it and the accuracy, the catches and the false alarms all move with it.

So ask a question that does not depend on the cut-off at all:

Does the model put schools in the right order? Are the struggling ones nearer the top of the list than the healthy ones?

A model that ranks well is useful whatever cut-off you eventually pick. That is what ROC and AUC measure.

What a ROC curve actually is

ROC stands for Receiver Operating Characteristic — a name that makes no sense until you know where it came from.

It was invented by radar engineers during the Second World War. The receiver was the radar set, and the problem was exactly ours: an operator watching a screen has to decide what is a real aircraft and what is noise. Catch more aircraft and you also raise more false alarms. Where you set that dial is the operating characteristic.

Same idea here. Slide the cut-off from 1 all the way down to 0, and at every setting write down two numbers:

  • How many struggling schools did we catch? (of those that really are below Good — we want this high)
  • How many healthy schools did we wrongly flag? (of those that are fine — we want this low)

Plot one against the other and you get the curve. Each red dot is one cut-off. The dashed diagonal is what pure guessing would give you — the further above it the curve sits, the more the model knows.

AUC — one number for the whole curve

A curve is awkward to quote in a sentence, so we summarise it by the area underneath it — which is all AUC stands for: Area Under the Curve. You will also see it written AUROC, area under the ROC curve; same thing.

AUC What it means
0.5 the diagonal — no better than tossing a coin
0.7–0.8 useful
0.9+ very strong
1.0 perfect — every struggling school ranked above every healthy one 
roc <- performance::performance_roc(logit_below, new_data = test)
bayestestR::area_under_curve(roc$Specificity, roc$Sensitivity)
[1] 0.8183283

What our AUC says

Ours is 0.818, which reads as:

Pick one school that really is below Good and one that really is not. The model gives the struggling one the higher predicted probability of being below Good 82% of the time.

“Predicted probability” means the model’s own output — the 0-to-1 number from predict(..., type = "response") that we have been cutting at 0.5. Not Progress 8, and not any of the predictors: a higher number means more likely to be below Good.

Remember accuracy said this model was barely better than useless. It ranks well; we were just reading it at the wrong cut-off.

Step 6: now choose your cut-off

The ROC curve showed you every cut-off at once. Choosing one just means deciding which point on that curve you want to sit at.

Here are the three red dots from that curve, written out as counts:

Flag if probability above Struggling schools caught …and missed False alarms
0.50 32 108 23
0.30 57 83 75
0.10 121 19 316

Lower the cut-off and you move up and to the right along the curve: more struggling schools caught, more healthy ones wrongly flagged. There is no point on the curve that is simply best.

The bit AUC cannot do for you

AUC told you the curve was worth picking from. A model with an AUC of 0.5 has nothing to offer at any cut-off. Ours ranks well, so some point on it is worth using.

But AUC cannot tell you where to sit on it, because that is not a statistical question:

  • If a missed school means children let down — catch more, and accept the false alarms.
  • If a false alarm means a wasted inspection — flag fewer, and accept the misses.

The model cannot tell you which of those matters more. You have to.

Wrapping up

The three you did not get today

Everything we have done generalises, and the hard part is behind you. Each of these is the same machinery with a different family:

If your outcome is… Use And the coefficients are
An ordered category — Ofsted’s four grades, a Likert scale MASS::polr() odds of landing in a higher category
An unordered category — which region, which destination nnet::multinom() odds against a reference category
A count — pupils excluded, schools below Good in an LA glm(family = poisson) rate ratios, read just like odds ratios

Two things to know if you reach for the count model:

  • An offset — offset(log(pupils)) — turns a model of counts into a model of rates, which is almost always what you want
  • Check for overdispersion (residual deviance ÷ residual df). Much above 1 and you want quasipoisson or MASS::glm.nb()

The practical covers these, and the reading goes further. You already have the two ideas that make them work: a link function, and coefficients you have to exponentiate before they mean anything.

The family, in one slide

Outcome family = Link Interpret coefficients as
Continuous gaussian identity Units of y per unit of x
Yes / no binomial logit Odds ratios
x out of n binomial logit Odds ratios
Ordered categories MASS::polr() logit Odds of being in a higher category
Unordered categories nnet::multinom() logit Odds vs a reference category
Count poisson log Rate ratios
Overdispersed count MASS::glm.nb() log Rate ratios

One framework. Change the family, change the link, and everything you learned about building and interpreting a linear model still applies.

What carried over, and what didn’t

Still true

  • The recipe: start simple, add complexity that earns its place
  • Confounding, mediation, interactions
  • Multicollinearity and VIF
  • Thinking hard about what belongs in the model

Different now

  • No R² - use AIC, deviance, or predictive performance
  • Coefficients need exponentiating before they mean anything
  • Residual plots look strange and need different diagnostics
  • Accuracy is a trap on imbalanced outcomes

Where you can go next

  • Generalised additive models (GAMs) - when the relationship is smooth but not a straight line
  • glmer() - the intersection of this week and last week: a GLM with random effects, for grouped non-normal data
  • Spatial and geographically weighted regression - when the relationship itself varies over space
  • Ridge and lasso - when you have too many predictors and need the model to choose

You now have the two big extensions to OLS: last week we relaxed independence, this week we relaxed normality. Almost everything else is a combination of the two.

Explaining, or predicting?

Worth noticing what changed across these four sessions — because it was not just the model.

Vols 1–3: explaining

R², residual plots, VIF, AIC, significance.

Every one of those is computed on the data you fitted to. The question is does this model describe the process I am looking at?

Today: predicting

Train/test split, accuracy against a baseline, AUC.

All computed on data the model has never seen. The question is will this work on the next school?

Both are legitimate. They are not the same job, and they can disagree — a model packed with significant coefficients can predict badly, and a model that predicts well can be uninterpretable.

Decide which one you are doing before you start, because it changes what counts as a good model. For the assessment, say which — “I am explaining why schools differ” and “I am predicting which schools will struggle” are different claims and need different evidence.

Shmueli, G. (2010) To Explain or to Predict? Statistical Science, 25(3), 289–310. doi:10.1214/10-STS330

So — can we predict School Ofsted ratings?

We asked that on slide one. The honest answer is partly, and the model is specific about which part.

Our model Verdict
Ranking — AUC 0.818 Good. It puts a struggling school above a healthy one 82% of the time
Labelling — accuracy at 0.5 86.2% vs a do-nothing baseline of 85.3% Almost worthless
Fit — Nagelkerke R² 0.320 (Cox & Snell 0.183) Respectable for cross-sectional social data

And on the variables: Progress 8 dominates everything else (OR 0.105), persistent absence adds a little per percentage point (OR 1.049), and most of the rest have intervals crossing 1 — doing nothing detectable.

In one line: a decent ranking device, and a poor labelling device. Those are different claims, and the deck spent an hour on why.

…and here is what it cannot tell you

  • Circularity. Inspectors read Progress 8 before they visit. We are not independently predicting a judgement — we are partly recovering one of its own inputs. Our strongest coefficient is the one we should trust least.
  • Timing. Ratings here span 2011–2024; the results are 2022-23. Some ratings were handed out a decade before the data we are “predicting” them with.
  • Causation. Nothing here says that raising a school’s Progress 8 would change its rating. We have an association, measured on log-odds.
  • Rarity. Only 15% of schools are below Good. The model inherits that, and so does every accuracy figure you quote from it.

So why does any of this matter in Brighton?

Volumes 1 to 3 asked what drives attainment, and the council’s answer was social mixing.

But families do not choose a school from a regression table. In Brighton they choose largely on two things: a school’s Ofsted rating and its raw Attainment 8.

So it is worth turning the question around and asking what drives the label — because the label is what actually moves people.

That is the question this whole lecture has been answering.

How much of the rating is just the intake?

Strip the model back and watch what survives:

What the model is allowed to see AUC
Everything, including Progress 8 0.818
Drop Progress 8, keep absence 0.747
Intake only — disadvantage, size, prior attainment, school type 0.656

Most of the ranking does come from Progress 8 — the thing inspectors read before they visit. But intake alone still ranks at 0.656, against 0.5 for a coin toss.

That is modest, and it is real: before a school does anything at all, who walks through its doors already shifts the odds on its Ofsted rating.

Which means choosing “a school with a Good rating” is partly choosing an intake. Advantaged intakes make favourable labels slightly more likely, the labels attract advantaged families, and the sorting repeats.

That is the mechanism the admissions policy in Volume 1 was trying to interrupt — and we have arrived at it from completely the other end.

Careful: this is association, not causation, and it is emphatically not a claim that inspectors are getting it wrong.

The bit the software will never do for you

Every number on the last two slides is correct. The code ran, the maths is right, the AUC is genuinely 0.82.

None of it tells you whether the model is any use. That verdict came from knowing something R does not: how an Ofsted inspection actually works.

These methods will carry you a long way, and they will show you things you could not have seen by eye. But they measure associations in the data you happened to collect — they cannot tell you where that data came from, who made it, or what it would mean to act on it.

That part is yours as an expert analyst.

Conclusions

  • Look at your outcome variable before you choose a model. Week 1’s data types were not an academic exercise - they determine what you are allowed to do
  • A GLM is not a new method, it is the same linear predictor with a different distribution and a link function bolted on
  • Always compare against a baseline. A model that beats “predict the majority every time” by half a percentage point has not learned much
  • Check your assumptions, and remember the threshold is a choice you make rather than one the model makes for you
  • And ask where your predictors came from. Our best one turned out to be something Ofsted had already read

Today’s Practical

  • Building and interpreting your own logistic and Poisson models on the schools data
  • Practice at the bit that matters most: interpreting odds ratios and rate ratios in plain English
  • And deciding, for an outcome of your own choosing, which member of the family you actually need