ewenharrison / finalfit
1
#' Hosmer-Lemeshow goodness of fit test
2
#'
3
#' Internal, not usually called directly
4
#'
5
#' @param y Observed y, usually of the form \code{fit$y}.
6
#' @param yhat Predicted y_hat, usually for the form \code{fit$fitted}
7
#' @param g Number of bins to calculate quantiles.
8
#' @param digits Number of decimal places of form \code{c(2,3)}, where \code{digits[1]} is
9
#'   for chi-sq estimate and \code{digits[2]} is for p-value.
10
#'
11
#' @return Character string of chi-sq result, df, and p-value. Significant
12
#'   p-value suggests poor fit.
13
#' @export
14
#' @importFrom stats pchisq xtabs
15
#'
16
#' @author Adapted from Peter Solymos.
17
#' @source https://github.com/psolymos/ResourceSelection/blob/master/R/hoslem.test.R
18
#'
19
#' @examples
20
#' fit = glm(mort_5yr~age.factor+extent.factor, data=colon_s, family="binomial")
21
#' metrics_hoslem(fit$y, fit$fitted)
22
metrics_hoslem <- function(y, yhat, g=10, digits = c(2,3)) {
23 1
  qq <- unique(quantile(yhat, probs=seq(0, 1, 1/g)))
24 1
  yhat_cut <- cut(yhat, breaks = qq, include.lowest = TRUE)
25 1
  observed <- xtabs(cbind("y0" = 1-y, "y1" = y) ~ yhat_cut)
26 1
  expected <- xtabs(cbind("yhat0" = 1-yhat, "yhat1" = yhat) ~ yhat_cut)
27 1
  chisq <- sum((observed - expected)^2 / expected)
28 1
  p = 1 - pchisq(chisq, g-2)
29 1
  par <- g-2
30 1
  out = paste0("Chi-sq(", par, ") ", round_tidy(chisq, digits = digits[1]), " (p", p_tidy(p, digits = digits[2]), ")")
31 1
  return(out)
32
}

Read our documentation on viewing source code .

Loading