component average functions.R

# functions to compute and graph component estimates of a batting average
# Jim Albert
# JQAS paper "Improved Component Predictions of Batting and Pitching Measures"

fit_comp_half <- function(d){
  # input d - a list with components playerID, AB, H, HR, SO
  # output - a list with components
  #     S - a data frame with all of the component and shrinkage estimates
  #     component -- values of K and eta for all the component fits
  #     shrinkage -- values of K and eta for the shrinkage fit
  require(dplyr)
  require(LearnBayes)

  ## ------------------------------------------------------------------------
  fit.model <- function(data){
    mode <- laplace(betabinexch, c(1, 1), 
                    cbind(data$y, data$n))$mode
    eta <- exp(mode[1]) / (1 + exp(mode[1]))
    K <- exp(mode[2])
    list(eta=eta, K=K, 
         d=data.frame(data, est=(data$y + K * eta) / (data$n + K)))
  }
  ## ------------------------------------------------------------------------
  S.SO <- fit.model(data.frame(playerID=d$playerID, 
                               y=d$SO, n=d$AB))
  S.HR <- fit.model(data.frame(playerID=d$playerID,
                               y=d$HR, n=d$AB - d$SO))
  S.H <- fit.model(data.frame(playerID=d$playerID,
                              y=d$H - d$HR, 
                              n=d$AB - d$HR - d$SO))
  S <- merge(S.SO$d, S.HR$d, by="playerID")
  S <- merge(S, S.H$d, by="playerID")
  names(S) <- c("playerID", "SO", "AB", "SO.Rate",
                "HR", "AB.SO", "HR.Rate",
                "H.HR", "AB.SO.HR", "H.Rate")
  
  component.fit <- data.frame(eta=c(S.SO$eta, S.HR$eta, S.H$eta),
                              K=c(S.SO$K, S.HR$K, S.H$K))
  row.names(component.fit) <- c("SO", "HR", "H")
  
  ## ------------------------------------------------------------------------
  S$Est <- with(S, 
                (1 - SO.Rate) * (HR.Rate + (1 - HR.Rate) * H.Rate))
  
  ## ------------------------------------------------------------------------
  S2 <- fit.model(data.frame(playerID=d$playerID, 
                             y=d$H, n=d$AB))
  shrinkage.fit <- c(eta=S2$eta, K=S2$K)
  
  S <- merge(S, S2$d, by="playerID")
  names(S)[c(11:14)] <- c("Comp.Est", "H", "AB1", "Shrinkage.Est")
  list(S=S, component=component.fit, shrinkage=shrinkage.fit)
}

plot_avg_results2 <- function(fitwork){
  # function to plot the results from fit_comp_half
  # will construct a graph of the observed One Minus SO.Rate and HIP Rate 
  # also constructs a graph of the estimates of One Minus SO.Rate and HIP Rate
  # returns the variables of the two plots (can use grid.arrange to put them together)
  ## ------------------------------------------------------------------------
  require(ggplot2)
  require(dplyr)
  S <- fitwork$S
  myf <- function(x, p) p / x
  AVG1 <- data.frame(P=.25, x=seq(.6, 1, .01), y=myf(seq(.6, 1, .01), .25))
  AVG2 <- data.frame(P=.3, x=seq(.6, 1, .01), y=myf(seq(.6, 1, .01), .3))
  AVG3 <- data.frame(P=.2, x=seq(.6, 1, .01), y=myf(seq(.6, 1, .01), .2))
  AVG4 <- data.frame(P=.35, x=seq(.6, 1, .01), y=myf(seq(.6, 1, .01), .35))
  AVG5 <- data.frame(P=.4, x=seq(.6, 1, .01), y=myf(seq(.6, 1, .01), .4))
  AVG <- rbind(AVG1, AVG2, AVG3, AVG4, AVG5)
  AVG <- mutate(AVG, BA=factor(P))
  d1 <- data.frame(Type="Observed BA",
                   One.Minus.SO.Rate=1 - S$SO / S$AB, 
                   H.Rate=S$H / S$AB.SO,
                   AB=S$AB)
  d2 <- data.frame(Type="Predicted BA",
                   One.Minus.SO.Rate=1 - S$SO.Rate,
                   H.Rate=S$HR.Rate + (1 - S$HR.Rate) * S$H.Rate,
                   AB=S$AB)
  Xlow <- range(c(d1$One.Minus.SO.Rate, d2$One.Minus.SO.Rate))
  Ylow <- range(c(d1$H.Rate, d2$H.Rate))
  p1 <- ggplot(d1, aes(One.Minus.SO.Rate, H.Rate)) + 
    geom_point() +
    geom_line(data=AVG, aes(x, y, color=BA)) +
    ylim(Ylow[1], Ylow[2]) + xlim(Xlow[1], Xlow[2]) +
    ggtitle("Observed BA") +
    labs(y = "Hits in Non-SO-AB Rate") 
  print(p1)
  p2 <- ggplot(d2, aes(One.Minus.SO.Rate, H.Rate)) + 
    geom_point(color="red") +
    geom_line(data=AVG, aes(x, y, color=BA)) +
    ylim(Ylow[1], Ylow[2]) + xlim(Xlow[1], Xlow[2]) +
    ggtitle("Predicted BA") +
    labs(y = "Hits in Non-SO-AB Rate") 
  print(p2)
  list(plot1=p1, plot2=p2)
}