I would have thought the main problem with any kind of $R^2$ measure for logistic regression is that you are dealing with a model which has a known noise value. This is unlike standard linear regression, where the noise level is usually treated as unknown. For we can write a glm probability density function as: $$f(y_i|\mu_i,\phi)=\exp\left(\frac{y_ib(\mu_i)-c(\mu_i)}{\phi}+d(y_i,\phi)\right)$$ Where $b(.),c(.),d(.;.)$ are known functions, and $\mu_i=g^{-1}(x_i^T\beta)$ for inverse link function $g^{-1}(.)$. If we define the usual glm deviance residuals as $$d_i^2=2\phi\left(log[f(y_i|\mu_i=y_i,\phi)]-log[f(y_i|\mu_i=\hat{\mu}_i,\phi)]\right)$$ $$=2\phi \left[y_ib(y_i)-y_ib(\hat{\mu}_i)-c(y_i)+c(\hat{\mu}_i)\right]$$ The we have (via likelihood ratio chi-square, $\chi^2=\frac{1}{\phi}\sum_{i=1}^{N}d_i^2$) $$E\left(\sum_{i=1}^{N}d_i^2\right)=E(\phi\chi^2)\approx (N-p)\phi$$ Where $p$ is the dimension of $\beta$. For logistic regression we have $\phi=1$, which is known. So we can use this to decide on a definite level of residual that is "acceptable" or "reasonable". This usually cannot be done for OLS regression (unless you have prior information about the noise). Namely, we expect each deviance residual to be about $1$. Too many $d_i^2>>1$ and it is likely that an important effects are missing from the model (under-fitting); too many $d_i^2<<1$ and it is likely that there are redundant or spurious effects in the model (over-fitting). (these could also mean model mispecification). Now this means that the problem for the pseudo-$R^2$ is that it fails to take into account that the level of binomial variation is predictable (provided the binomial error structure isn't being questioned). Thus even though Nagelkerke ranges from $0$ to $1$, it is still not scaled properly. Additionally, I can't see why these are called pseudo $R^2$ if they aren't equal to the usual $R^2$ when you fit a "GLM" with an identity link and normal error. For example, the equivalent cox-snell R-squared for normal error (using REML estimate of variance) is given by: $$R^2_{CS}=1-\exp\left(-\frac{N-p}{N}\frac{R^2_{OLS}}{1-R^2_{OLS}}\right)$$ Which certainly looks strange. I think the better "Goodness of Fit" measure is the sum of the deviance residuals, $\chi^2$. This is mainly because we have a target to aim for.