[R] Hosmer-Lemeshow 'goodness of fit'

[Michael Dewey]

At 19:34 29/03/2012, JimeBoHe wrote:
>I am a new user in R, so I am sorry if this is a basic question but I am kind
>of lost here and I really would appreciatte your help.

Dear Jimena

Comments in-line


>I have behavioral information of sea lions and I've done a binomial
>Generalized Linear model using Rcmdr, to understand what variables are
>affecting the place where Male's fights are being made (water or land). One
>of my possible models includes Agression Type, number of females around the
>fight and temperature in the moment.
>
>*******************
> > GLM.2 <- glm(Place ~ AgType + Females + Temp, family=binomial(cloglog),
>+   data=June)
>
> > summary(GLM.2)
>
>Call:
>glm(formula = Place ~ AgType + Females + Temp, family = binomial(cloglog),
>     data = June)
>
>Deviance Residuals:
>     Min       1Q   Median       3Q      Max
>-2.4448  -0.7248   0.4220   0.6741   1.8798
>
>Coefficients:
>             Estimate Std. Error z value Pr(>|z|)
>(Intercept)  3.43071    0.82164   4.175 2.97e-05 ***
>AgType      -1.23463    0.37166  -3.322 0.000894 ***
>Females      0.05812    0.01320   4.404 1.06e-05 ***
>Temp        -0.08614    0.02453  -3.512 0.000444 ***
>---
>Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>(Dispersion parameter for binomial family taken to be 1)
>
>     Null deviance: 175.79  on 141  degrees of freedom
>Residual deviance: 125.33  on 138  degrees of freedom
>AIC: 133.33
>
>Number of Fisher Scoring iterations: 8
>****************************************
>
>
>When I ran hosmerlem test it seems that I have a good fit to the data.
>
>************************
> > hosmerlem(Place, fitted(GLM.2))
>       X^2        Df                  P(>Chi)
>7.5123428 8.0000000 0.4824925
>*****************************
>
>However if I run only AgType, or Agtype + Temp, when running Hosmerlem test;
>it gives me an error that says: * ERROR:  'breaks' are not unique*.
>
>I've been having the same error when I try to run other models or when I use
>different link functions like logit.
>
>
>So... right now my questions are:
>
>1. What is this error means? as far as I have seen in the Hosmerlem test
>code that I am using, it seems to depend on the yhat probs, but to be honest
>I don't really understand what is this "yhat". The code I am using is:

yhat is the predicted value of y (the outcome variable). You do not 
provide us with a reproducible example so this is only a guess but I 
suspect if you look at the values of yhat for the failing models you 
will find you have many replicated values and so at least two values 
of the breakpoints are identical. But, as I say, this is at best a guess.


>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)
>}
>
>
>2. I've seen in older posts that Cessie - van Houwelingen - Copas - Hosmer
>unweighted sum of squares test for global goodness of fit test; seems to be
>a better alternative. However I am probably doing something wrong because I
>can't get it work.
>
>I installed and loaded rms package and I ran the following command:
>
> > resid(GLM.2, 'gof')
>
>As far as I understand 'gof' is another package with this description:
>Implementation of model-checking technique based on cumulative residuals.
>So, for the command I also installed and loaded 'gof' package

You seem to be so far away from success here that I think remote 
advice is not going to help.


>However I get this: *ERROR: 'arg' should be one of "deviance", "pearson",
>"working", "response", "partial"*
>
>If I use any of the arguments that it suggests, it gives me a list of
>residuals, but I doesn't show the Z or the P.

That is because it is doing what it is documented to. I am not 
familiar with the gof function nor do I see it in the index for rms 
(as you implied) so it is a bit hard to know what you are trying to do here.



>So, if someone could help me understand what I am doing wrong I really would
>appreciate it.
>
>
>
>
>--
>View this message in context: 
>http://r.789695.n4.nabble.com/Hosmer-Lemeshow-goodness-of-fit-tp3508127p4516441.html
>Sent from the R help mailing list archive at Nabble.com.

Michael Dewey
info at aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.html