High-level Findings

• There are ~2300 observations with 16 variables included
• Five of the variables have a number of missing values
• Two of the variables have outliers
• Multicollinearity present in seven variables
• Untransformed (raw) data is not highly correlated with response

# read in data.  
mb_train <- paste0("https://raw.githubusercontent.com/jbrnbrg/",
                   "msda-datasets/master/moneyball-training-data.csv")
mb_eval <- paste0("https://raw.githubusercontent.com/jbrnbrg/",
                   "msda-datasets/master/moneyball-evaluation-data.csv")

bsb_train_raw <- read_csv(mb_train)

## 
## -- Column specification --------------------------------------------------------
## cols(
##   INDEX = col_double(),
##   TARGET_WINS = col_double(),
##   TEAM_BATTING_H = col_double(),
##   TEAM_BATTING_2B = col_double(),
##   TEAM_BATTING_3B = col_double(),
##   TEAM_BATTING_HR = col_double(),
##   TEAM_BATTING_BB = col_double(),
##   TEAM_BATTING_SO = col_double(),
##   TEAM_BASERUN_SB = col_double(),
##   TEAM_BASERUN_CS = col_double(),
##   TEAM_BATTING_HBP = col_double(),
##   TEAM_PITCHING_H = col_double(),
##   TEAM_PITCHING_HR = col_double(),
##   TEAM_PITCHING_BB = col_double(),
##   TEAM_PITCHING_SO = col_double(),
##   TEAM_FIELDING_E = col_double(),
##   TEAM_FIELDING_DP = col_double()
## )

mini_tbl<-psych::describe(bsb_train_raw[,-1], IQR=T)[,c(1:5,8:10,12,14)]

bsb_train <- data.table::setnames(bsb_train_raw, 
                                  tolower(names(bsb_train_raw[1:17])))
bsb_train[is.na(bsb_train)] <- 0  # used for simple plots - not used for model fitting
                                  # see data preperation section for more details.   
col_name_updater <- function(x, text_to_updt, replc_text){
  colnames(x) = gsub(text_to_updt, replc_text, colnames(x))
  return(x)
}
# shorten and clean colnames 
c1 <- col_name_updater(bsb_train, "batting", "bat")
c2 <- col_name_updater(c1, "pitching", "pitch")
c3 <- col_name_updater(c2, "fielding", "field")
c4 <- col_name_updater(c3, "baserun", "baser")
bsb_data <- col_name_updater(c4, "team", "t")
rm(c1,c2,c3, c4)

c1 <- col_name_updater(bsb_train_raw, "batting", "bat")
c2 <- col_name_updater(c1, "pitching", "pitch")
c3 <- col_name_updater(c2, "fielding", "field")
c4 <- col_name_updater(c3, "baserun", "baser")
bsb_train_raw <- col_name_updater(c4, "team", "t")
rm(c1,c2,c3, c4)

Findings Detail

The below investigation was completed prior to model-fitting in order to get a sense for the density and quality of the training data set and to preview what types of actions may need to be taken in the model-fitting steps.

rm(mini_tbl)

p1 <- bsb_data %>%
  gather(-index, key = "vars", value = "vals") %>% 
  select_("vars", "vals") %>% 
  ggplot() + 
  geom_boxplot(aes(x = vars, y=vals)) +
  coord_flip() + xlab("")

p2 <- bsb_data %>%
  gather(-index, key = "vars", value = "vals") %>% 
  select_("vars", "vals") %>% 
  ggplot() + 
  geom_boxplot(aes(x = vars, y=vals)) +
  scale_y_log10() + 
  coord_flip() + ylab("log(vals)") + xlab("")

suppressMessages(gridExtra::grid.arrange(p1, p2, nrow=1))

# check correlations between prediction vars:  
# GGally::ggcorr(bsb_data[,c(6,8,12,13,15:17)], 
#                geom = "circle", nbreaks = 6, 
#                label=T, label_size = 3)

# see if there's cor btween target wins and others... 
sm_tab <- cor(bsb_data[,-1]) %>% as.data.frame() %>% 
  tibble::rownames_to_column(var = "var_name") %>% 
  mutate(cor_val = round(target_wins , 4), 
         abs_cor_val = abs(cor_val)) %>% 
  select_("var_name", "abs_cor_val", "cor_val") %>% 
  arrange(desc(abs_cor_val)) %>% 
  na.omit()

bsb_data[,-1] %>%
  gather(-target_wins, key = "vars", value = "vals") %>%
  ggplot(aes(x = vals, y = target_wins)) +
    geom_point(alpha=.1) +
    facet_wrap(~ vars, scales = "free", ncol = 3)

Model 1 Data:

Motivation: Based on the most complete set of records with minimal imputation:

• For TEAM_BATTING_SO and TEAM_PITCHING_SO, NA-values replaced with the mean and median respectively to address normal and skewed distributions respectively - these variables were only missing 4% of their records.
• Removed TEAM_BATTING_HBP because of its numerous missing values (92% of records missing values for this variable.)

f1_bsb_raw <- bsb_train_raw %>% 
  mutate(t_bat_so = na_ifelse(t_bat_so, "mean"),
         t_pitch_so = na_ifelse(t_pitch_so, "med")) %>% 
  select(-t_bat_hbp)
# sapply(new_bsb_raw, function(x) sum(is.na(x)))   # count na's per column!
# View(new_bsb_raw)

Model 2 Data:

Motivation: a refinement of the first model data-set that excludes variables that may be doing ‘double-duty’:

• Replace NA with mean for: TEAM_FIELD_DP to cover 13% of missing values
• Replace NA with median for TEAM_PITCH_SO to cover 4% missing values
• Remove TEAM_BATTING_SO because it’s highly correlated with TEAM_PITCHING_SO
• Remove TEAM_BASERUN_CS and TEAM_BATTING_HBP because 44% and 92% of records are missing respectively.

f2_bsb_raw <- bsb_train_raw %>% 
  mutate(t_field_dp = na_ifelse(t_field_dp, "mean"), 
         t_pitch_so = na_ifelse(t_pitch_so, "med")) %>% 
  select(-t_bat_so, -t_baser_cs, -t_bat_hbp)

Model3 Data:

Motivation: Address all/most multicollinearity encountered and reduce the number of predictors overall to simplify the model. Variable selection between correlated predictors was based on ‘hitting-variables’ because those intuitively may result in more runs, and therefore, more wins.

• Remove TEAM_BASERUN_CS and TEAM_BATTING_HBP because 44% and 92% of records are missing respectively.
• Remove TEAM_BATTING_2B/3B/HR as TEAM_BATTING_H covers all hit types - this is an attempt to simplify the model fit
• Remove one side of the correlated pairs focusing on batting-related variables such that:
  • TEAM_BATTING_SO kept in, removed TEAM_PITCHING_SO
  • TEAM_BATTING_H kept in, removed TEAM_PITCHING_H
  • TEAM_PITCHING_HR kept in, removed TEAM_BATTING_HR as TEAM_BATTING_H covers all hits
  • TEAM_FIELDING_E kept in, removed TEAM_FIELDING_DP

f3_bsb_raw <- bsb_train_raw %>% 
  select(-t_baser_cs, -t_bat_hbp,
         -t_bat_2b, -t_bat_3b, -t_bat_hr,
         -t_pitch_so, -t_pitch_h, -t_field_dp)


## A function to pull the coefficients from a linear model and present them nicely. 
get_coefs <- function(an_lm){
  temp = an_lm$coefficients %>% as.data.frame()
  temp$Vars = rownames(temp)
  colnames(temp)[1] = "Coefs."
  rownames(temp) = NULL
  return(temp[,2:1])
}

Model Fitting

Models were fit with the respective data-sets created in the previous section. Each section will display the resulting fit’s coefficients but the discussion of the coefficients will be saved for the Model Selection section.

######## MODEL 1 FIT ########################
fitA <- lm(target_wins ~., f1_bsb_raw[,-1])
fitA_stepped <- MASS::stepAIC(fitA, trace=F)

#summary(fitA_stepped)
par(mfrow=c(1,2))
plot(fitA_stepped,which=c(2,1))

par(mfrow=c(1,1))

######## MODEL 2 FIT ########################
fitB <- lm(target_wins ~., f2_bsb_raw[,-1])
#fitB1 <- lm(target_wins+.01 ~., f2_bsb_raw[,-1])
#fitB_bc <- car::powerTransform(fitB1, family="bcPower") # 1.193639
#fitB_bc$roundlam # provides maximum estimator
#df <- cbind(lambda = fitB_bc$lambda, log_likli = fitB_bc$llik) %>% as.data.frame()
#df <- df %>% ungroup() %>% arrange(desc(log_likli)) 

#summary(fitB)

par(mfrow=c(1,2))
plot(fitB,which=c(2,1))

par(mfrow=c(1,1))

Criteria:

My selection criteria is based on comparing \(R^2\) and F-Statistics. If the performance is similar across these measures, I will select the model that is the easiest to describe of the two.

Model Selection:

I utilized the stargazer package to create the table below that compares Models 1, 2, and 3. Note that at the bottom of the table the models are compared by \(R^2\), adjusted \(R^2\), and the \(F\)-statistic.

stargazer(fitA_stepped,fitB, fitC,
          type='html',
          column.sep.width = "1pt",
          single.row = T,
          omit.stat=c('ser'),
          no.space=TRUE)

Observations from these 3 models:

• None of the 3 models have very high \(R^2\) values.
• Each model’s F-statistic indicates that there is a relationship between predictor and the response variables included in the fits

Coefficients:

Model1: The intercept (‘Constant’) is too low for the mean value of TARGET_WINS but it’s not likely that every baseball team will get ~58 wins a season. The signs of the coefficients are problematic for several cases e.g. TEAM_BASERUN_CS (caught stealing) is positive but that would be detrimental to the wins and TEAM_FIELDING_E (fielding errors) are negative but fielding errors would lead to more on-bases and therefore, more runs.

Model2: More-intuitive intercept value but similarly un-intuitive coefficient signs that mirror the issues pointed out in Model1.

Model3: This model contains the fewest variables but 2 of them still have problematic signs including TEAM_FIELDING_E and TEAM_PITCHING_BB which I would expect to be positively correlated with wins.

Based on my criteria, since the performance of Model 1 is superior compared to the three models, I have selected Model 1.

Model Evaluation:

The summary output from Model 1, above, reveals a number of challenges:

• The Intercept value of ~ 60 could be too high given that TARGET_WINS can take values as low as 0 (team loses all games) but it could also be too low given that the mean and median is \(approx\) 80.
• The coefficients have very low p-values so the observed data is extremely unlikely to to have arisen if the the null hypothesis were true (i.e. that \(H_0: \beta_i = 0\) so we can reject the null hypothesis)
• 790 observations related to the missing values in columns TEAM_BATTING_HBP
• The \(R^2\) and multiple \(R^2\) values are low (\(R^2<.70\)) at \(\approx.43\) - this indicates to me that the data isn’t very close to the fitted regression and may provide poor predictions
• The negative coefficients under TEAM_BATTING_2B and other base-gaining activities is counter-intuitive as getting on-base intuitively would result in more wins.

Below, the residual and QQ plots of Model1 appear to support a normal classification with some reservations. The residuals vs. the fitted plot appear to have a slight ‘kink’ when I expected an even, straight line. The normal QQ plot doesn’t deviate from the diagnol except in a few outlying cases.

par(mfrow=c(1,2))
plot(fitA_stepped,which=c(2,1))

par(mfrow=c(1,1))

Model CI:

Prediction:

#READ IN EVALUATION DATA
bsb_eval_raw <- read_csv(mb_eval) 

## 
## -- Column specification --------------------------------------------------------
## cols(
##   INDEX = col_double(),
##   TEAM_BATTING_H = col_double(),
##   TEAM_BATTING_2B = col_double(),
##   TEAM_BATTING_3B = col_double(),
##   TEAM_BATTING_HR = col_double(),
##   TEAM_BATTING_BB = col_double(),
##   TEAM_BATTING_SO = col_double(),
##   TEAM_BASERUN_SB = col_double(),
##   TEAM_BASERUN_CS = col_double(),
##   TEAM_BATTING_HBP = col_double(),
##   TEAM_PITCHING_H = col_double(),
##   TEAM_PITCHING_HR = col_double(),
##   TEAM_PITCHING_BB = col_double(),
##   TEAM_PITCHING_SO = col_double(),
##   TEAM_FIELDING_E = col_double(),
##   TEAM_FIELDING_DP = col_double()
## )

bsb_eval_raw <- data.table::setnames(bsb_eval_raw, 
                                  tolower(names(bsb_eval_raw[1:16])))

c1 <- col_name_updater(bsb_eval_raw, "batting", "bat")
c2 <- col_name_updater(c1, "pitching", "pitch")
c3 <- col_name_updater(c2, "fielding", "field")
c4 <- col_name_updater(c3, "baserun", "baser")
bsb_eval_raw <- col_name_updater(c4, "team", "t")

bsb_eval_raw <- bsb_eval_raw %>% 
  mutate(t_bat_so = na_ifelse(t_bat_so, "mean"),
         t_pitch_so = na_ifelse(t_pitch_so, "med")) %>% 
  select(-t_bat_hbp)

Below, I’ve used the evaluation data set, predicted the wins based on its ~200 records, and plotted the TARGET_WINS vs TEAM_BATTING_H. You can see the values wins predicted based on model 1 red in color.

# load test data into prediction
target_wins <- predict.lm(fitA_stepped, newdata=bsb_eval_raw[,-1])

# complete the test set
eval_set <- cbind(as.data.frame(target_wins),bsb_eval_raw[,-1]) %>% 
  mutate(tag='eval')

tot_set <- f1_bsb_raw[,-1] %>% 
  mutate(tag='test')

tot_set <- rbind(tot_set, eval_set)

x<-ggplot(tot_set, aes(x=t_bat_h, y=target_wins, color=tag)) +
  geom_point(alpha=.4)
x