I encountered a strange phenomenon when calculating pseudo R2 for logistic models when using aggregated files: the results are simply too good to be true. An example (but as far as I can see, every aggregated file offers similar problems):
library(pscl) cuse <- read.table("http://data.princeton.edu/wws509/datasets/cuse.dat", header=TRUE) head(cuse) cuse.fit <- glm( cbind(using, notUsing) ~ age + education + wantsMore, family = binomial, data=cuse) summary(cuse.fit) pR2(cuse.fit)
The results are:
> summary(cuse.fit) Call: glm(formula = cbind(using, notUsing) ~ age + education + wantsMore, family = binomial, data = cuse) Deviance Residuals: Min 1Q Median 3Q Max -2.5148 -0.9376 0.2408 0.9822 1.7333 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.8082 0.1590 -5.083 3.71e-07 *** age25-29 0.3894 0.1759 2.214 0.02681 * age30-39 0.9086 0.1646 5.519 3.40e-08 *** age40-49 1.1892 0.2144 5.546 2.92e-08 *** educationlow -0.3250 0.1240 -2.620 0.00879 ** wantsMoreyes -0.8330 0.1175 -7.091 1.33e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 165.772 on 15 degrees of freedom Residual deviance: 29.917 on 10 degrees of freedom AIC: 113.43 Number of Fisher Scoring iterations: 4 > pR2(cuse.fit) llh llhNull G2 McFadden r2ML -50.7125647 -118.6401419 135.8551544 0.5725514 0.9997947 r2CU 0.9997950
The last three outcomes from pscl function pR2 present McFadden's pseudo r-squared, Maximum likelihood pseudo r-squared (Cox & Snell) and Cragg and Uhler's or Nagelkerke's pseudo r-squared. The calculation seems to be flawless, but the outcomes close to 1 seem to good to be true.
Using weight instead of cbind:
cuse2 = rbind(cuse,cuse) cuse2$using.contraceptive=1 cuse2$using.contraceptive[1:nrow(cuse)]=0 cuse2$freq = cuse2$notUsing cuse2$freq[1:nrow(cuse)] = cuse2$using[1:nrow(cuse)] cuse.fit2 = glm(using.contraceptive ~ age + education + wantsMore, weight=freq, family = binomial, data = cuse2) summary(cuse.fit2) round(pR2(cuse.fit2),5)
produces different logistic regression coefficients, and slightly different pseudo R2's for r2ml and r2CU and a large difference for McFadden R2:
> round(pR2(cuse.fit2),5) llh llhNull G2 McFadden r2ML -933.91920 -1001.84677 135.85515 0.06780 0.98567 r2CU 0.98567
Full expansion results in very different estimates from pR2:
cuse3 = rbind(cuse[rep(1:nrow(cuse), cuse[["notUsing"]]), ], cuse[rep(1:nrow(cuse), cuse[["using"]]), ]) cuse3$using.contraceptive=1 cuse3$using.contraceptive[1:sum(cuse$notUsing)]=0 summary(cuse3) cuse.fit3 = glm(using.contraceptive ~ age + education + wantsMore, family = binomial, data = cuse3) summary(cuse.fit3) round(pR2(cuse.fit3),5) > round(pR2(cuse.fit3),5) llh llhNull G2 McFadden r2ML -933.91920 -1001.84677 135.85515 0.06780 0.08106 r2CU 0.11376
This indicates a logistic model which explains very little, which is a little bit more believable than the near perfect results from the aggregated files. Is there a more correct, and preferably more consistent, way to calculate the pseudo R2's?