Chapter 9 showed that for a linear model the Cox-Snell formula recovers ordinary \(R^2\) exactly, while the naive proportional reduction in deviance does not. This appendix follows that up: what happens to these measures once the outcome is no longer continuous, and why the competing versions disagree.
For models fitted by maximum likelihood there is no sum of squared error to reduce, so \(R^2\) has to be rebuilt out of deviance. Different people rebuilt it differently.
Proportional reduction in deviance (McFadden’s), which mimics PRE directly:
\[\text{PRD} = \frac{D_0 - D_1}{D_0}\]
Cox and Snell’s, which comes from the likelihood ratio and rescales by sample size:
\[R^2_{CS} = 1 - \exp\left(\frac{-(D_0 - D_1)}{n}\right)\]
Before asking how these behave on a binary outcome, it helps to see them where we already know the right answer.
states <- as.data.frame(state.x77) |>
rename(life_exp = `Life Exp`, hs_grad = `HS Grad`)
mc <- lm(life_exp ~ 1, data = states)
ma <- lm(life_exp ~ hs_grad, data = states)
D0 <- -2 * as.numeric(logLik(mc))
D1 <- -2 * as.numeric(logLik(ma))
c(PRD = (D0 - D1) / D0,
Cox_Snell = 1 - exp(-(D0 - D1) / nobs(ma)),
ordinary_R2 = summary(ma)$r.squared) PRD Cox_Snell ordinary_R2
0.1215190 0.3389757 0.3389757 Cox-Snell is exactly \(R^2\). PRD is not close.
Cox-Snell is a genuine generalization of \(R^2\), and it reduces to it when the model is linear. PRD is a different construction that happens to resemble PRE, and for a linear model it is simply wrong. As Chapter 9 explains, the trouble is the \(\sigma\) estimate: as RSS goes to zero, \(\sigma\) goes to zero and the likelihood diverges, so the ratio of deviances doesn’t behave.
Neither, though, will recover the \(R^2\) you would have obtained had a binary outcome not been dichotomized in the first place. That’s what the rest of this appendix shows.
We can construct data where the truth is known. Generate two variables with an exact correlation, so the latent \(R^2\) is fixed at a value we choose, then dichotomize the outcome and see what the measures report.
# create dataset with exact known correlation
# convert outcome to binary and see if PRD = Cox-Snell
get_stats <- function(n, r) {
Sigma <- matrix(c(1, r,
r, 1), 2, 2)
mat <- MASS::mvrnorm(n, mu = c(0, 0), Sigma = Sigma, empirical = TRUE)
x = mat[,1]
ystar = mat[,2]
y = if_else(ystar > 0, 1, 0) # binom transform
m0 <- glm(y ~ 1, family = binomial())
m1 <- glm(y ~ x, family = binomial())
d0 <- as.numeric(-2 * logLik(m0))
d1 <- as.numeric(-2 * logLik(m1))
delta <- d0 - d1
prd <- (d0 - d1) / d0
cs <- 1 - exp(-delta / length(x))
tmp <- tibble::tibble(
cs = cs,
prd = prd
)
return(tmp)
}
set.seed(522)
# for a given PRE (r2), get the data and compare stats for increasing N
mytests <- expand_grid(
n = seq(50, 1000, 50),
r = c(sqrt(.1), sqrt(.5), sqrt(.8))) |>
rowwise() |>
mutate(stats = list(get_stats(n, r))) |>
unnest_wider(stats) |>
ungroup() |>
mutate(latent = factor(round(r^2, 1)))
mytests |>
select(n, latent, prd, cs) |>
head(4) |>
mutate(across(c(prd, cs), \(z) round(z, 4)))# A tibble: 4 x 4
n latent prd cs
<dbl> <fct> <dbl> <dbl>
1 50 0.1 0.0581 0.077
2 50 0.5 0.507 0.504
3 50 0.8 0.699 0.614
4 100 0.1 0.0265 0.036plt(cs ~ n | latent, data = mytests, type = "l", lwd = 2,
ylim = c(0, 0.85),
xlab = "Sample size", ylab = "Cox-Snell pseudo-R²")
abline(h = c(0.1, 0.5, 0.8), lty = 3, col = "gray60")
ggplot(mytests, aes(x = n, y = cs, color = latent)) +
geom_line(linewidth = 0.8) +
geom_hline(yintercept = c(0.1, 0.5, 0.8), linetype = "dotted",
color = "gray60") +
ylim(0, 0.85) +
labs(x = "Sample size", y = "Cox-Snell pseudo-R²", color = "Latent R²")
mytests |>
filter(n %in% c(50, 200, 500, 1000)) |>
select(latent, n, prd, cs) |>
mutate(across(c(prd, cs), \(z) round(z, 3)))# A tibble: 12 x 4
latent n prd cs
<fct> <dbl> <dbl> <dbl>
1 0.1 50 0.058 0.077
2 0.5 50 0.508 0.504
3 0.8 50 0.699 0.614
4 0.1 200 0.05 0.067
5 0.5 200 0.326 0.363
6 0.8 200 0.473 0.481
7 0.1 500 0.06 0.08
8 0.5 500 0.285 0.326
9 0.8 500 0.543 0.529
10 0.1 1000 0.045 0.06
11 0.5 1000 0.279 0.32
12 0.8 1000 0.553 0.535The two measures disagree with each other. At a latent \(R^2\) of 0.1, Cox-Snell runs about a third larger than PRD. At 0.8 the ordering reverses for most sample sizes, with PRD the larger of the two. There is no fixed relationship you can use to convert one into the other, and not even a consistent direction.
Neither recovers the latent \(R^2\). When the truth is 0.8, both report values in the neighbourhood of 0.5, well short of it. Dichotomizing an outcome throws information away, and no summary statistic can retrieve it. Note that this is a different failure from the linear-model case above: there, Cox-Snell was exactly right and PRD was wrong. Here Cox-Snell is behaving correctly as a likelihood-based measure; it is the dichotomization that has destroyed the information.
Cox-Snell cannot reach 1. Its maximum is bounded below one for a binary outcome, a defect that motivated yet another version (Nagelkerke’s), which rescales Cox-Snell by its own maximum. That fixes the ceiling and introduces its own oddities.
Report a pseudo-\(R^2\) if a reader expects one, and say which one it is. Never compare a pseudo-\(R^2\) against an ordinary \(R^2\) from a different model, and never compare pseudo-\(R^2\) values computed by different formulas.
For comparing models fitted to the same data, the information criteria of Chapter 10 are better behaved and easier to defend. For describing how much a predictor matters, the average marginal effects of Chapter 16 are more informative than any single-number summary of fit.