Lab 1: Baseball

Author

Varun Srinvasa Murthy

Section B: Data Exploration

Setup

  1. Make a code chunk below and run 1+1, to make sure it is working.
1+1
[1] 2

  1. Run the chunk below to install the libraries you will need.

    # R
install.packages("ISLR")
install.packages("janitor")

# Tidy
install.packages("tidyverse")
install.packages("tidymodels")
install.packages("GGally")

  1. Make a code chunk below, and label it “libraries”. Load your libraries.

    #|: libraries-r

library("ISLR")
library("janitor")
    
Attaching package: 'janitor'
    The following objects are masked from 'package:stats':

    chisq.test, fisher.test
    #Tidy
library("tidyverse")
    ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
    ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
    library("tidymodels")
    ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
✔ broom        1.0.6      ✔ rsample      1.2.1 
✔ dials        1.2.1      ✔ tune         1.2.1 
✔ infer        1.0.7      ✔ workflows    1.1.4 
✔ modeldata    1.4.0      ✔ workflowsets 1.1.0 
✔ parsnip      1.2.1      ✔ yardstick    1.3.1 
✔ recipes      1.0.10     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()
• Dig deeper into tidy modeling with R at https://www.tmwr.org
    library("GGally")
    Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

  2. Make a code chunk below and label it “data-load”. Load the data from the ISLR/ISLP package.

    #|: data-load

data("Hitters")

  3. Look at a quick summary of the dataset. What kind of information does it contain? (Documentation of the dataset can be found here.)

    #|: summary-r

names(Hitters)
     [1] "AtBat"     "Hits"      "HmRun"     "Runs"      "RBI"       "Walks"    
 [7] "Years"     "CAtBat"    "CHits"     "CHmRun"    "CRuns"     "CRBI"     
[13] "CWalks"    "League"    "Division"  "PutOuts"   "Assists"   "Errors"   
[19] "Salary"    "NewLeague"
    summary(Hitters)
         AtBat            Hits         HmRun            Runs       
 Min.   : 16.0   Min.   :  1   Min.   : 0.00   Min.   :  0.00  
 1st Qu.:255.2   1st Qu.: 64   1st Qu.: 4.00   1st Qu.: 30.25  
 Median :379.5   Median : 96   Median : 8.00   Median : 48.00  
 Mean   :380.9   Mean   :101   Mean   :10.77   Mean   : 50.91  
 3rd Qu.:512.0   3rd Qu.:137   3rd Qu.:16.00   3rd Qu.: 69.00  
 Max.   :687.0   Max.   :238   Max.   :40.00   Max.   :130.00  

      RBI             Walks            Years            CAtBat       
 Min.   :  0.00   Min.   :  0.00   Min.   : 1.000   Min.   :   19.0  
 1st Qu.: 28.00   1st Qu.: 22.00   1st Qu.: 4.000   1st Qu.:  816.8  
 Median : 44.00   Median : 35.00   Median : 6.000   Median : 1928.0  
 Mean   : 48.03   Mean   : 38.74   Mean   : 7.444   Mean   : 2648.7  
 3rd Qu.: 64.75   3rd Qu.: 53.00   3rd Qu.:11.000   3rd Qu.: 3924.2  
 Max.   :121.00   Max.   :105.00   Max.   :24.000   Max.   :14053.0  

     CHits            CHmRun           CRuns             CRBI        
 Min.   :   4.0   Min.   :  0.00   Min.   :   1.0   Min.   :   0.00  
 1st Qu.: 209.0   1st Qu.: 14.00   1st Qu.: 100.2   1st Qu.:  88.75  
 Median : 508.0   Median : 37.50   Median : 247.0   Median : 220.50  
 Mean   : 717.6   Mean   : 69.49   Mean   : 358.8   Mean   : 330.12  
 3rd Qu.:1059.2   3rd Qu.: 90.00   3rd Qu.: 526.2   3rd Qu.: 426.25  
 Max.   :4256.0   Max.   :548.00   Max.   :2165.0   Max.   :1659.00  

     CWalks        League  Division    PutOuts          Assists     
 Min.   :   0.00   A:175   E:157    Min.   :   0.0   Min.   :  0.0  
 1st Qu.:  67.25   N:147   W:165    1st Qu.: 109.2   1st Qu.:  7.0  
 Median : 170.50                    Median : 212.0   Median : 39.5  
 Mean   : 260.24                    Mean   : 288.9   Mean   :106.9  
 3rd Qu.: 339.25                    3rd Qu.: 325.0   3rd Qu.:166.0  
 Max.   :1566.00                    Max.   :1378.0   Max.   :492.0  

     Errors          Salary       NewLeague
 Min.   : 0.00   Min.   :  67.5   A:176    
 1st Qu.: 3.00   1st Qu.: 190.0   N:146    
 Median : 6.00   Median : 425.0            
 Mean   : 8.04   Mean   : 535.9            
 3rd Qu.:11.00   3rd Qu.: 750.0            
 Max.   :32.00   Max.   :2460.0            
                 NA's   :59                
    tabyl(Hitters, League, Division)
     League  E  W
      A 85 90
      N 72 75

  4. In this lab, we will study what factors lead to a player having a higher salary. Let’s look at a plot of the distribution of the Salary variable. Make a chunk below labeled “salary-hist”, with the figure caption “Distribution of Baseball Player Salaries (1986 and 1987)” and use it to plot a histogram.

    hist(Hitters$Salary, xlab= "Annual Salary($)")

    Distribution of Baseball Player Salaries in 1986 and 1987

  5. Render your Quarto document now, to make sure everything works.

Exploratory Data Analysis

  1. Produce a scatterplot matrix of the first five columns or variables of the data, as well as the Salary variable. Do any of the variables seem to be strongly associated with each other? With Salary?

    #|: label: pairs-r

pairs(Hitters[,c(1:5, 19)])

    Hitters %>%
  select(1:5, Salary) %>%
  ggpairs()
    Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
Removed 59 rows containing missing values
Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
Removed 59 rows containing missing values
Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
Removed 59 rows containing missing values
Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
Removed 59 rows containing missing values
Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
Removed 59 rows containing missing values
    Warning: Removed 59 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 59 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 59 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 59 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 59 rows containing missing values or values outside the scale range
(`geom_point()`).
    Warning: Removed 59 rows containing non-finite outside the scale range
(`stat_density()`).

  2. Create a new qualitative variable, called BigHitter, by binning the HmRun variable. We are going to divide hitters into two groups based on whether or not they hit more than 20 home runs in a season.

    Hitters$BigHitter <- cut(Hitters$HmRun, 
                 breaks = c(-Inf, 20, +Inf), 
                 labels = c("<20 Home Runs", ">20 Home Runs"), 
                 right = FALSE)

  3. Do “Big Hitters” tend to get paid more? Plot side-by-side box-plots of BigHitter versus Salary.

    BigHitters <- Hitters[Hitters$BigHitter == ">20 Home Runs", ]
SmallHitters <- Hitters[Hitters$BigHitter == "<20 Home Runs", ]

boxplot(BigHitters$Salary, main = "Big Hitters' Salaries")

    boxplot(SmallHitters$Salary, main = "Small Hitters' Salaries")

  4. Are there more “Big Hitters” in one conference than the other? Make bar-plots to answer the question.

    barplot(table(BigHitters$League), 
        names.arg=c("A", "N"), 
        xlab = "Conference", 
        ylab = "Count", 
        main = "Big Hitters by League")

    barplot(table(SmallHitters$League), 
        names.arg=c("A", "N"), 
        xlab = "Conference", 
        ylab = "Count",
        main = "Small Hitters by League")

  5. Are there more “Big Hitters” in one division of one conference than the others? Make colored bar-plots to answer the question.

    BigHittersConferenceA = BigHitters[BigHitters$League == "A",]
barplot(table(BigHittersConferenceA$Division),
        xlab = "Division", 
        ylab = "Count", 
        main = "Big Hitters (A) by Division",
        col = rgb(0.2,0.4, 0.6, 0.6))

    BigHittersConferenceN = BigHitters[BigHitters$League == "N",]
barplot(table(BigHittersConferenceN$Division), 
        xlab = "Division", 
        ylab = "Count",
        main = "Big Hitters (N) by Division",
        col = rgb(0.2,0.4, 0.6, 0.6))

    # There are more big hitters in division E of conference A.
# There are more big hitters in division W of conference N.

“Prediction” of Salaries

  1. Run the following code to divide the dataset into a test set with 25% of the players and a training set with the remaining 75%.
## Base

set.seed(1234)
rand <- sample(1:nrow(Hitters), 100)

testing_data <- Hitters[rand,]
training_data <- Hitters[-rand,]

## Tidy

splits <- initial_split(Hitters, prop = 0.25)
testing_data <- testing(splits)
training_data <- training(splits)

  1. Find the mean salary in each league. Make a new column of the dataset called LeagueMean that contains the corresponding mean salaries.

    meanSalaryA = mean(na.omit(Hitters[Hitters$League == "A", ]$Salary))
print(meanSalaryA)
    [1] 541.9995
    Hitters[Hitters$League == "A", ]$LeagueMean = meanSalaryA
    Warning in `[<-.data.frame`(`*tmp*`, Hitters$League == "A", , value =
structure(list(: provided 22 variables to replace 21 variables
    meanSalaryN = mean(na.omit(Hitters[Hitters$League == "N", ]$Salary))
print(meanSalaryN)
    [1] 529.1175
    Hitters[Hitters$League == "N", ]$LeagueMean = meanSalaryN
    Warning in `[<-.data.frame`(`*tmp*`, Hitters$League == "N", , value =
structure(list(: provided 22 variables to replace 21 variables

  2. In the training set, make a new variable called PredError that contains the difference between each player’s salary, and the average salary for his League; i.e., the “predicted” salary.

    training_data$LeagueMean = 0
training_data[training_data$League == "A", ]$LeagueMean = meanSalaryA
training_data[training_data$League == "N", ]$LeagueMean = meanSalaryN

training_data$PredError = training_data$Salary - training_data$LeagueMean

  3. Find the average squared training error for each League

    squaredTrainingError_A = mean(na.omit(training_data[training_data$League == "A", ]$PredError)^2)
print(squaredTrainingError_A)
    [1] 267637.4
    squaredTrainingError_N = mean(na.omit(training_data[training_data$League == "N", ]$PredError)^2)
print(squaredTrainingError_N)
    [1] 80848.68

  4. In the test set, make a new variable called PredError that contains the difference between each player’s salary, and the average salary for his League from the training data; i.e., the “predicted” salary.

    testing_data$PredError = 0
testing_data[testing_data$League == "A", ]$PredError = testing_data[testing_data$League == "A", ]$Salary - meanSalaryA
testing_data[testing_data$League == "N", ]$PredError = testing_data[testing_data$League == "N", ]$Salary - meanSalaryN

  5. Find the average squared test error for each League

    squaredTestingError_A = mean(na.omit(testing_data[testing_data$League == "A", ]$PredError)^2)
print(squaredTestingError_A)
    [1] 195182.3
    squaredTestingError_N = mean(na.omit(testing_data[testing_data$League == "N", ]$PredError)^2)
print(squaredTestingError_N)
    [1] 219733.5

  6. Which was larger, the training error or the test error? Why does this make sense?

    The test error is much larger (sum of errors for both Leagues). This makes sense because the test set is 75% of the data (3x the training set).

  7. Find the average squared test error for Big Hitters versus not. Based on this, do you think it would be a good idea to include home runs in the prediction process?

    print(mean(na.omit(testing_data[testing_data$BigHitter == ">20 Home Runs", ]$PredError)^2))
    [1] 445001.8
    print(mean(na.omit(testing_data[testing_data$BigHitter == "<20 Home Runs", ]$PredError)^2))
    [1] 149172.9

Because the average squared test error is very high for Big Hitters, I don’t think it would be a good idea to include this feature. The high error might be a signal of overfitting and lead to weak generalization.

AI Statement

I used Bard to gain more perspective on the na.omit() function and to understand how to think about a high test error when considering a particular feature. I used Bard as a resource to ask questions and understand answers, and then wrote my own answer based of what I just learned. Bard was simply a useful tool, not a replacement for my learning or work in this assignment.

Section C: Linear modeling

Simple linear regression

  1. Salary by league.

    1. Fit a linear regression on your training set to predict salary from league. Report and interpret the slope value.

    2. Calculate the mean squared training error and the mean squared test error for this model. How does this compare to your analysis in Section D?

    3. Find and interpret the R-Squared value for this model.

  2. Salary by league and division.

    1. Fit a linear regression on your training set to predict salary from the interaction between league and division. Report and interpret the slope value.

    2. Why did we include the interaction term, but not each term separately?

    3. Calculate the mean squared training error and the mean squared test error for this model. How does this compare to the model with just league?

    4. Find and interpret the R-Squared value for this model.

  3. Salary by home runs.

    1. Fit a linear regression on your training set to predict salary from number of home runs. Report and interpret the slope value.

    2. What would be the advantages or disadvantages of using HmRun as the predictor, as we did here, rather than BigHitter?

    3. Calculate the mean squared training error and the mean squared test error for this model.

    4. Find and interpret the R-Squared value for this model.

    5. Plot and comment on the residuals of this model.

  4. Which model do you think is better: Q1, Q2, or Q3? Give two different reasons for your answer.

Overfitting

  1. Make four new columns in your original Hitters dataset: HmRun_2, HmRun_3, HmRun_4, and HmRun_5.

  2. Fit a linear regression using HmRun and all four new columns as predictors. Report (but do not interpret) the resulting coefficients.

  3. Report and interpret the R-squared value of this model.

  4. Report the training error and test error of this model. Give an explanation for why this makes sense.

  5. Which model do you think is better, this one or the one using just HmRun?

AI Statement

In this section, make brief notes about which of your code and/or interpretations was assisted by genAI, and what information you got from the genAI.

Section D: Multiple Regression

  1. Fit and report metrics for the following models:

  • a model using the first 6 columns of the Hitters dataset as predictors, plus the interaction of league and division.

  • a model using the next 6 columns (i.e. CAtBat through CWalks) as predictors, plus the interaction of league and division.

  • a model using all of the first 12 columns as predictors, plus the interaction of league and division.

  • (optional, no grade points) a model using more interactions, different predictor combinations, or more predictors that you create.

  1. Based on the above, which model best predicts salary?

  2. Plot and comment on the residuals of your best model.

  3. Report and interpret the three most significant coefficients from your best model, including appropriate plots of each predictor with Salary.

AI Statement

In this section, make brief notes about which of your code and/or interpretations was assisted by genAI, and what information you got from the genAI.

Final Report

Make a copy of this Quarto document and save it as Lab1Report_YOURLASTNAME_YOURFIRSTNAME.qmd. Then delete, edit, and rearrange the contents of Sections A-D to construct a brief but clean Final Report on the question of what factors influence baseball player salaries. Do this by:

  1. Setting code-fold: true in your YAML.
  2. Choosing a few exploratory plots from Section B that help inform your final model.
  3. Briefly addressing the question of whether different conferences/divisions have different salaries, using visuals from Section B and a model from Section C.
  4. Giving a full report and interpretation of your best model from Section D.
  5. Writing a (very short) conclusion that states the main real-world takeaways from the analysis.
  6. Adding a (very short) ethics statement, giving any concerns or cautions you have about this analysis.
  7. Adding a (very short) AI disclosure statement, attributing help you got from AI.

This report does not need to be formal like an academic paper, and it does not have to contain any more written text than is necessary to explain your analysis and report your results. It should, however, have informative section headers and titles; e.g. “Factors that contribute to baseball salaries in 1987” not “Lab 1”.

In other words, you do not need to add extra content or analysis to make this report; you are simply cleaning everything up into a nice looking summary.