# Truncate predictions at .0001 so that a predicted probability of 0 doesn't give log(0)
pred_loglik <- function(p) sum(dbinom(oos$y_2022, oos$n_2022, prob = pmin(pmax(p*scale_2022, .0001), .9999), log=T))
strengths <- c(50, 200, 500, 1000)
pull_loglik <- sapply(strengths, function(k) pred_loglik((oos$y + k*state_rate)/(oos$n + k)))
best_k <- strengths[which.max(pull_loglik)] # The best-predicting prior strength (we'll map it below)
tibble(estimator = c("MLE", "Uniform prior", "Statewide rate only", paste("Pull to statewide rate, prior strength", strengths)),
loglik = c(pred_loglik(oos$mle), pred_loglik((oos$y + 1)/(oos$n + 2)), pred_loglik(rep(state_rate, nrow(oos))), pull_loglik))