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

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?

-

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.

-

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.

-
Why did you post this question and then immediately answer it? –  Macro Jun 1 '12 at 18:56
I agree with Macro. If you wanted answers fro the group you should have waited to hear what others have to say first. Anything you might get now is biased by seeing your answer. Also, if you know the answer what are you trying to prove by doing this? –  Michael Chernick Jun 1 '12 at 19:21
@Macro - Firefeather has posted four questions to this site (including this one) and answered three of them (including this one) him/herself, accepting one of his/her own answers once. A few more like this and I might start seeing a pattern! –  jbowman Jun 1 '12 at 20:42
@Macro, please see this official link: blog.stackoverflow.com/2011/07/… (it's linked on the Ask a Question page in the checkbox label at the bottom: "Answer your own question – share your knowledge, Q&A-style"). I had this question while I was doing homework (having chosen to use R instead of Minitab, though Minitab was demonstrated in class), but I didn't have enough time to type up the question and wait for a response. I figured out this workaround, and decided to share it with the community. –  Firefeather Jun 4 '12 at 14:30
@Macro, you're very welcome! I wish I could ask more questions where I don't provide the answer, and answer more questions that I didn't ask. But jbowman's right about a pattern: my contributions to the community are tending toward talking to myself. :) (At least I'm contributing somehow to the community, right?) –  Firefeather Jun 4 '12 at 14:43

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.

EDIT: last sentence removed as no longer relevant.

-
Is this a direct reply or do you want to ask a new question, with reference to this thread (cf. your last sentence)? –  chl Aug 3 '12 at 8:15