burke.R

############################################################
# Hosmer-Lemeshow test
############################################################

hosmerlem = function(observed, predicted, g = 10) {
  cutyhat = cut(predicted,
                breaks = quantile(predicted,
                                  probs=seq(0, 1, 1/g)), 
                include.lowest=TRUE)
  
  obs = xtabs(cbind(1 - observed, observed) ~ cutyhat)
  expect = xtabs(cbind(1 - predicted, predicted) ~ cutyhat)
    
  chisq = sum((obs - expect)^2/expect)
  P = 1 - pchisq(chisq, g - 2)
  
  return(list(cHat = chisq, p.value = P))
}

############################################################
# Pigeon-Heyse test
############################################################

pigeonheyse <- function(observed, predicted, g = 10) {
  cutyhat <- as.character(cut(predicted,
                              breaks = quantile(predicted,
                                                probs=seq(0, 1, 1/g)), 
                              include.lowest=TRUE))
  data$cutyhat <- cutyhat
  unq <- unique(cutyhat)
  for(i in 1:g) {
    data[data$cutyhat == unq[i],4] <- i
  }
  
  Phi = NULL
  for(i in 1:g) {
    Phi <- rbind(Phi, (sum(data[data$observed==1 & data$cutyhat==i, 2] *
                             (1 - data[data$observed==1 & data$cutyhat==i, 2]))) /
                   (length(data[data$observed==1 & data$cutyhat==i, 2]) *
                      mean(data[data$observed==1 & data$cutyhat==i, 2]) *
                      (1 - mean(data[data$observed==1 & data$cutyhat==i, 2])))
    )
  }
  
  Phi <- cbind(Phi, Phi)
  cutyhat = cut(predicted,
                breaks = quantile(predicted,
                                  probs=seq(0, 1, 1/g)), 
                include.lowest=TRUE)
  
  obs = xtabs(cbind(1 - observed, observed) ~ cutyhat)
  expect = xtabs(cbind(1 - predicted, predicted) ~ cutyhat)
  
  chisq = sum((obs - expect)^2/(Phi*expect))
  P = 1 - pchisq(chisq, g - 2)
  
  return(list(J2 = chisq, p.value = P))
}

############################################################
# Basic data manipulation
############################################################

year <- c(as.character(nflprobdata$YEAR.WEEK))
year <- strsplit(year, "-")
year <- as.numeric(lapply(year, "[[",1))

data <- data.frame(year = year,
                   predicted = nflprobdata$Home.Team.Probability,
                   observed = nflprobdata$Did.Home.Team.Win.)
data <- data[order(data[ ,2]), ]
rm(year)

############################################################
# Find our C-hat and J^2!
############################################################

# Overall
with(data, 
     hosmerlem(observed,
               predicted, g = 10))

with(data, 
     pigeonheyse(observed,
               predicted, g = 10))

years <- c(unique(data$year))
for(i in 1:length(years)) {
  print(years[i])
  with(data[data$year == years[i], ],
       print(hosmerlem(observed,
                       predicted, g = 10)))
  with(data[data$year == years[i], ],
       print(pigeonheyse(observed,
                         predicted, g = 10)))
}

nflprob.csv