Hi to everybody. I use the following routine (i found it in the internet)
to compute the Hosmer-Lemeshow test in the framework of logistic regression.

hosmerlemeshow = function(obj, g=10) {
# first, check to see if we fed in the right kind of object
   stopifnot(family(obj)$family=="binomial" && family(obj)$link=="logit")
   y = obj$model[[1]]
   # the double bracket (above) gets the index of items within an object
   if (is.factor(y))
      y = as.numeric(y)==2
   yhat = obj$fitted.values
   cutyhat=cut(yhat,quantile(yhat,0:g/g),include.lowest=TRUE)

   obs = xtabs(cbind(1 - y, y) ~ cutyhat)
   expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
   if (any(expect < 5))
   #   warning("Some expected counts are less than 5. Use smaller number of
groups")
   chisq = sum((obs - expect)^2/expect)
   P = 1 - pchisq(chisq, g - 2)
   cat("H-L2 statistic",round(chisq,4)," and its p value",P,"\n")
   # by returning an object of class "htest", the function will perform
like the
   # built-in hypothesis tests
   return(structure(list(
      method = c(paste("Hosmer and Lemeshow goodness-of-fit test with", g,
"bins", sep=" ")),
      data.name = deparse(substitute(obj)),
      statistic = c(X2=chisq),
      parameter = c(df=g-2),
      p.value = P
   ), class='htest'))
   return(list(chisq=chisq,p.value=P))
}

However when i run it with the data set below i get the value NaN. Using
the same data set with SPSS i get a specific value.

FlightNo Temp ThermalDisast
1 66 0
2 70 1
3 69 0
4 68 0
5 67 0
6 72 0
7 73 0
8 70 0
9 57 1
10 63 1
11 70 1
12 78 0
13 67 0
14 53 1
15 67 0
16 75 0
17 70 0
18 81 0
19 76 0
20 79 0
21 75 1
22 76 0
23 58 1

Any suggestions will greatly appreciated.
With regards
Endy

	[[alternative HTML version deleted]]

