Query regarding the calculation to Nagelkerke's R2

[ianchlee]

(Moved from issue #852 to discussion)

I am trying to use Nagelkerke's R2 for some general linear models (negbinom) using performance::r2_nagelkerke.
The formula to Cox & Snell's R2 ($R^2_{CS}$) and Nagelkerke's R2 ($R^2_N$) are as follows (ref):

$R^2_{CS} = 1 - \left( \frac {\text{L}(M_{full})} {\text{L}(M_{null})} \right) ^{\frac {2} {n}} = 1 -\text{exp} \left(\frac{2}{n} \times \left(\text{logLik}(M_{full}) - \text{logLik}(M_{null})\right)\right)$

$R^2_N = \frac {1 - \left( \frac {\text{L}(M_{full})} {\text{L}(M_{null})} \right) ^{\frac {2} {n}}} {1 - \text{L}(M_{null}) ^ {\frac {2} {n}} } = \frac {R^2_{CS}} {1 - \text{exp} \left(\frac{2}{n} \times \text{logLik}(M_{null})\right) } $

Based on the source code to r2_nagelkerke.glm, I can benchmark the following calculation:

library(performance)
library(insight)

# create model glm
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm_dt <- data.frame(treatment, outcome, counts) # showing data

glm_model <- glm(data = glm_dt, counts ~ outcome, family = poisson())
# null model, same as null_model(glm_model)
model0 <- glm(data = glm_dt, counts ~ 1, family = poisson())

n <- n_obs(glm_model, disaggregate = TRUE)
model_dev <- glm_model$deviance
null_dev <- glm_model$null.deviance

cs_r2 <- ( 1 - exp( (model_dev - null_dev) / n ) )
r2_coxsnell(glm_model) == cs_r2
# > true

r2_nag_func <- r2_nagelkerke(glm_model)
r2_nag_calc <- cs_r2 / ( 1 - exp(-null_dev / n) )

r2_nag_func == r2_nag_calc
# > true

I understand that the deviance is different from log likelihood and can be calculated from difference of log-likelihood (see link), with
$\text{dev} = 2 \times \left( \text{logLik}(M_{saturated}) - \text{logLik}(M_{proposed}) \right)$

null_dev == -2 * logLik(model0) 
# > false
model_dev == -2 * logLik(glm_model) 
# > false

For Cox & Snell's R2 (CS R2), the result should be the same, as the common term of $\text{logLik}(M_{saturated})$ is cancelled:

D.LL_glm <- as.numeric(-2 * logLik(glm_model))
D.LL_null <- as.numeric(-2 * logLik(model0))

all.equal(cs_r2,  1 - exp( (D.LL_glm - D.LL_null) / n) )
# > true

However, with Nagelkerke's R2, wouldn't dividing the$R^2_{CS}$ with $1 - \text{exp} \left(- \frac{\text{dev} (M_{null}) } {n} \right)$ be different from $1 - \text{exp} \left(\frac{2}{n} \times \text{logLik}(M_{null})\right)$?

[ianchlee]

I've made some digging in other threads, and from #165 and #300, it seems the key issue is in the calculation of the log likelihood with frequency models (such as log binomial, poisson, negative binomial...):

set.seed(101)
n <- 200
counts <- sample(50:90, size = n, replace = T)
class <- gl(5,1,n)
treatment <- gl(4,50)
freq_dt <- data.frame(treatment, class, counts)
freq_dt$total <- sample(100:110, size = n, replace = T)
freq_dt$fail <- freq_dt$total - freq_dt$counts

### log-binomial model formulation 1 (cbind)
logbin.mod <- glm(data = freq_dt, cbind(counts, fail) ~ treatment + class, family = binomial(link = "log"))
logbin.null <- insight::null_model(logbin.mod)
n <- insight::n_obs(logbin.mod)

r2_nagelkerke(logbin.mod)
#> null (can't calculate accurate R2 for binomial models that are not Bernoulli models)
(1 - exp(-2/n * (logLik(logbin.mod) - logLik(logbin.null)))) / (1 - exp(2/n * logLik(logbin.null)))
#> 0.3896144

### log-binomial model formulation 2 (frequencies + total)
logbin.mod2 <- glm(data = freq_dt, counts/total ~ treatment + class, weights = total, family = binomial(link = "log"))
logbin.null2 <- insight::null_model(logbin.mod2)

r2_nagelkerke(logbin.mod2)
#> null (can't calculate accurate R2 for binomial models that are not Bernoulli models)
(1 - exp(-2/n * (logLik(logbin.mod2) - logLik(logbin.null2)))) / (1 - exp(2/n * logLik(logbin.null2)))
#> 0.3896144

### poisson model (counts + offset)
poisson.mod <- glm(data = freq_dt, counts ~ treatment + class + offset(log(total)), family = poisson)
poisson.null <- insight::null_model(poisson.mod)

r2_nagelkerke(poisson.mod)
#> 0.1807906
(1 - exp(-2/n * (logLik(poisson.mod) - logLik(poisson.null)))) / (1 - exp(2/n * logLik(poisson.null)))
#> 0.1584426

### Negative binomial model (counts + offset)
nb.mod <- MASS::glm.nb(data = freq_dt, counts ~ treatment + class + offset(log(total)))
nb.null <- MASS::glm.nb(data = freq_dt, counts ~ offset(log(total))) # note insight::null_model(nb.mod) throws an error

r2_nagelkerke(nb.mod)
#> 0.1304831
(1 - exp(-2/n * (logLik(nb.mod) - logLik(nb.null)))) / (1 - exp(2/n * logLik(nb.null)))
#> 0.08321774

r2_nagelkerke gives the same result when it is a Bernoulli's model:

set.seed(101)
n <- 200
x <- rnorm(n)
a <- 1
b <- -2
p <- exp(a + b * x) / (1 + exp(a + b * x))
y <- factor(ifelse(runif(n) < p, 1, 0), levels = 0:1)
d <- data.frame(x, y)

mod.bern <- glm(y ~ x, data = d, family = binomial())
bern.null <- insight::null_model(mod.bern)
n <- insight::n_obs(mod.bern)

r2_nagelkerke(mod.bern)
#> 0.5852906
(1 - exp(-2/n * (logLik(mod.bern) - logLik(bern.null)))) / (1 - exp(2/n * logLik(bern.null)))
#> 0.5852906

So can I ask for frequency based model (in particular negative binomial), how should the log likelihood and R2 be calculated? Should I use the suggested offset of log (N choose k) in #165 ?