Some background: A number of papers in the literature (various ones by Hosmer and Lemeshow; Copas; le Cessie and van Houwelingen; Cressie and Read; Osius and Rojek; J. R. Dale) discuss a family of power-divergence statistics $SD_\lambda$ for binomial regression whose distributions are asymptotically normal; see in particular this paper by Osius and Rojek. For binomial regression (rather than the more general multinomial case), with the small restriction that $\hat \pi$ is never equal to 0 or 1, the statistic $SD_\lambda$ is defined by $$ SD_\lambda = \sum \frac{2}{\lambda+1} \left[ \frac{Y}{\lambda} \left( \frac{Y}{\hat\pi} \right)^{\lambda} + \frac{1-Y}{\lambda} \left( \frac{1-Y}{1-\hat\pi} \right)^\lambda - \frac{1}{\lambda} \right] $$ This family includes the deviance $$ D = G^2 = \lim_{\lambda \to 0} SD_\lambda = \sum \left( Y \log \frac{Y}{\hat \pi} + (1-Y) \log \frac{1-Y}{1-\hat \pi} \right),$$ and the Pearson "chi-square" statistic $$ X^2 = SD_1 = \sum \frac{(Y-\hat\pi)^2}{\hat\pi(1-\hat\pi)}, $$ A similar quantity is the Copas unweighted sum of squares $$ S = \sum (Y-\pi)^2, $$ which is not a power divergence statistic but shares many of the same properties, including that it asymptotically follows a normal distribution. Since each asymptotically follows a normal distribution, if the regression model is valid, each can be used to produce a goodness-of-fit test from a standard normal distribution. The last one, notably, is used in R's rms package, under the name "le Cessie–van Houwelingen–Copas–Hosmer test."
I am working on a binomial regression tool in C++ that allows several possible link functions: the standard logit as well as probit, cauchit, and the identity (not advisable because probabilities could become negative or greater than one, but needed for historical reasons). I've also implemented the tests for the three power divergence statistics mentioned above. However, I've discovered something strange: while the test based on $S$ is stable, the test based on the deviance is unstable when using a logit link while the test based on $X^2$ is unstable for the identity link. (All work for the probit and cauchit links). By "unstable," I mean that in toy studies their distributions are not standard-normal. I believe this is largely because of numerical problems. Here are the distributions of the tests (from toy studies generating 1000 trials with 500 observations each; there are correlations in the statistics because they are applied to the same data for each toy):
As you can see, the tests fail to be standard normal in the deviance-logit and Pearson-identity combinations. In each case, I find that the variance of the statistic (before rescaling) is extraordinarily tiny. For example, in the Pearson/identity case, the expected value of the statistic is exactly the number of events, 500, while the observed value agrees with this to a large number of decimal places, and the variance might be below machine precision because the value calculated is often ~ 1e-15 and even negative.
Is this behavior expected? Is it perhaps the case that these tests sometimes have extraordinarily small variances, depending on the link function used?
SECOND EDIT
To keep the question concise I've removed the C++ code that I use to calculate these statistics from the tex and uploaded it to to GitHub: gof.cpp. There are some references to outside classes, but these are mostly self-explanatory. Table
is a container of structs, whose values include pi
, the regressed probability; dpi
, the gradient of the regression function at the fit point; correct
, a boolean corresponding to Y above; and W
, a sample weight which is identically equal to 1 in all these tests. selfCalibrated
is true in all these tests (it would be false in a validation setting).
I've also (mostly) successfully ported the code to R. This code is completely self-contained and only uses the base R packages, and includes code to reproduce the toy studies I've shown here. You can find it here: GitHub: gof.R. This code mostly reproduces my C++ results, but seems to have an additional problem with the Pearson statistic: wherever the C++ code finds a good standard normal distribution, the R code shows a positively skewed distribution. But if I apply the C++ and R implementations of the tests to the same toy data, I get nearly the same results -- which is good but doesn't explain the difference in the distributions. This is odd, but probably just a bug from my port.
In both these implementations, there is only one obvious place where the link function shows up (hidden in the C++ case). For the regression model $y = g(\beta_0 + \beta_1 x)$, where $g$ is the link function, the quantities
$$
\hat \pi(x) = g( \hat \beta_0 + \hat \beta_1 x)
$$
and
$$
\frac{d \hat \pi(x)}{d \beta} = g'(\hat \pi(x)) \left[ 1, x \right],
$$
represented by pi
and dpi
respectively, are the only obvious place where the link function shows up in these calculations.