Re: [R] Different goodness of fit tests leads to contradictory conclusions

[Kiyoshi Sasaki]

Dear Joris,
 
Thank you for your prompt reply! I have a binary dependent variable (whether a 
snake is pregnant or not pregnant). Independent/predictor variable is the 
snake's body size. Each observation (row) of the data represents each snake. 
One column of the data contain '0' or '1' to indicate whether a snake is 
pregnant. Another column contain body size for each snake. So, if I understand 
correctly, I can use only Hosmer-Lemeshow test. Am I correct?
 
Thank you very much!
 
Kiyoshi

--- On Wed, 7/7/10, Joris Meys <jorism...@gmail.com> wrote:


From: Joris Meys <jorism...@gmail.com>
Subject: Re: [R] Different goodness of fit tests leads to contradictory 
conclusions
To: "Kiyoshi Sasaki" <skiyoshi2...@yahoo.com>
Cc: r-help@r-project.org
Date: Wednesday, July 7, 2010, 10:08 AM


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 <skiyoshi2...@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@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.m...@ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php



      
        [[alternative HTML version deleted]]


______________________________________________
R-help@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.