How can I compute Pearson's $\chi^2$ test statistic for lack of fit on a logistic regression model in R?

Question

The likelihood ratio (a.k.a. deviance) $G^2$ statistic and lack-of-fit (or goodness-of-fit) test is fairly straightforward to obtain for a logistic regression model (fit using the glm(..., family = binomial) function) in R. However, it can be easy to have some cell counts end up low enough that the test is unreliable. One way to verify the reliability of the likelihood ratio test for lack of fit is to compare its test statistic and P-value to those of Pearson's chi square (or $\chi^2$) lack-of-fit test.

Neither the glm object nor its summary() method report the test statistic for Pearson's chi square test for lack of fit. In my search, the only thing I came up with is the chisq.test() function (in the stats package): its documentation says "chisq.test performs chi-squared contingency table tests and goodness-of-fit tests." However, the documentation is sparse on how to perform such tests:

If x is a matrix with one row or column, or if x is a vector and y is not given, then a goodness-of-fit test is performed (x is treated as a one-dimensional contingency table). The entries of x must be non-negative integers. In this case, the hypothesis tested is whether the population probabilities equal those in p, or are all equal if p is not given.

I'd imagine that you could use the y component of the glm object for the x argument of chisq.test. However, you can't use the fitted.values component of the glm object for the p argument of chisq.test, because you'll get an error: "probabilities must sum to 1."

How can I (in R) at least calculate the Pearson $\chi^2$ test statistic for lack of fit without having to run through the steps manually?

Answer 1

The sum of the squared Pearson residuals is exactly equal to the Pearson $\chi^2$ test statistic for lack of fit. So if your fitted model (i.e., the glm object) is called logistic.fit, the following code would return the test statistic:

sum(residuals(logistic.fit, type = "pearson")^2)

See the documentation on residuals.glm for more information, including what other residuals are available. For example, the code

sum(residuals(logistic.fit, type = "deviance")^2)

will get you the $G^2$ test statistic, just the same as deviance(logistic.fit) provides.

Answer 2

The Pearson statistic has a degenerate distribution so is not recommended in general for logistic model goodness-of-fit. I prefer structured tests (linearity, additivity). If you want an omnibus test see the single degree of freedom le Cessie - van Houwelingen - Copas - Hosmer unweighted sum of squares test as implemented in the R rms package residuals.lrm function.

Answer 3

Thanks, I didn't realize it was as simple as: sum(residuals(f1, type="pearson")^2) However please note that Pearsons residual varies depending on whether it is calculated by covariate group or by individual. A simple example:

m1 is a matrix (this one is the head of a larger matrix):

m1[1:4,1:8]

    x1 x2 x3 obs    pi   lev  yhat y
obs  1  1 44   5 0.359 0.131 1.795 2
obs  0  1 43  27 0.176 0.053 4.752 4
obs  0  1 53  15 0.219 0.062 3.285 1
obs  0  1 33  22 0.140 0.069 3.080 3

Where x1-3 are predictors, obs is no. observations in each group, pi is probability of group membership (predicted from regression equation), lev is leverage, the diagonal of the hat matrix, yhat the predicted no. (of y=1) in the group and y the actual no.

This will give you Pearson's by group. Note how it's different if y==0: '$ fun1 <- function(j){ if (m1[j,"y"] ==0){ # y=0 for this covariate pattern Pr1 <- sqrt( m1[i,"pi"] / (1-m1[i,"pi"])) Pr2 <- -sqrt (m1[i,"obs"]) res <- round( Pr1 * Pr2, 3) return(res) } else { Pr1 <- m1[j,"y"] - m1[j,"yhat"] Pr2 <- sqrt( m1[j,"yhat"] * ( 1-(m1[j,"pi"]) ) ) res <- round( Pr1/Pr2, 3) return(res) } } $'

Thus

nr <-nrow(m1)
nr1 <- seq(1,nr)
Pr <- sapply(1:nrow[m1], FUN=fun1)
PrSj <- sum(Pr^2)

If there are large numbers of subjects with y=0 covariate patterns, then Pearons residual will be much larger when calculated using the 'by group' rather than the 'by indiviual' method.

See e.g. Hosmer & Lemeshow "Applied Logistic Regression", Wiley, 200.

Answer 4

You can also use c_hat(mod) that will give the same output as sum(residuals(mod, type = "pearson")^2).