Kiyoshi Sasaki
2010-Jul-07 14:00 UTC
[R] Different goodness of fit tests leads to contradictory conclusions
I am trying to test goodness of fit for my legalistic regression using several options as shown below. Hosmer-Lemeshow test (whose function I borrowed from a previous post), Hosmer–le Cessie omnibus lack of fit test (also borrowed from a previous post), Pearson chi-square test, and deviance test. All the tests, except the deviance tests, produced p-values well above 0.05. Would anyone please teach me why deviance test produce p-values such a small value (0.001687886) that suggest lack of fit, while other tests suggest good fit? Did I do something wrong? Thank you for your time and help! Kiyoshi mod.fit <- glm(formula = no.NA$repcnd ~ no.NA$svl, family = binomial(link = logit), data = no.NA, na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = F))> # Option 1: Hosmer-Lemeshow test > mod.fit <- glm(formula = no.NA$repcnd ~ no.NA$svl, family = binomial(link = logit), data = no.NA, na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = F))> hosmerlem <- function (y, yhat, g = 10) { cutyhat <- cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), include.lowest = T) obs <- xtabs(cbind(1 - y, y) ~ cutyhat) expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat) chisq <- sum((obs - expect)^2/expect) P <- 1 - pchisq(chisq, g - 2) c("X^2" = chisq, Df = g - 2, "P(>Chi)" = P) }> hosmerlem(no.NA$repcnd, fitted(mod.fit))X^2 Df P(>Chi) 7.8320107 8.0000000 0.4500497> # Option 2 - Hosmer–le Cessie omnibus lack of fit test: > library(Design) > lrm.GOF <- lrm(formula = no.NA$repcnd ~ no.NA$svl, data = no.NA, y = T, x = T) > resid(lrm.GOF,type = "gof")Sum of squared errors Expected value|H0 SD Z P 48.3487115 48.3017399 0.1060826 0.4427829 0.6579228> # Option 3 - Pearson chi-square p-value: > pp <- sum(resid(lrm.GOF,typ="pearson")^2) > 1-pchisq(pp, mod.fit$df.resid)[1] 0.4623282> # Option 4 - Deviance (G^2) test: > 1-pchisq(mod.fit$deviance, mod.fit$df.resid)[1] 0.001687886 [[alternative HTML version deleted]]
Joris Meys
2010-Jul-07 15:08 UTC
[R] Different goodness of fit tests leads to contradictory conclusions
The first two options are GOF-tests for ungrouped data, the latter two can only be used for grouped data. According to my experience, the G^2 test is more influenced by this than the X^2 test (gives -wrongly- significant deviations more easily when used for ungrouped data). If you started from binary data, you can only use the Hosmer tests. Cheers Joris On Wed, Jul 7, 2010 at 4:00 PM, Kiyoshi Sasaki <skiyoshi2001 at yahoo.com> wrote:> I am trying to test goodness of fit for my legalistic regression using several options as shown below. ?Hosmer-Lemeshow test (whose function I borrowed from a previous post), Hosmer?le Cessie omnibus lack of fit test (also borrowed from a previous post), Pearson chi-square test, and deviance test. ?All the tests, except the deviance tests, produced p-values well above 0.05. ?Would anyone please teach me why deviance test produce p-values such a small value (0.001687886) that suggest lack of fit, while other tests suggest good fit? Did I do something wrong? > > Thank you for your time and help! > > Kiyoshi > > > mod.fit <- glm(formula = no.NA$repcnd ~? no.NA$svl, family = binomial(link = logit), data =? no.NA, na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = F)) > >> # Option 1: Hosmer-Lemeshow test >> mod.fit <- glm(formula = no.NA$repcnd ~? no.NA$svl, family = binomial(link = logit), data =? no.NA, na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = F)) > > ?> ?hosmerlem <- function (y, yhat, g = 10) > { > cutyhat <- cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), include.lowest = T) > obs <- xtabs(cbind(1 - y, y) ~ cutyhat) > expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat) > chisq <- sum((obs - expect)^2/expect) > P <- 1 - pchisq(chisq, g - 2) > c("X^2" = chisq, Df = g - 2, "P(>Chi)" = P) > } > >> hosmerlem(no.NA$repcnd, fitted(mod.fit)) > ?X^2??????? ?????????????? Df?? ??????????????????????? P(>Chi) > 7.8320107 ?????????? 8.0000000??????????? 0.4500497 > > >> # Option 2 - Hosmer?le Cessie omnibus lack of fit test: >> library(Design) >> lrm.GOF <- lrm(formula = no.NA$repcnd ~? no.NA$svl, data =? no.NA, y = T, x = T) >> resid(lrm.GOF,type = "gof") > Sum of squared errors???? Expected value|H0??????? SD???????????????????? Z???????????????????? P > ????????? ????? ?48.3487115?????????? ????????????? ?48.3017399??????????? ???????????? ?0.1060826???? 0.4427829???? 0.6579228 > >> # Option 3 - Pearson chi-square p-value: >> pp <- sum(resid(lrm.GOF,typ="pearson")^2) >> 1-pchisq(pp, mod.fit$df.resid) > [1] 0.4623282 > > >> # Option 4 - Deviance (G^2) test: >> 1-pchisq(mod.fit$deviance, mod.fit$df.resid) > [1] 0.001687886 > > > > ? ? ? ?[[alternative HTML version deleted]] > > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- Joris Meys Statistical consultant Ghent University Faculty of Bioscience Engineering Department of Applied mathematics, biometrics and process control tel : +32 9 264 59 87 Joris.Meys at Ugent.be ------------------------------- Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php