HL.gof <-
function (y, yhat, g=10) {
  ret     <- list(method="Hosmer and Lemeshow Goodness of Fit test")
  endpts  <- quantile(yhat, probs=seq(0, 1, 1/g))
  bins    <- cut(yhat, breaks=endpts, include.lowest=TRUE)
  ob      <- xtabs(cbind(1 - y, y) ~ bins)
  ex      <- xtabs(cbind(1 - yhat, yhat) ~ bins)
  H       <- sum((ob - ex)^2/ex)
  p       <- 1 - pchisq(H, g-2)

  ret$data.name <- deparse(substitute(y))

  ret$statistic <- H
    names(ret$statistic) <- "H"
  ret$parameter <- g-2
    names(ret$parameter) <- "df"
  ret$p.value <- p
  ret$alternative <- "Model does not fit the data"

  if( any(table(bins)<5) ) {
    warning("The HS statistic assumes no bin has fewer than five elements. \n\tPlease reduce the number of groups (g).")
  }

  if( H<0 ) {
    warning("This procedure only works on y (and yhat) between 0 and 1.")
  }

  class(ret) <- "htest"
return(ret)
}