Pseudo R2 Values For Negative Binomial Regression Model in R Yields Inconsistent, Sometimes Negative Values

I am getting some odd results and I feel like I have been banging my head against the wall for a while, so here I am. I am running a negative binomial regression analysis using MASS:glm.nb:

summary(m1 <- glm.nb(NumberOfAppUses ~
MainCountryIncome +
offset(log(SpanDays/30)),
data=df))


The results are quite promising and consistent with visual inspection of the data, which demonstrates a clear relationship.

    Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                           3.17752    0.04399  72.239  < 2e-16 ***
MainCountryIncomeLower middle income -0.31888    0.04612  -6.913 4.73e-12 ***
MainCountryIncomeUpper middle income -0.43847    0.04588  -9.557  < 2e-16 ***
MainCountryIncomeHigh income         -0.67500    0.04608 -14.647  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.9664) family taken to be 1)

Null deviance: 21867  on 18311  degrees of freedom
Residual deviance: 21427  on 18308  degrees of freedom
(12861 observations deleted due to missingness)
AIC: 170052

Number of Fisher Scoring iterations: 1

Theta:  0.96640
Std. Err.:  0.00936

2 x log-likelihood:  -170041.80900


I am trying to report out pseudo-R2 values for this analysis and I am getting the following results using the packages "modEvA" and "pscl":

> modEvA::RsqGLM(model=m1)
NOTE: Tjur R-squared applies only to binomial GLMs
$CoxSnell [1] -0.2660243$Nagelkerke
[1] -0.2660555

$McFadden [1] -0.02606445$Tjur
[1] NA

$sqPearson [1] 0.09906378  And: > write.table(pscl::pR2(m1),sep="\t") "x" "llh" -85020.9044052324 "llhNull" -123236.91762709 "G2" 76432.0264437156 "McFadden" 0.310101988573731 "r2ML" 0.984607524589421 "r2CU" 0.98460892997297  When I do a linear regression on the frequency data against my outcome, I get sane results. While the data itself is has a nice normal distribution after log-normalization, the data is count/rate data and therefore we went for a Poisson model. The data is overdispersed, hence the usage of the negative binomial approach. > (summary(fit<-glm(LogNatFrequency_UsesPer30Days ~ + MainCountryIncome, data = df))) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.79566 0.04660 -17.073 < 2e-16 *** MainCountryIncomeLower middle income -0.30098 0.04886 -6.160 7.41e-10 *** MainCountryIncomeUpper middle income -0.45663 0.04860 -9.396 < 2e-16 *** MainCountryIncomeHigh income -0.69500 0.04881 -14.237 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 1.211936) Null deviance: 22764 on 18311 degrees of freedom Residual deviance: 22188 on 18308 degrees of freedom (12861 observations deleted due to missingness) AIC: 55493 Number of Fisher Scoring iterations: 2 > modEvA::RsqGLM(model=fit) NOTE: Tjur R-squared applies only to binomial GLMs$CoxSnell
[1] 0.02530993

$Nagelkerke [1] 0.02656091$McFadden
[1] 0.008389995

$Tjur [1] NA$sqPearson
[1] 0.02530993

> write.table(pscl::pR2(fit),sep="\t")
"x"
"llh"   -27741.5566772792
"llhNull"   -28004.8880168262
"G2"    526.662679093963
"r2ML"  0.0283508748160685
"r2CU"  0.0297475710989573


I'm not a statistician and feel that I am well underwater already. Any thoughts?

Particularly:

1. Why do you think negative pseudo R2 values are being calculated?

2. As I understand it, a negative value is meaningless, so how should I interpret this?

3. Is this a limitation of the packages I am using to calculate or is there something I do not understand about negative binomial modelling itself such that I am using the wrong approach to goodness-of-fit evaluation?

• So what exactly is your question? – mdewey Dec 11 '16 at 21:56
• 1. Why do you think negative pseudo R2 values are being calculated? 2. As I understand it, a negative value is meaningless, so how should I interpret this? 3. Is this a limitation of the packages I am using to calculate or is there something I do not understand about negative binomial modelling itself such that I am using the wrong approach to goodness-of-fit evaluation? - Added to question also – heights1976 Dec 11 '16 at 22:01
• It's a limitation of pseudo-R2, I think. There is more than one, and none are entirely satisfactory. – Jeremy Miles Dec 12 '16 at 0:30
• In this article, Martin & Hall note that the pseudo-R2 calculations can lead to slightly negative values. None of the pseudo-R2 calculations are "actual" R2. They propose R2 equations, and propose an adjustment to the formula based on the parameters that would eliminate negative values, calling it "Adjusted-R2", much like in regular regression, the Adjusted-R2 compensates for the number of predictor variables. tandfonline.com/doi/abs/10.1080/… – jeramy townsley Nov 2 '17 at 14:48
• btw, you can find the formula for the McFaddens, etc, and should be able to figure out how they are being calculated, and thus find the math behind the negative values. For example, this Oxford Handbook has those formula. Back to your original question #2, how should you interpret these pseudo R2 values-- presumably the same way as a regular R2, and since these are approaching zero, your model seems to be explaining very little of the variance. global.oup.com/academic/product/… – jeramy townsley Nov 2 '17 at 14:53