5
$\begingroup$

I am working on my Master's Thesis and I was fitting various glms on my data; and since I can't calculate adjusted R² values for my models, I opted for Nagelkerke's pseudo R². I used the rcompanion package and the nagelkerke() function to calculate these values for my various models, that I had stored inside a list. I used a loop to calculate the Nagelkerke Pseudo R² for every model inside the list. I didn't notice anything weird there, no errors, no warnings, the values looked fine, so I trusted the results and proceeded.

Now I have a slightly different dataset (a different explanatory variable), I used the same loop and the Nagelkerke Pseudo-R² are ridiculously high, up to 1.8. I've read that this value should be between 0 and 1. So, I picked one of these models (the one that had a Pseudo R² of 1.8) and "hand-recalculated" the nagelkerke R² to see what went wrong. I fit the glm and saved it as a single object, not in a list of glms. When I put it into nagelkerke(), the value it gave me -0.0960369. Its better, but again out of the range(for comparison purposes, I used glm error distribution: gaussian, link function: identity, so I could compare that model and the Pseudo R² to the lm() and the adjusted R², which was -0.003449).

So I thought, well, maybe this package rcompanion also adjusts the Pseudo R² somehow? Hence, I tried to calculate the Nagelkerke Pseudo R² by hand:

testdata <-na.omit(testdata)

testmodel_lm <- lm(CSRatio21 ~ FML, data = testdata) #for comparison
testmodel <- glm(CSRatio21 ~ FML, data = testdata, family = gaussian(link="identity"))
summary(testmodel)
nullmodel <- glm(CSRatio21 ~ 1, data = testdata, family = gaussian(link="identity"))

# Number of data points
n <- nrow(testdata)

# residuals testmodel
residuals <- testdata$CSRatio21 - (0.0031229 * testdata$FML + 1.336591) #numbers are intercept and slope estimates from testmodel

# sum of squared residuals
sum_squared_residuals <- sum(residuals^2)

# residual standard deviation (sigma)
sigma <- sqrt(sum_squared_residuals / (n - 2))

# likelihood testmodel
likelihood_linear <- prod(1 / (sigma * sqrt(2 * pi)) * exp(-residuals^2 / (2 * sigma^2)))

loglikmodel <- logLik(testmodel) #second variant

# likelihood null model
null_residuals <- testdata$CSRatio21 - 1.432629
null_sum_squared_residuals <- sum(null_residuals^2)
null_sigma <- sqrt(null_sum_squared_residuals / (n - 1))
likelihood_null <- prod(1/(null_sigma * sqrt(2 * pi)) * exp(-null_residuals^2/(2*null_sigma^2)))

logliknull <- logLik(nullmodel) # second variant
                        
likelihood_ratio <- likelihood_null/likelihood_linear

# Nagelkerke's Pseudo-R²
pseudo_r2 <- ((1-(likelihood_ratio^(2/n)))/(1-(likelihood_null^(2/n))))
 
pseudo_r2_loglik <- (1 - exp(-2/n * (loglikmodel - logliknull)))/((1 - exp(2/n * logliknull))) #second variant

nagelkerke(testmodel)

Now, my results are:

  • n: 118
  • likelihood_linear: 28.89061
  • likelihood_null: 19.76952
  • logLik(testmodel): 3.372088 (df=3) [why 3 df here though?]
  • logLik(nullmodel): 3.068793 (df=2) [and why here 2 and not 1?]

The Pseudo-R²:

  • nagelkerke() in loop: 1.823810
  • nagelkerke() singular: -0.0960369
  • handcalculated with "normal" likelihood (pseudo_r2): -0.123545
  • handcalculated with log Likelihood (pseudo_r2_loglik): -0.09603694
  • For comparison: the adjusted R² of the lm was -0.003449

So... Still negative...Why? Can I trust these negative values? I thought that Nagelkerke's Pseudo-R² should range between 0 and 1? Looks like I have to re-calculate all of these values, also the ones which I talked about in the beginning? Which approach should I use now?

Thanks so much in advance for your help!

Dataset:

testdata <- data.frame(
  FML = c(22.29128, 22.39992, 22.50847, 22.83361, 23.37387, 23.69708, 24.01961, 24.12697, 24.44863, 24.98331, 25.09005, 25.30332, 25.51633, 26.25991, 26.3659, 26.57769, 26.78926, 27.10618, 27.52798, 27.94895, 28.8, 28.8932, 28.8932, 29.98, 30.3, 30.33, 30.4, 30.42, 30.5, 30.73, 30.9, 31.31, 31.31, 31.33, 31.39, 31.8, 32, 32.13, 32.14, 32.25, 32.47, 32.5, 32.53, 32.6, 32.62, 32.69, 32.8, 32.85, 33, 33.1, 33.22, 33.84, 33.92, 33.99, 34, 34.06, 34.4, 34.45, 34.6, 34.6, 34.61, 34.68, 34.82, 34.84, 34.87, 34.87, 34.92, 35.03, 35.07, 35.12, 35.19, 35.21, 35.35, 35.44, 35.46, 35.79, 35.89, 36.18, 36.2, 36.29, 36.33, 36.35, 36.4, 36.46, 36.5, 36.6, 36.92, 36.92, 36.96, 36.97, 37.02, 37.11, 37.21, 37.24, 37.3, 37.4, 37.44, 37.49, 37.5, 37.62, 37.78, 37.78, 37.87, 38.02, 38.05, 38.17, 38.19, 38.44, 38.6, 38.64, 38.77, 39.53, 41, 42.18, 42.5, 43.3, 43.89, 60.6),
  CSRatio21 = c(1.4774311, 1.22585, 1.3784473, 1.6238005, 1.3249576, 1.7084772, 1.5943105, 1.9541404, 1.1171914, 1.3236172, 1.2929554, 1.0386303, 1.561736, 1.3239787, 1.1792404, 1.9197001, 1.4155338, 1.543906, 1.6820963, 1.6484193, 1.4603812, 1.5447676, 1.3982134, 1.3560707, 1.684361, 1.2844309, 1.5050369, 1.246091, 1.5280047, 1.4810207, 1.8660375, 1.1933631, 1.2836483, 1.0328543, 1.7887168, 1.4201204, 1.0225164, 1.3032925, 1.2286226, 0.9198779, 1.6436054, 1.5337189, 1.537951, 1.1403919, 1.6007158, 1.337508, 1.5251319, 0.9975092, 1.7674451, 1.7825995, 1.4685762, 1.7518935, 1.5658627, 1.2009758, 1.1491171, 1.6162482, 1.8831728, 1.6755523, 1.7600169, 1.448914, 1.5329522, 1.1955977, 1.5366577, 1.4699221, 1.3493276, 1.3176134, 1.4754639, 1.4144239, 1.5749245, 1.4713806, 1.4195688, 1.371345, 1.3518266, 1.0300401, 1.4767193, 1.4706224, 1.2066236, 1.1956105, 1.3105297, 1.1632983, 1.2868495, 1.3720136, 1.1776412, 1.3619443, 1.3287843, 1.7111001, 1.2242829, 1.5002895, 1.9023916, 1.398345, 1.1438687, 1.5444296, 1.5875994, 1.5138578, 1.8444355, 1.4297542, 1.5532971, 1.0265456, 1.6225182, 1.9630474, 1.2549328, 1.2016863, 1.1794137, 1.0864437, 1.2412811, 1.3124565, 1.6559418, 1.6567076, 1.4177803, 1.9896252, 1.112104, 1.5502013, 1.384694, 1.662408, 1.326554, 1.5228765, 1.5544942, 1.8852055)
)

Also, some additional info:

  • nullmodel and linear model have the same length and the same family (see code)
  • NAs are ommited (see code)
  • my code here is not a loop anymore, my problem are the now negative values (without the loop, and one singular model + nullmodel). I did this "individual approach" to see if the Nagelkerke R² > 1 problem that I had in the loop persists - it did not, but now I have negative values... x)
  • I calculate a lm as a glm to have an adjusted R² value that I can compare the Nagelkerke values to. I will have glms with the families Gaussian, Gamma, Inverse Gaussian and the link functions identity, inverse, log, (1/mu²). The code I provided is basically at a method testing stage.
$\endgroup$
18
  • 2
    $\begingroup$ I can't remember what error checking these functions do, but do you have any missing values for FML? That would cause weird results as the data sets wouldn't be the same size for your two models. $\endgroup$ Commented Oct 13, 2023 at 16:44
  • 1
    $\begingroup$ It's hard to see how this issue can be resolved without the data. @Pauline Can you share the dataset or provide a simulated dataset which exhibits the same problem? $\endgroup$
    – dipetkov
    Commented Oct 13, 2023 at 17:39
  • 2
    $\begingroup$ All this said, I think the first thing to do is, as @dipetkov, provide a reproducible example that shows this behavior. ... I suspect your problem has to do with the loop, and maybe with ending up with models that aren't properly nested, probably due to NAs sneaking in, or maybe not having the families match on the model and null model, or something similar. $\endgroup$ Commented Oct 13, 2023 at 18:19
  • 1
    $\begingroup$ Yes, so, as @dipetkov explains, the negative Nagelkerke pseudo r-square makes sense in this case. Just as point of interest, it's similar to the McFadden value. In contrast, the Cox and Snell and Efron pseudo r-squared are slightly positive (c. 0.005). I wouldn't worry about the negative values; they make sense in these cases. Use whichever pseudo r-squared measure you think is most useful or appropriate for your purpose. FYI: You can get Efron's pseudo r-squared with library(rcompanion); efronRSquared(testmodel) or with library(DescTools); PseudoR2(testmodel, which="Efron") $\endgroup$ Commented Oct 14, 2023 at 12:40
  • 1
    $\begingroup$ Basically, the pseudo r-square is negative because the fitted model predicts the values worse than the null model. But you can describe this in terms of likelihood values if you want a more technical explanation. $\endgroup$ Commented Oct 17, 2023 at 11:25

2 Answers 2

4
$\begingroup$

There is a single highly influential data point in your dataset (an outlier?) As a result the full model y ~ x fits the data worse than the null (intercept-only) model y ~ 1. That's what the negative R-squared is telling you.

This is straightforward to see in a scatter plot of the outcome (y = CSRatio21) against the predictor (x = FML) with the two regression lines superimposed. Graphing the data should probably be the first analysis step: plot(y ~ x) is more informative than any single number such as Nagelkerke pseudo R-squared.

About the by-hand calculations not matching the rcompanion calculations. You've hard-coded the coefficient estimates and the intercept coefficient for the null model is copied incorrectly. Your calculations are also all on the original scale, not on the log scale, so there may be loss of numerical precision as well.

# hard-coded residual calculations
# residuals <- testdata$CSRatio21 - (0.0031229 * testdata$FML + 1.336591)
# null_residuals <- testdata$CSRatio21 - 1.432629

# coefficient estimates from the full and null models
coef(lm(CSRatio21 ~ FML, data = testdata))
#> (Intercept)         FML 
#>  1.33659065  0.00312292
coef(lm(CSRatio21 ~ 1, data = testdata))
#> (Intercept) 
#>    1.441449
$\endgroup$
4
  • $\begingroup$ Thank you very much! This is indeed a datapoint that is likely a real outlier / mismeasurement, I definetly will throw it out. However, I don't quite understand how this mathematically can lead to the negative value? Anyway this is helpful, thanks! $\endgroup$
    – Pauline
    Commented Oct 16, 2023 at 20:34
  • $\begingroup$ Start with the formula for Nagelkerke's R squared and compute it twice: once assuming that the null model is better (has higher log-likelihood) and then assuming that the full model is better (has higher log-likelihood). $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 20:55
  • $\begingroup$ Spoke too soon. I think with your dataset it's the case that logLik(full_model) > logLik(null_model) but the difference is small and the arithmetic works such that the numerator is positive and denominator is negative. (Plug the values in to the formula to try it out.) $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 21:12
  • $\begingroup$ So in addition we learn not to use (pseudo) R-squared to compare to bad and/or ineffective models. $\endgroup$
    – dipetkov
    Commented Oct 16, 2023 at 21:13
3
$\begingroup$

This may be more of a comment. At least for the example below, the Nagelkerke pseudo r-squared from the rcompanion package, from the DescTools package, and from your # variant 2 calculation all return the same result.

Note that in your question, your handcalculated with log Likelihood (pseudo_r2_loglik) and nagelkerke() singular return the same result.

Likely one issue is with the loop, perhaps with what you are using as the null model somewhere.

It may also be that the with the family of glm in your orginial model, that the Nagelkerke calculation behaves strangely. What was the original family ?

N = 30
CSRatio21 = rnorm(N, 0, 1)
FML = c(rep("A", N/2), rep("B", N/2))
testdata = data.frame(CSRatio21=CSRatio21, FML=FML)

testmodel <- glm(CSRatio21 ~ FML, data = testdata, family = gaussian(link="identity"))
nullmodel <- glm(CSRatio21 ~ 1, data = testdata, family = gaussian(link="identity"))

library(rcompanion)
nagelkerke(testmodel)$Pseudo.R.squared.for.model.vs.null[3]

library(DescTools)
PseudoR2(testmodel, which="Nagelkerke")

n=length(testdata[,1])
loglikmodel <- logLik(testmodel)
logliknull <- logLik(nullmodel)
pseudo_r2_loglik <- (1 - exp(-2/n * (loglikmodel - logliknull)))/((1 - exp(2/n * logliknull)))
as.numeric(pseudo_r2_loglik)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.