Relationship between binomial regression link function and goodness-of-fit tests [now with link to R code]

Question

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):

LOGIT LINK

PROBIT LINK

CAUCHIT LINK

IDENTITY LINK

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.

Answer 1

I've been able to prove both effects shown here.

Let the model matrix be $X$, an $N \times (p+1)$ matrix whose first column is the intercept column (all ones) and whose $1 \times (p+1)$ columns are $x_k^T$. The fitted value from the regression is $p_k = g(x_k^T \hat\theta)$, with the link function $g(\eta)$.

Pearson test incompatibility with identity link

First I'll consider the collision between the Pearson test and the identity link. According to Osius and Rojek (citing McCullagh and Nelder) the expected variance of the Pearson statistic is $$ \hat{\sigma}^2 = \sum_k \left[ \frac{1}{p_k(1-p_k)} - 4 \right] - c^T I^{-1} c $$ where $I$ is the information matrix at $\hat\theta$ $$ I = \sum_k\frac{g'(g^{-1}(p_k))^2}{p_k(1-p_k)} x_k x_k^T $$ and $c$ is equal to $$ c = \sum_k\frac{(1-2p_k)}{p_k(1-p_k)} g'(g^{-1}(p_k)) x_k $$ Going further with helpful matrix notation, the $N \times (p+1)$ matrix $X$ can be written $\left[ x_k^T \right]$, i.e. a column vector of row vectors. Also define an $N \times N$ matrix $B$ and an $N \times 1$ column vector $C$: \begin{align} B &= \mathrm{diag}\left(\frac{g'^2_k}{p_k(1-p_k)}\right) \\ C &= \left[\frac{(1-2 p_k) g'_k}{p_k(1-p_k)}\right] \end{align} where $g'_k = g'(g^{-1}(p_k))$. Thus \begin{align} \hat{\sigma}^2 &= C^T B^{-1} C - C^T X (X^T B X)^{-1} X^T C \\ &= C^T \left[ B^{-1} - X (X^T B X)^{-1} X^T \right] C \end{align} Defining $X'= B^{1/2} X$ and $C' = B^{-1/2} C$, this can be rewritten as $$ \hat{\sigma}^2 = C'^T \left[ I - X' (X'^T X')^{-1} X'^T \right] C' $$ This equation is of the form of quadratic form involving the hat matrix (e.g. projection matrix) based on $X'$: $$ H' = X' (X'^T X')^{-1} X'^T $$ and the vector $C'$. The matrices $X'$ and $C'$ can be written \begin{align} X' = B^{1/2} X &= \left[ \begin{array}{c} \frac{g'_k x_k^T}{\sqrt{p_k(1-p_k)}} \end{array} \right] \\ C' = B^{-1/2} C &= \left[ \frac{1-2p_k}{\sqrt{p_k(1-p_k)}}\right] \end{align} Now, when a non-identity link is used $g'_k \neq 1$ and is not constant and $p_k$ is not linearly related to $x_k$. However, using the identity link where $g'_k=1$, $C'$ becomes a vector in the column space of $X'$. Explicitly: Since in this case $p_k = x_k^T \hat\theta$, and the first element of $x_k$ is 1 for the intercept, $1 - 2 p_k = x_k^T \hat\theta'$ where $\hat\theta' = \left[1-2\hat\theta_0,-2\hat\theta_1,-2\hat\theta_2,\ldots \right]$, and thus $C' = X' \hat\theta'$. Since $C'$ is in the column space of $X'$, $H' C' = C'$. Therefore, the expected variance is $$ \hat{\sigma}^2 = C'^T \left[ I - H' \right] C = 0 $$ This proves that there is a collision between the Pearson test and the identity link.

Deviance test incompatibility with logit link

Things are very similar for the deviance statistic (as it is a part of the same power-divergence family as the Pearson statistic). Everything above still holds, adding the subscript 2 and with the replacement that $$ C_2 = \left[-2\log \left(\frac{p_k}{1-p_k} \right) g'_k\right] $$ and thus $$ C'_2 = \left[-2\log \left(\frac{p_k}{1-p_k} \right) \sqrt{p_k(1-p_k)}\right] $$ Now $C'_2$ isn't in the column space of $X'_2$ for the identity link $g'=1$ nor indeed for most links. However, with the logit link, $g'_k = p_k(1-p_k)$ and $\log \left( \frac{p_k}{1-p_k} \right) = x_k^T \hat\theta$, and thus \begin{align} X'_2 &= \left[ \sqrt{p_k(1-p_k)} x_k^T \right] \\ C'_2 &= \left[-2\sqrt{p_k(1-p_k)} x_k^T \hat\theta\right] \end{align} It immediately follows that $C'_2 = -2 X'_2 \hat\theta$, that $H'_2 C'_2 = C'_2$, and therefore $\hat{\sigma}^2_2 = 0$.

The above proves that the Pearson test does not work in the identity link and that the deviance test does not work with the logit link. The power divergence family is a continuous family of tests indexed by $\lambda$, where the Pearson statistic and deviance statistic. It seems possible that for each $\lambda$ there is a corresponding link function for which the expected variance (according to a literal use of the formula provided in the literature) is null.

General case

In the general case, for the statistic $SD_\lambda$, the vector $C$ is given by $$ C_\lambda = \frac{2}{\lambda(\lambda+1)} \left[ \left( p_k^{-\lambda } - (1-p_k)^{-\lambda} \right) g'_k \right] $$ This is derived from the summand in the definition of $SD_\lambda$, taking the difference between the $\lambda=1$ and $\lambda=0$ cases, and including a $g'$ factor. You can check that $\lambda = 1$ and $\lambda \to 0$ give the Pearson and deviance cases above. With a general link $g$, \begin{align} X'_\lambda &= \left[ \begin{array}{c} \frac{g'_k x_k^T}{\sqrt{p_k(1-p_k)}} \end{array} \right] \\ C'_\lambda &= \frac{2}{\lambda(\lambda+1)} \left[ \left( p_k^{-\lambda } - (1-p_k)^{-\lambda} \right) \sqrt{p_k(1-p_k)} \right] \end{align} In order for $C'_\lambda$ to be in the column space of $X'_\lambda$, there must be some vector $v$ such that $$ g'_k x_k^T v = \left( p_k^{-\lambda } - (1-p_k)^{-\lambda} \right) p_k (1-p_k) $$ Letting $v = \alpha[1,0,...] + \beta \theta$, $x_k^T v = v_0 + g^{-1}(p_k)$. Since the derivative of $h = g^{-1}$ is $h'(p) = 1/g'(g^{-1}(p))$, this is satisfied by a function $h$ obeying the differential equation $$ \frac{\alpha+\beta h(p)}{h'(p)} = \frac{2}{\lambda(\lambda+1)} \left( p_k^{-\lambda } - (1-p_k)^{-\lambda} \right) p_k (1-p_k) $$ Any solution to this differential equation will specify an (inverse) link function $h(p)$ incompatible with the $SD_\lambda$ test. For example, for the Pearson test ($\lambda = 1$), the entire family of such functions is $$ h_1(p) = -\frac{\alpha}{\beta}+\gamma (1-2p)^{-\beta/2} $$ which includes the identity for $\alpha=1$, $\beta=-2$, and $\gamma = -1/2$. For the deviance test, this d.e. in the limit $\lambda \to 0$ is $$ \frac{\alpha+\beta h(p)}{h'(p)} = -2 \log \left(\frac{p}{1-p} \right) p_k (1-p_k) $$ The family of solutions is $$ h_0(p) = -\frac{\alpha}{\beta}+\gamma \log \left(\frac{p}{1-p} \right)^{-\beta/2} $$ which includes the logit for $\alpha=0$, $\beta=-2$, and $\gamma=1$.

The differential equation is difficult for Mathematica to solve for arbitrary $\lambda > 0$, even with $\alpha=0$ and $\beta=1$. Trying out specific values like $\lambda=2/3$ (a choice suggested by Cressie and Reed), the solution is quite horrendous, involving the exponent of a term involving logs, powers, and two separate hypergeometric functions, so the inability to find a solution for general $\lambda$ is understandable.