In Chapter 9 we saw that for a linear model, the Cox-Snell formula gives exactly the same answer as ordinary \(R^2\), but the naive proportional reduction in deviance doesn’t. In this appendix we’ll look at what happens to these measures when the outcome isn’t continuous, and why the different versions don’t agree.
For models fit by maximum likelihood, there’s no SSE, so we have to build something like \(R^2\) out of the deviance instead. Different people have done this in different ways.
Proportional reduction in deviance (McFadden’s) works just like PRE:
\[\text{PRD} = \frac{D_0 - D_1}{D_0}\]
Cox and Snell’s is based on the likelihood ratio and adjusts for the sample size:
\[R^2_{CS} = 1 - \exp\left(\frac{-(D_0 - D_1)}{n}\right)\]
Before we look at a binary outcome, let’s see what these do when 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 the same as \(R^2\), but PRD isn’t even close.
So Cox-Snell is a real generalization of \(R^2\). When the model is linear, it gives you \(R^2\). PRD just looks like PRE, and for a linear model it’s wrong. As we saw in Chapter 9, the problem is the \(\sigma\) estimate: as RSS goes to zero, \(\sigma\) goes to zero and the likelihood goes to infinity, so the ratio of deviances doesn’t work.
But neither one can tell you what the \(R^2\) would have been if a binary outcome hadn’t been dichotomized in the first place. That’s what the rest of this appendix shows.
We can make up data where we know the truth. We’ll generate two variables with an exact correlation, so we know the “latent” \(R^2\). Then we’ll dichotomize the outcome and see what the two measures say.
# 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 don’t agree. When the latent \(R^2\) is .1, Cox-Snell is about a third bigger than PRD. At .8, PRD is usually the bigger one. So you can’t convert one into the other, and one isn’t even consistently bigger than the other.
Neither one recovers the latent \(R^2\). When the truth is .8, both give values between about .5 and .7. When you dichotomize an outcome, you throw information away, and no summary statistic can get it back. (This is a different problem from the linear model above, where Cox-Snell was exactly right and PRD was wrong. Here, Cox-Snell is doing what it’s supposed to do. It’s the dichotomization that lost the information.)
Cox-Snell can’t reach 1. For a binary outcome, its maximum is less than 1. That’s why there’s yet another version (Nagelkerke’s), which rescales Cox-Snell by its maximum. But that has its own problems.
If your readers expect a pseudo-\(R^2\), go ahead and report one, but say which one it is. Don’t compare a pseudo-\(R^2\) to an ordinary \(R^2\) from a different model, and don’t compare pseudo-\(R^2\) values from different formulas.
To compare models fit to the same data, I’d use the information criteria from Chapter 10 instead. And to describe how much a predictor matters, the average marginal effects from Chapter 16 tell you a lot more than any single measure of fit.