R/hosmerlem.R

In IRBLleida/UdBRpackage: Miscellaneous R Functions

#' Funció que calcula l'estadístic de calibració de Hosmer-Lemeshow
#' 
#' @param y a numerical vector with observed outcomes.
#' @param yhat a numerical vector with fitted values.
#' @param g number of groups. Degrees of freedom for the test is g-2.
#' @return A list with complete results: observed, expected, chi-square statistic, degree of freedom and p-value.
#' @export
hosmerlem <- function(y, yhat, g = 10){
  # cutyhat = cut(yhat, breaks = unique(quantile(yhat, probs=seq(0, 1, 1/g))), include.lowest = F, right = T)
  cutyhat = Hmisc::cut2(yhat, g = g)

  g = length(levels(cutyhat))

  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)

  return(list(observed = obs, expected = expect, chisq = chisq, df = g-2, p.value = p))
}