taronesZStat <- function(model){

  #Extracting model residuals
  res <- residuals(model)

  #Number of residuals
  n <- length(res)

  #Calculating Tarone's Z-statistic
  p_hat <- sum(res)/n

  S <- sum((res - n * p_hat)^2 / (p_hat * (1 - p_hat)))

  Z_score <- (S - sum(n)) / sqrt(2 * sum(n * (n - 1)))

  return(Z_score)
}