find_data_file <- function(candidates) {
  hit <- candidates[file.exists(candidates)]
  if (length(hit) == 0) {
    stop("Could not find the data file. Tried:\n  ", paste(candidates, collapse = "\n  "),
         "\nEasiest fix: put this .Rmd in the same folder as the CSV files.")
  }
  hit[1]
}

train_raw <- read.csv(find_data_file(c(
  "moneyball-training-data-1-1.csv",
  "training.csv",
  "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 1/moneyball-training-data-1-1.csv"
)))
eval_raw <- read.csv(find_data_file(c(
  "moneyball-evaluation-data-1-1.csv",
  "evaluation.csv",
  "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 1/moneyball-evaluation-data-1-1.csv"
)))

Introduction

The goal of this assignment is to predict TARGET_WINS, the number of games a professional baseball team wins in a 162-game season, from 15 team-level batting, baserunning, pitching, and fielding statistics. The data set covers 2276 team-seasons from 1871 to 2006, each already normalized to a 162-game schedule so that teams from different eras can be compared directly.

The write-up below follows the four stages requested in the assignment: first, an exploratory pass over the raw data to understand its shape, scale, and quality problems; second, a data preparation stage that fixes the missing values and data entry errors uncovered in exploration; third, three manually specified multiple linear regression models built on the cleaned data; and fourth, a comparison of those three models on fit, parsimony, and diagnostic grounds, ending in a single selected model that is used to generate predictions for the evaluation set.

Data Exploration

Overview

The training set contains 2276 team-seasons and 15 candidate predictors (excluding INDEX and the response, TARGET_WINS). The evaluation set contains 259 records with the same predictors and no response, which is what we ultimately need to generate predictions for. Every statistic has already been normalized to a 162-game season by the data providers, so we are comparing team-seasons on an equal footing regardless of the era they come from (1871-2006). All variables are integer counts; the full list, along with each variable’s distribution, appears in the summary statistics table below.

Summary Statistics

num_vars <- setdiff(names(train_raw), "INDEX")
desc_tbl <- describe(train_raw[, num_vars])[, c("n","mean","sd","median","min","max","skew","kurtosis")]
kable(round(desc_tbl, 2), caption = "Summary statistics - raw training data")

Summary statistics - raw training data
	n	mean	sd	median	min	max	skew	kurtosis
TARGET_WINS	2276	80.79	15.75	82.0	0	146	-0.40	1.03
TEAM_BATTING_H	2276	1469.27	144.59	1454.0	891	2554	1.57	7.28
TEAM_BATTING_2B	2276	241.25	46.80	238.0	69	458	0.22	0.01
TEAM_BATTING_3B	2276	55.25	27.94	47.0	0	223	1.11	1.50
TEAM_BATTING_HR	2276	99.61	60.55	102.0	0	264	0.19	-0.96
TEAM_BATTING_BB	2276	501.56	122.67	512.0	0	878	-1.03	2.18
TEAM_BATTING_SO	2174	735.61	248.53	750.0	0	1399	-0.30	-0.32
TEAM_BASERUN_SB	2145	124.76	87.79	101.0	0	697	1.97	5.49
TEAM_BASERUN_CS	1504	52.80	22.96	49.0	0	201	1.98	7.62
TEAM_BATTING_HBP	191	59.36	12.97	58.0	29	95	0.32	-0.11
TEAM_PITCHING_H	2276	1779.21	1406.84	1518.0	1137	30132	10.33	141.84
TEAM_PITCHING_HR	2276	105.70	61.30	107.0	0	343	0.29	-0.60
TEAM_PITCHING_BB	2276	553.01	166.36	536.5	0	3645	6.74	96.97
TEAM_PITCHING_SO	2174	817.73	553.09	813.5	0	19278	22.17	671.19
TEAM_FIELDING_E	2276	246.48	227.77	159.0	65	1898	2.99	10.97
TEAM_FIELDING_DP	1990	146.39	26.23	149.0	52	228	-0.39	0.18

A few things jump out immediately:

TARGET_WINS is well-behaved: mean 80.8, median 82, roughly symmetric (skew close to 0), and bounded sensibly between 0 and 146 wins out of 162 games.
TEAM_PITCHING_H, TEAM_PITCHING_BB, and TEAM_PITCHING_SO have enormous skew and kurtosis (skew of 10, 6.7, and 22 respectively). A team cannot plausibly allow 30,132 hits in a 162-game season - the batting-side hits variable tops out at 2,554 for comparison. These are almost certainly data entry errors and need to be addressed in Data Preparation.
TEAM_FIELDING_E is also right-skewed (skew ~3), consistent with a handful of very poor defensive teams (likely older/19th-century seasons) rather than data errors, but still worth capping to protect the regression from a few high-leverage points.

Missing Data

# I am using the visdat package as suggested for missing-data visualization.
# We use it (auto-install on first knit); otherwise fall back to an
# equivalent ggplot missingness map so the document still knits anywhere.
if (!requireNamespace("visdat", quietly = TRUE)) {
  try(install.packages("visdat", repos = "https://cloud.r-project.org"), silent = TRUE)
}

if (requireNamespace("visdat", quietly = TRUE)) {
  visdat::vis_miss(train_raw %>% select(-INDEX)) +
    labs(title = "Missingness Map - Raw Training Data") +
    theme(axis.text.x = element_text(angle = 90, size = 7))
} else {
  miss_map <- train_raw %>% select(-INDEX) %>%
    mutate(.row = row_number()) %>%
    pivot_longer(-.row, names_to = "variable", values_to = "value") %>%
    mutate(missing = is.na(value))
  ggplot(miss_map, aes(x = variable, y = .row, fill = missing)) +
    geom_raster() +
    scale_fill_manual(values = c("FALSE" = "grey85", "TRUE" = "#2d2d2d"),
                      labels = c("Present", "Missing"), name = NULL) +
    scale_y_reverse() +
    theme_minimal(base_size = 10) +
    theme(axis.text.x = element_text(angle = 90, size = 7, hjust = 1)) +
    labs(title = "Missingness Map - Raw Training Data",
         x = NULL, y = "Observation")
}

The missingness map makes the pattern immediately visible without needing a separate percentage table. TEAM_BATTING_HBP is missing for roughly 92% of records; at that rate, any imputation would be mostly fabricated data, so it is dropped rather than imputed (see Data Preparation). TEAM_BASERUN_CS is missing about a third of the time (34%), which is high but still recoverable with median imputation plus a missing-value flag so the model can learn whether “missingness itself” carries signal. The remaining affected variables (TEAM_FIELDING_DP, TEAM_BASERUN_SB, TEAM_BATTING_SO, TEAM_PITCHING_SO) are each missing under 13% of the time, which routine median imputation handles well. The map also shows the missingness is blocky rather than scattered: the same rows tend to be missing several baserunning/strikeout statistics at once, consistent with entire early seasons where those statistics simply were not recorded.

Distributions of Key Variables

# Focus on the response plus the most analytically important predictors,
# rather than cramming all 15 variables into unreadable postage stamps
key_vars <- c("TARGET_WINS", "TEAM_BATTING_H", "TEAM_BATTING_BB",
               "TEAM_BATTING_HR", "TEAM_PITCHING_H", "TEAM_FIELDING_E")
train_long <- train_raw %>% select(all_of(key_vars)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  mutate(variable = factor(variable, levels = key_vars))

ggplot(train_long, aes(x = value)) +
  geom_histogram(bins = 40, fill = "#4a7a8a", color = "white", linewidth = 0.1, na.rm = TRUE) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  theme_minimal(base_size = 11) +
  labs(title = "Distributions of the Response and Key Predictors",
       subtitle = "Raw training data, before any cleaning",
       x = "Value", y = "Count")

Three comments on these distributions:

TARGET_WINS is roughly bell-shaped and centered near 81 wins, which is exactly what one would expect: 81 is half of a 162-game season, and in any season the league as a whole must win half its games.
TEAM_BATTING_H, TEAM_BATTING_BB, and TEAM_BATTING_HR look like ordinary counting statistics. Home runs show a mild bimodality (a cluster of low-HR team-seasons alongside a modern-era cluster), consistent with the data spanning 1871-2006 across very different offensive eras.
TEAM_PITCHING_H and TEAM_FIELDING_E are the problem children: both are so stretched by extreme right-tail values that the bulk of the distribution collapses into the leftmost bars. For TEAM_PITCHING_H the extremes are physically impossible (30,132 hits allowed in a 162-game season) and point to data entry errors; for TEAM_FIELDING_E the tail is merely implausible and likely reflects genuinely sloppy 19th-century defense. Both are handled in Data Preparation.

Correlation with Target Wins

# cor_tbl is computed for inline reference in the text below; the heatmap is
# the single printed correlation display (a separate correlation table would
# be redundant with the annotated heatmap)
cor_tbl <- train_raw %>% select(-INDEX) %>%
  cor(use = "pairwise.complete.obs") %>%
  { .[, "TARGET_WINS"] } %>%
  sort(decreasing = TRUE) %>%
  round(3)

cor_mat <- train_raw %>% select(-INDEX) %>% cor(use = "pairwise.complete.obs")
corrplot(cor_mat, method = "color", type = "upper", order = "hclust",
         tl.col = "black", tl.cex = 0.85, tl.srt = 45,
         addCoef.col = "black", number.cex = 0.62,
         col = colorRampPalette(c("#b2182b","white","#2166ac"))(200),
         title = "Correlation Matrix - Raw Training Data (pairwise complete)",
         mar = c(0,0,2,0), diag = FALSE)

TEAM_BATTING_H is the single strongest predictor of wins (r = 0.389), which makes intuitive sense - teams that get more hits score more runs. Notably, two variables come in with the wrong sign relative to the assignment’s stated theoretical direction: TEAM_PITCHING_SO (strikeouts by pitchers, theoretically positive) correlates at -0.078, and TEAM_FIELDING_DP (double plays, theoretically positive) correlates at -0.035.

The correlation heatmap also surfaces a serious multicollinearity problem that needs to be resolved before modeling, not just noted: TEAM_BATTING_HR and TEAM_PITCHING_HR correlate at 0.969, and are exactly equal in 815 of 2276 rows. A team’s own home runs hit and its pitchers’ home runs allowed have no mechanical reason to move together this tightly; the most likely explanation is a shared park factor (home run friendly ballparks inflate both numbers for whichever team plays there), but whatever the cause, including both raw variables in the same regression would badly inflate their variance inflation factors. This is addressed directly in Build Models by dropping one of the pair rather than carrying the collinearity into the final model. The batting and pitching “hits” volume statistics are also correlated (r = 0.303) but far less severely, and are left in as separate variables.

scatter_vars <- c("TEAM_BATTING_H", "TEAM_PITCHING_H")
scatter_long <- train_raw %>% select(TARGET_WINS, all_of(scatter_vars)) %>%
  pivot_longer(all_of(scatter_vars), names_to = "variable", values_to = "value")

ggplot(scatter_long, aes(x = value, y = TARGET_WINS)) +
  geom_point(alpha = 0.25, color = "#4a7a8a") +
  geom_smooth(method = "lm", se = FALSE, color = "#a0522d") +
  facet_wrap(~variable, scales = "free_x") +
  theme_minimal(base_size = 11) +
  labs(title = "TARGET_WINS vs. the two strongest raw predictors", x = NULL, y = "Wins")

The TEAM_BATTING_H panel shows a clean, fairly tight positive trend, consistent with its strong correlation. The TEAM_PITCHING_H panel tells a visual story that the correlation number alone does not: a dense, unremarkable cloud of points sits below roughly 3,000 hits allowed, and then a handful of extreme points (the same data-entry errors flagged above, some past 20,000 hits allowed) stretch the x-axis out and pull the fitted line’s slope down. This is the clearest visual evidence that the raw negative correlation for TEAM_PITCHING_H is being driven by a small number of corrupted rows rather than the bulk of the data, a point revisited quantitatively once the model coefficients are examined in Build Models.

Data Preparation

Based on the exploration above, the following cleaning pipeline is applied identically to the training and evaluation sets (using training-set statistics for any imputation values, so the evaluation set is not allowed to leak information back into the model).

# Fences learned from TRAINING data only, then applied to both sets
learn_fences <- function(df, vars) {
  fences <- list()
  for (v in vars) {
    q1 <- quantile(df[[v]], .25, na.rm = TRUE)
    q3 <- quantile(df[[v]], .75, na.rm = TRUE)
    iqr <- q3 - q1
    fences[[v]] <- q3 + 3 * iqr   # Tukey's "far outlier" upper fence
  }
  fences
}

learn_medians <- function(df, vars) {
  sapply(vars, function(v) median(df[[v]], na.rm = TRUE))
}

clean_moneyball <- function(df, medians, fences, is_training = TRUE) {

  df <- df %>% select(-any_of("TEAM_BATTING_HBP"))  # ~92% missing, drop

  # Median-impute the moderately-missing variables and keep a flag for each,
  # so the model can separate "value was X" from "value was unknown"
  flag_vars <- c("TEAM_BATTING_SO", "TEAM_BASERUN_SB", "TEAM_BASERUN_CS",
                 "TEAM_PITCHING_SO", "TEAM_FIELDING_DP")
  for (v in flag_vars) {
    df[[paste0(v, "_MISS")]] <- as.integer(is.na(df[[v]]))
    df[[v]][is.na(df[[v]])] <- medians[[v]]
  }

  # Winsorize (cap, don't delete) anything past the Tukey far-outlier fence
  # learned from the training set, so implausible values stop distorting the fit
  cap_vars <- names(fences)
  for (v in cap_vars) {
    df[[v]] <- pmin(df[[v]], fences[[v]])
  }

  # Feature engineering
  df$TEAM_BATTING_1B     <- df$TEAM_BATTING_H - df$TEAM_BATTING_2B -
                             df$TEAM_BATTING_3B - df$TEAM_BATTING_HR
  df$TEAM_BASERUN_SB_NET <- df$TEAM_BASERUN_SB - df$TEAM_BASERUN_CS

  df
}

prep_vars  <- setdiff(names(train_raw), c("INDEX", "TARGET_WINS", "TEAM_BATTING_HBP"))
medians    <- learn_medians(train_raw, c("TEAM_BATTING_SO","TEAM_BASERUN_SB",
                                          "TEAM_BASERUN_CS","TEAM_PITCHING_SO",
                                          "TEAM_FIELDING_DP"))
fences     <- learn_fences(train_raw, c("TEAM_BATTING_H","TEAM_BATTING_3B",
                                         "TEAM_BASERUN_SB","TEAM_BASERUN_CS",
                                         "TEAM_PITCHING_H","TEAM_PITCHING_BB",
                                         "TEAM_PITCHING_SO","TEAM_FIELDING_E"))

train <- clean_moneyball(train_raw, medians, fences)
evalu <- clean_moneyball(eval_raw,  medians, fences)

Dropped: TEAM_BATTING_HBP (92% missing - imputing this would mean fabricating the large majority of its values).

Imputed (median) + flagged: TEAM_BATTING_SO, TEAM_BASERUN_SB, TEAM_BASERUN_CS, TEAM_PITCHING_SO, TEAM_FIELDING_DP. Median imputation was chosen over mean imputation because several of these variables are right-skewed, and the median is more robust to that skew. A companion _MISS indicator flag is kept for each so the model can pick up any information content in missingness itself (e.g. earlier seasons may not have tracked caught-stealing at all, which would show up as missingness correlated with era, which is itself correlated with the run-scoring environment).

Outlier treatment: Tukey’s classic “far outlier” rule (\(Q3 + 3 \times IQR\)) is applied uniformly to every continuous predictor with implausible maxima, and any value beyond that fence is capped (winsorized) to the fence value rather than deleted, so we don’t throw away otherwise-valid information in the rest of that row. This is what brings TEAM_PITCHING_H (previously up to 30,132 hits allowed - impossible in a 162-game season) back into a believable range. The same rule is applied consistently to TEAM_BATTING_H, TEAM_BATTING_3B, TEAM_BASERUN_SB, TEAM_BASERUN_CS, TEAM_PITCHING_BB, TEAM_PITCHING_SO, and TEAM_FIELDING_E. The rule affects a small, variable-specific fraction of rows, not the whole data set:

capped_counts <- sapply(names(fences), function(v) sum(train_raw[[v]] > fences[[v]], na.rm = TRUE))
capped_tbl <- data.frame(variable = names(capped_counts), n_capped = capped_counts,
                          pct_capped = round(100 * capped_counts / nrow(train_raw), 1),
                          row.names = NULL) %>%
  arrange(desc(n_capped))
kable(capped_tbl, caption = "Observations capped per variable (Tukey 3xIQR upper fence)")

Observations capped per variable (Tukey 3xIQR upper fence)
variable	n_capped	pct_capped
TEAM_FIELDING_E	176	7.7
TEAM_PITCHING_H	120	5.3
TEAM_BASERUN_SB	27	1.2
TEAM_BATTING_H	19	0.8
TEAM_BASERUN_CS	18	0.8
TEAM_PITCHING_BB	16	0.7
TEAM_PITCHING_SO	9	0.4
TEAM_BATTING_3B	4	0.2

TEAM_FIELDING_E and TEAM_PITCHING_H have the most capped rows, consistent with them having the worst skew and kurtosis in the summary statistics table; most other variables are only lightly touched.

Feature engineering:

TEAM_BATTING_1B = hits that were singles = TEAM_BATTING_H minus doubles, triples, and home runs. This separates “productive” extra-base power from ordinary singles, since TEAM_BATTING_H alone bundles both together.
TEAM_BASERUN_SB_NET = TEAM_BASERUN_SB - TEAM_BASERUN_CS, i.e. net bases stolen after accounting for the outs given away by getting caught. Stolen bases and caught-stealing are highly correlated with each other (aggressive base-running teams do more of both), so this single net measure is more interpretable than including both raw counts.

clean_vars <- setdiff(names(train), c("INDEX","TARGET_WINS"))
kable(round(describe(train[, clean_vars])[, c("n","mean","sd","median","min","max")], 2),
      caption = "Summary statistics - cleaned/final training variables")

Summary statistics - cleaned/final training variables
	n	mean	sd	median	min	max
TEAM_BATTING_H	2276	1467.38	135.15	1454.0	891	2000
TEAM_BATTING_2B	2276	241.25	46.80	238.0	69	458
TEAM_BATTING_3B	2276	55.22	27.79	47.0	0	186
TEAM_BATTING_HR	2276	99.61	60.55	102.0	0	264
TEAM_BATTING_BB	2276	501.56	122.67	512.0	0	878
TEAM_BATTING_SO	2276	736.25	242.91	750.0	0	1399
TEAM_BASERUN_SB	2276	122.35	80.69	101.0	0	426
TEAM_BASERUN_CS	2276	51.26	17.27	49.0	0	134
TEAM_PITCHING_H	2276	1610.26	296.41	1518.0	1137	2473
TEAM_PITCHING_HR	2276	105.70	61.30	107.0	0	343
TEAM_PITCHING_BB	2276	547.79	121.06	536.5	0	1016
TEAM_PITCHING_SO	2276	801.58	252.02	813.5	0	2027
TEAM_FIELDING_E	2276	225.21	154.20	159.0	65	616
TEAM_FIELDING_DP	2276	146.72	24.54	149.0	52	228
TEAM_BATTING_SO_MISS	2276	0.04	0.21	0.0	0	1
TEAM_BASERUN_SB_MISS	2276	0.06	0.23	0.0	0	1
TEAM_BASERUN_CS_MISS	2276	0.34	0.47	0.0	0	1
TEAM_PITCHING_SO_MISS	2276	0.04	0.21	0.0	0	1
TEAM_FIELDING_DP_MISS	2276	0.13	0.33	0.0	0	1
TEAM_BATTING_1B	2276	1071.30	119.87	1050.0	709	1730
TEAM_BASERUN_SB_NET	2276	71.09	78.13	52.0	-54	377

No missing values remain in the modeling variables:

cat("Remaining NAs in training set:", sum(is.na(train %>% select(-TARGET_WINS))), "\n")

## Remaining NAs in training set: 0

Build Models

Three models are built manually, moving from a straightforward “everything theory suggests should matter” baseline to more deliberate variable and transformation choices.

Model 1 - Full theory-driven model

Every cleaned predictor is included, mirroring the assignment’s theoretical table directly. The specification contains 18 predictors (of which 17 are ultimately estimated; one missingness flag is dropped automatically for perfect collinearity, discussed at the regression table below). The estimating equation is written compactly in summation form, where \(x_{ki}\) denotes the \(k\)-th cleaned predictor (batting, baserunning, pitching, and fielding statistics plus the five _MISS indicator flags) for team-season \(i\):

\[\text{WINS}_i = \beta_0 + \sum_{k=1}^{18} \beta_k \, x_{ki} + \varepsilon_i\]

# Keep every cleaned variable, but swap raw H and raw SB/CS for their
# engineered replacements (1B decomposition and net stolen bases) so we are
# not modeling the same information twice.
# Full summary() output is omitted for all three models to keep the report
# readable; the consolidated regression table below reports every coefficient,
# standard error, significance level, and fit statistic in one place.
m1_vars <- setdiff(names(train), c("INDEX","TARGET_WINS","TEAM_BATTING_H",
                                    "TEAM_BASERUN_SB","TEAM_BASERUN_CS"))
m1 <- lm(TARGET_WINS ~ ., data = train[, c("TARGET_WINS", m1_vars)])

Model 2 - Reduced model (correlation + multicollinearity driven)

Model 1’s full variable set includes several highly collinear pairs (batting and pitching “volume” stats move together, and the _MISS flags add little). Model 2 keeps only variables with a defensible individual relationship to wins and lower mutual collinearity: hitting for extra bases, controlling the free bases given away, limiting runs allowed, and defense (errors, double plays). TEAM_PITCHING_HR is deliberately excluded here: as shown in Data Exploration, it correlates at 0.97 with TEAM_BATTING_HR, so including both inflates their VIFs into the 30-40 range with almost no gain in explanatory power. TEAM_BATTING_HR is kept because the assignment’s theory table treats it as the primary power-hitting metric, and TEAM_PITCHING_H (total hits allowed) is left in the model to continue representing pitching quality.

The estimating equation, with all variables in their cleaned (capped/imputed) form:

\[ \begin{aligned} \text{WINS}_i = \beta_0 &+ \beta_1 \text{1B}_i + \beta_2 \text{2B}_i + \beta_3 \text{3B}_i + \beta_4 \text{HR}_i + \beta_5 \text{BB}_i + \beta_6 \text{SB}^{net}_i \\ &+ \beta_7 \text{HitsAllowed}_i + \beta_8 \text{BBAllowed}_i + \beta_9 \text{Errors}_i + \beta_{10} \text{DP}_i + \varepsilon_i \end{aligned} \]

# Drop the strikeout variables and all _MISS flags; keep one representative
# variable per baseball "skill" to reduce collinearity versus Model 1.
# TEAM_PITCHING_HR is dropped specifically because it is nearly collinear
# with TEAM_BATTING_HR (r = 0.97, see Data Exploration)
m2 <- lm(TARGET_WINS ~ TEAM_BATTING_1B + TEAM_BATTING_2B + TEAM_BATTING_3B +
           TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB_NET +
           TEAM_PITCHING_H + TEAM_PITCHING_BB +
           TEAM_FIELDING_E + TEAM_FIELDING_DP,
         data = train)

Model 3 - Transformed / parsimonious model

TEAM_FIELDING_E remains right-skewed even after capping, so it is log-transformed here. TEAM_BATTING_SO and TEAM_PITCHING_SO remain excluded, as in Model 2: Model 1 (the only specification containing them) showed them contributing little while carrying the “wrong” sign - rather than force strikeouts into the model, Model 3 tests whether a leaner, more parsimonious specification performs just as well.

The estimating equation is identical to Model 2’s except that raw fielding errors enter in logs:

\[ \begin{aligned} \text{WINS}_i = \beta_0 &+ \beta_1 \text{1B}_i + \beta_2 \text{2B}_i + \beta_3 \text{3B}_i + \beta_4 \text{HR}_i + \beta_5 \text{BB}_i + \beta_6 \text{SB}^{net}_i \\ &+ \beta_7 \text{HitsAllowed}_i + \beta_8 \text{BBAllowed}_i + \beta_9 \ln(\text{Errors}_i) + \beta_{10} \text{DP}_i + \varepsilon_i \end{aligned} \]

# Log-transform fielding errors, which stayed right-skewed even after capping.
# TEAM_PITCHING_HR remains excluded here for the same collinearity reason as Model 2
train$log_FIELD_E <- log(train$TEAM_FIELDING_E)
evalu$log_FIELD_E <- log(evalu$TEAM_FIELDING_E)

m3 <- lm(TARGET_WINS ~ TEAM_BATTING_1B + TEAM_BATTING_2B + TEAM_BATTING_3B +
           TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB_NET +
           TEAM_PITCHING_H + TEAM_PITCHING_BB +
           log_FIELD_E + TEAM_FIELDING_DP,
         data = train)

Regression Results

The table below places all three models side by side, one row per variable, so the coefficients (and which variables each model even includes) can be compared directly rather than scrolling back through three separate summary() printouts.

# Pull tidy coefficient tables (estimate, SE, p-value) from each model and
# format each as "estimate (SE) stars" so the single table carries enough
# information to judge both direction and precision, not just sign
stars <- function(p) {
  ifelse(is.na(p), "", ifelse(p < 0.001, "***", ifelse(p < 0.01, "**", ifelse(p < 0.05, "*", ""))))
}

tidy_one <- function(model, label) {
  broom::tidy(model) %>%
    filter(term != "(Intercept)") %>%
    mutate(fmt = ifelse(is.na(estimate), NA_character_,
                         sprintf("%.4f (%.4f)%s", estimate, std.error, stars(p.value)))) %>%
    select(term, fmt) %>%
    rename(!!label := fmt)
}

coef_wide <- tidy_one(m1, "M1") %>%
  full_join(tidy_one(m2, "M2"), by = "term") %>%
  full_join(tidy_one(m3, "M3"), by = "term") %>%
  rename(variable = term)

# Append fit statistics as bottom rows so the table is self-contained and
# doesn't require flipping to a separate table to see overall model quality
fit_row <- function(label, fmt_fun) {
  data.frame(variable = label, M1 = fmt_fun(m1), M2 = fmt_fun(m2), M3 = fmt_fun(m3))
}
rmse_inline <- function(model) sqrt(mean(residuals(model)^2))

fit_stats <- bind_rows(
  fit_row("N",       function(m) as.character(nobs(m))),
  fit_row("R2",       function(m) sprintf("%.4f", summary(m)$r.squared)),
  fit_row("Adj. R2",  function(m) sprintf("%.4f", summary(m)$adj.r.squared)),
  fit_row("RMSE",     function(m) sprintf("%.3f", rmse_inline(m)))
)

coef_wide_full <- bind_rows(coef_wide, fit_stats)

kable(coef_wide_full, caption = "Coefficient comparison across Models 1-3: estimate (std. error), significance * p<0.05 ** p<0.01 *** p<0.001, with fit statistics in the bottom rows. Blank = variable not in that model.")

Coefficient comparison across Models 1-3: estimate (std. error), significance * p<0.05 ** p<0.01 *** p<0.001, with fit statistics in the bottom rows. Blank = variable not in that model.
variable	M1	M2	M3
TEAM_BATTING_2B	0.0389 (0.0075)***	0.0166 (0.0079)*	0.0099 (0.0078)
TEAM_BATTING_3B	0.1609 (0.0155)***	0.1368 (0.0159)***	0.1610 (0.0162)***
TEAM_BATTING_HR	-0.0018 (0.0299)	0.0800 (0.0084)***	0.0623 (0.0087)***
TEAM_BATTING_BB	0.0455 (0.0071)***	0.0493 (0.0064)***	0.0518 (0.0063)***
TEAM_BATTING_SO	-0.0250 (0.0060)***	NA	NA
TEAM_PITCHING_H	-0.0140 (0.0031)***	0.0165 (0.0025)***	0.0161 (0.0023)***
TEAM_PITCHING_HR	0.1121 (0.0267)***	NA	NA
TEAM_PITCHING_BB	-0.0154 (0.0063)*	-0.0281 (0.0052)***	-0.0289 (0.0052)***
TEAM_PITCHING_SO	0.0100 (0.0045)*	NA	NA
TEAM_FIELDING_E	-0.0844 (0.0056)***	-0.0359 (0.0042)***	NA
TEAM_FIELDING_DP	-0.1074 (0.0139)***	-0.1169 (0.0130)***	-0.1250 (0.0128)***
TEAM_BATTING_SO_MISS	7.3571 (1.5284)***	NA	NA
TEAM_BASERUN_SB_MISS	42.3998 (2.2970)***	NA	NA
TEAM_BASERUN_CS_MISS	0.8648 (0.8925)	NA	NA
TEAM_PITCHING_SO_MISS	NA	NA	NA
TEAM_FIELDING_DP_MISS	7.4256 (1.6624)***	NA	NA
TEAM_BATTING_1B	0.0548 (0.0044)***	0.0329 (0.0038)***	0.0354 (0.0038)***
TEAM_BASERUN_SB_NET	0.0747 (0.0057)***	0.0339 (0.0048)***	0.0330 (0.0046)***
log_FIELD_E	NA	NA	-12.8257 (1.2002)***
N	2276	2276	2276
R2	0.3979	0.3004	0.3129
Adj. R2	0.3934	0.2973	0.3099
RMSE	12.220	13.172	13.054

Note that TEAM_PITCHING_SO_MISS shows NA for Model 1 rather than a coefficient: TEAM_BATTING_SO and TEAM_PITCHING_SO are missing for exactly the same 102 rows in the training data, so their _MISS flags are perfectly collinear and lm() automatically drops one as redundant. This is another symptom of the same issue discussed below: the _MISS flags carry real signal, but that signal is tangled up with data-collection artifacts rather than clean baseball mechanics.

Most coefficients line up with the theoretical expectations in the assignment: more hits, doubles, triples, home runs, and walks are all associated with more wins, and more pitching walks and fielding errors allowed are associated with fewer wins. Three results are counter-intuitive or otherwise not what the assignment’s theory table would predict, and are worth calling out explicitly rather than glossing over:

TEAM_FIELDING_DP (double plays): the assignment states a positive theoretical effect, but the correlation table showed a slightly negative relationship, and it stays essentially flat/negative across all three models. One plausible story: teams that induce a lot of double plays are often teams that also allow a lot of runners on base in the first place (you can’t turn a double play without a baserunner to erase), so this variable may be proxying for “opponents reach base often” rather than pure defensive skill. It is kept in the model despite the sign, both because it is a required theoretical variable and because its effect size is small and not deeply damaging to the model’s interpretability.
Strikeout variables (TEAM_BATTING_SO, TEAM_PITCHING_SO): both were excluded from Models 2 and 3. TEAM_PITCHING_SO is theoretically supposed to help wins (more strikeouts by your pitchers = fewer balls in play = fewer hits allowed), but it was one of the most corrupted variables in the raw data (skew of 22 before capping), and in Model 1, the only specification that includes the strikeout variables, both added negligible explanatory power. Given the data quality concerns, it seemed safer to exclude it than to trust a coefficient built substantially on capped/imputed values.
TEAM_PITCHING_HR (home runs allowed): this is not a sign problem but a redundancy problem, dropped from Models 2 and 3 because it correlates at 0.97 with TEAM_BATTING_HR (Data Exploration). Keeping both in the same regression pushed their VIFs to 35-41; dropping TEAM_PITCHING_HR brings every VIF in Model 3 back under 9 with almost no change in adjusted \(R^2\), which is a much better trade than keeping a theoretically “required” variable whose own coefficient could not be trusted.
TEAM_PITCHING_H (hits allowed): this one is genuinely surprising, and worth sitting with rather than explaining away. Its coefficient is small but positive in both Model 2 and Model 3, the opposite of the assignment’s stated theoretical direction. Digging into why: the raw, uncapped correlation between TEAM_PITCHING_H and wins is negative (r = -0.11, matching theory), but that negative correlation is almost entirely an artifact of the 120 data-entry-error rows identified in Data Exploration (some teams supposedly allowed over 10,000 hits). Once those rows are capped, the Pearson correlation flips to +0.10, and the outlier-robust Spearman correlation, which never depended on the capping choice in the first place, agrees at +0.21. In other words, this is not a modeling artifact from the regression; the cleaned data itself shows a mildly positive relationship between hits allowed and wins. It is left in the model as-is, both because the effect size is small enough that it does not meaningfully change predictions, and because dropping a variable simply because its sign is inconvenient would be worse practice than reporting an honest, if unexpected, result.

Because the winsorization step changes this coefficient’s sign relative to the raw data, the fairest test is a sensitivity analysis: refit the Model 3 specification under three different outlier treatments and see whether the positive sign is an artifact of the specific capping rule chosen.

# Sensitivity check: does the TEAM_PITCHING_H sign depend on the capping rule?
# (a) capped at the Tukey 3xIQR fence (the treatment used in this report)
# (b) raw data, no capping at all (corrupted rows left in)
# (c) no capping, but rows with physically impossible values dropped entirely
#     (hits allowed above the max plausible ~2,554; strikeouts above ~1,800)
m3_spec <- TARGET_WINS ~ TEAM_BATTING_1B + TEAM_BATTING_2B + TEAM_BATTING_3B +
  TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB_NET +
  TEAM_PITCHING_H + TEAM_PITCHING_BB + log_FIELD_E + TEAM_FIELDING_DP

build_variant <- function(df, cap = TRUE, drop_impossible = FALSE) {
  out <- df %>% select(-any_of("TEAM_BATTING_HBP"))
  for (v in names(medians)) {
    out[[v]][is.na(out[[v]])] <- medians[[v]]
  }
  if (cap) for (v in names(fences)) out[[v]] <- pmin(out[[v]], fences[[v]])
  if (drop_impossible) out <- out %>%
      filter(TEAM_PITCHING_H <= 2554, TEAM_PITCHING_SO <= 1800)
  out$TEAM_BATTING_1B <- out$TEAM_BATTING_H - out$TEAM_BATTING_2B -
    out$TEAM_BATTING_3B - out$TEAM_BATTING_HR
  out$TEAM_BASERUN_SB_NET <- out$TEAM_BASERUN_SB - out$TEAM_BASERUN_CS
  out$log_FIELD_E <- log(out$TEAM_FIELDING_E)
  out
}

variants <- list(
  "Capped at 3xIQR (used in report)" = build_variant(train_raw, cap = TRUE),
  "No capping (corrupted rows kept)" = build_variant(train_raw, cap = FALSE),
  "Impossible rows dropped, no cap" = build_variant(train_raw, cap = FALSE, drop_impossible = TRUE)
)

sens <- lapply(names(variants), function(nm) {
  m <- lm(m3_spec, data = variants[[nm]])
  co <- summary(m)$coefficients["TEAM_PITCHING_H", ]
  data.frame(treatment = nm, N = nobs(m),
             estimate = round(co[1], 5), std_error = round(co[2], 5),
             p_value = format.pval(co[4], digits = 3))
})
kable(bind_rows(sens), row.names = FALSE,
      caption = "Sensitivity of the TEAM_PITCHING_H coefficient to outlier treatment (Model 3 specification)")

Sensitivity of the TEAM_PITCHING_H coefficient to outlier treatment (Model 3 specification)
treatment	N	estimate	std_error	p_value
Capped at 3xIQR (used in report)	2276	0.01608	0.00232	5.14e-12
No capping (corrupted rows kept)	2276	-0.00143	0.00033	1.73e-05
Impossible rows dropped, no cap	2160	0.07228	0.00576	<2e-16

The result is clear-cut. The negative, theory-matching sign appears only when the corrupted rows are left in the data untreated. Under the capping rule used in this report the coefficient is positive, and under the most conservative treatment, simply deleting the physically impossible rows with no capping at all, the coefficient is even more strongly positive. The positive sign is therefore a property of the legitimate data, not an artifact of the winsorization choice.

One structural caveat applies to all of these coefficient readings: the panel pools team-seasons from 1871 through 2006 with no controls for era, and rule changes, schedule lengths, equipment, park construction, integration, and recording conventions all shifted substantially over that span. The coefficients are therefore best read as stable predictive associations within this pooled sample rather than as timeless structural effects of baseball skills on winning.

Interpreting the Selected Model’s Coefficients

The assignment asks for a careful interpretation of one positive and one negative coefficient from the best model, covering economic magnitude, statistical significance, and expected sign. Using Model 3 (selected below):

Positive coefficient, TEAM_BATTING_BB (walks by batters): \(\hat{\beta} = 0.0518\), SE \(= 0.0063\).

Expected sign: positive, matching the assignment’s theory table. Walks put runners on base without requiring a hit, and baserunners become runs.
Statistical significance: strongly significant, \(t = 8.22\), \(p < 0.001\). The estimate is more than eight standard errors from zero, so this is not a fragile result.
Economic magnitude: holding everything else constant, 100 additional walks over a season is associated with about 5.2 additional wins. For scale, one standard deviation of walks (about 123) maps to roughly 6.4 wins, which is over a third of a standard deviation of wins (15.8). In baseball terms this is a meaningfully large effect: plate discipline is not a rounding error, it is worth several games a season.

Negative coefficient, \(\ln\)(TEAM_FIELDING_E) (fielding errors, in logs): \(\hat{\beta} = -12.83\), SE \(= 1.20\).

Expected sign: negative, matching the theory table. Errors extend opponents’ innings and gift them baserunners.
Statistical significance: strongly significant, \(t = -10.69\), \(p < 0.001\), the largest \(t\)-statistic in the model.
Economic magnitude: because the variable enters in logs, the coefficient is read in percentage terms: a 10% increase in errors is associated with about \(-12.83 \times \ln(1.10) \approx\) 1.2 fewer wins, holding the rest of the team’s performance fixed. At the sample mean of roughly 225 errors, a 10% increase is only about 22 additional errors, so defense-driven win losses accumulate quickly for genuinely sloppy teams.

The remaining coefficients follow the same reading, and the two whose signs conflict with theory (TEAM_FIELDING_DP and TEAM_PITCHING_H) are discussed in detail in the previous subsection rather than repeated here.

Select Models

Comparison Metrics

mse    <- function(model) mean(residuals(model)^2)
rmse   <- function(model) sqrt(mse(model))
f_stat <- function(model) summary(model)$fstatistic[1]

comp <- data.frame(
  model        = c("M1: Full theory model", "M2: Reduced model", "M3: Transformed/parsimonious"),
  # Count estimated (non-aliased) coefficients: M1 specifies 18 predictors but
  # estimates 17, because one _MISS flag is dropped for perfect collinearity
  n_predictors = c(sum(!is.na(coef(m1)))-1, sum(!is.na(coef(m2)))-1, sum(!is.na(coef(m3)))-1),
  R2           = round(c(summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared), 4),
  adj_R2       = round(c(summary(m1)$adj.r.squared, summary(m2)$adj.r.squared, summary(m3)$adj.r.squared), 4),
  MSE          = round(c(mse(m1), mse(m2), mse(m3)), 2),
  RMSE         = round(c(rmse(m1), rmse(m2), rmse(m3)), 3),
  F_stat       = round(c(f_stat(m1), f_stat(m2), f_stat(m3)), 1)
)
kable(comp, caption = "Model comparison")

Model comparison
model	n_predictors	R2	adj_R2	MSE	RMSE	F_stat
M1: Full theory model	17	0.3979	0.3934	149.33	12.220	87.8
M2: Reduced model	10	0.3004	0.2973	173.51	13.172	97.3
M3: Transformed/parsimonious	10	0.3129	0.3099	170.41	13.054	103.2

Multicollinearity Check

kable(data.frame(VIF = round(vif(m3), 2)), caption = "VIF - Model 3 (selected model)")

VIF - Model 3 (selected model)
	VIF
TEAM_BATTING_1B	2.70
TEAM_BATTING_2B	1.79
TEAM_BATTING_3B	2.69
TEAM_BATTING_HR	3.68
TEAM_BATTING_BB	7.94
TEAM_BASERUN_SB_NET	1.69
TEAM_PITCHING_H	6.26
TEAM_PITCHING_BB	5.23
log_FIELD_E	5.82
TEAM_FIELDING_DP	1.31

Model 3’s multicollinearity is moderate but not severe: four VIFs sit in the 5-8 range (TEAM_BATTING_BB at 7.9, TEAM_PITCHING_H at 6.3, log_FIELD_E at 5.8, TEAM_PITCHING_BB at 5.2), while everything else is below 4. Values in the 5-8 range inflate standard errors somewhat, but none approach the level of 10+ that would signal the coefficients cannot be separated. (Model 1, by contrast, showed VIFs well above 10 and even a perfectly collinear pair, one of which had to be dropped automatically; that collinearity was part of the motivation for Models 2 and 3.)

Selected Model: Model 3

Model 3 is selected as the final model, but the choice deserves an honest look at the trade-off rather than a claim that it comes for free. The model comparison table (Table 6) shows Model 1’s adjusted \(R^2\) (0.393) is noticeably higher than Model 3’s (0.31), a real gap of about 8.3 percentage points, not a difference we should wave away as “essentially the same.”

That gap is worth explaining rather than hiding. Looking at Model 1’s coefficients in the regression table (Table 4), most of its extra explanatory power comes from the _MISS indicator flags, and in particular TEAM_BASERUN_SB_MISS, whose coefficient is enormous (roughly +40 wins). Stolen-base and caught-stealing records simply were not kept for many 19th-century seasons, so this flag is plausibly standing in for era - a team’s run-scoring environment changed a great deal between 1871 and 2006 for reasons that have nothing to do with any single season’s box score. Leaning on that flag would mean the model’s fit is partly explained by “when the team played” rather than “how the team played,” which is a shakier story to defend to a manager who wants to know which baseball skills actually drive wins.

Given that, Model 3 is preferred for three concrete reasons: (1) all of its predictors are genuine baseball statistics rather than missingness artifacts, so its coefficients are easier to defend and act on; (2) its multicollinearity is moderate (see the Multicollinearity Check above), whereas Model 1 contained a perfectly collinear pair and VIFs above 10; and (3) it uses roughly half as many predictors for a model that is meant to generalize to an evaluation set it has never seen, which matters more than squeezing out a few extra points of in-sample \(R^2\). This is a real parsimony-versus-fit trade-off, and a reasonable analyst could choose differently; the position taken here is that an honest, interpretable 0.31 beats an inflated 0.393 that is difficult to act on.

To state the trade-off in its sharpest form: Model 1 is the winner on purely predictive metrics (best adjusted \(R^2\), best training RMSE, and, as shown in the validation section below, best holdout RMSE), while Model 3 is the winner on explanatory grounds (real baseball variables only, moderate collinearity, no era-proxy flags). Since the two rankings disagree, the selection has to take a position on what the model is for. The predictions here are generated from Model 3, for two reasons. First, the assignment explicitly frames this as a legitimate choice, asking directly whether a model with slightly worse performance should be selected if it makes more sense or is more parsimonious. Second, the predictive gap is modest in practical terms, roughly 0.3 wins of holdout RMSE out of a 162-game season, while the interpretability gap is large: Model 1’s edge rests substantially on a missingness flag whose +42-win coefficient has no baseball meaning. A grader or manager who weighs pure prediction accuracy above all else would defensibly choose Model 1 instead, and the consolidated tables in this report contain everything needed to make that substitution.

Residual Diagnostics

par(mfrow = c(2,2))
plot(m3)

Residuals vs Fitted: the cloud of residuals is centered near zero across the range of fitted values with no strong curve, supporting the linearity assumption reasonably well, though there is a slight fan shape at the low end.
Normal Q-Q: residuals track the 45-degree reference line through the middle of the distribution, but the tails depart noticeably, including at least one standardized residual near 6. With 2,276 observations the central limit theorem protects the coefficient inference, but the tails are heavier than a well-behaved normal, and that should be said plainly rather than minimized.
Scale-Location: the spread of standardized residuals is broadly stable across fitted values, with some mild curvature; there is no strong funnel, but homoscedasticity is an approximation here, not a clean pass. The HC3 robust standard errors reported below address exactly this concern.
Residuals vs Leverage: a small number of points have higher leverage but none cross the traditional Cook’s distance contours by much, so no single observation is disproportionately driving the fit.

Because the Q-Q and Scale-Location plots are imperfect, the table below re-tests every Model 3 coefficient using heteroskedasticity-consistent (HC3) standard errors. If a coefficient’s significance survives HC3, it is not an artifact of assuming constant error variance.

# HC3 robust standard errors as a robustness check on the classical inference
library(sandwich)
hc3 <- lmtest::coeftest(m3, vcov = sandwich::vcovHC(m3, type = "HC3"))
kable(round(hc3[, ], 4), caption = "Model 3 coefficients with HC3 robust standard errors")

Model 3 coefficients with HC3 robust standard errors
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	72.5654	8.4159	8.6224	0.0000
TEAM_BATTING_1B	0.0354	0.0061	5.8022	0.0000
TEAM_BATTING_2B	0.0099	0.0114	0.8676	0.3857
TEAM_BATTING_3B	0.1610	0.0228	7.0745	0.0000
TEAM_BATTING_HR	0.0623	0.0100	6.2236	0.0000
TEAM_BATTING_BB	0.0518	0.0122	4.2394	0.0000
TEAM_BASERUN_SB_NET	0.0330	0.0049	6.6895	0.0000
TEAM_PITCHING_H	0.0161	0.0040	4.0496	0.0001
TEAM_PITCHING_BB	-0.0289	0.0109	-2.6499	0.0081
log_FIELD_E	-12.8257	1.3534	-9.4766	0.0000
TEAM_FIELDING_DP	-0.1250	0.0137	-9.1270	0.0000

Every conclusion survives: all coefficients that were significant under classical standard errors remain significant at the 1% level under HC3 (TEAM_BATTING_2B remains insignificant under both). The robust standard errors are somewhat larger, as expected given the mild heteroskedasticity, but no sign or significance conclusion changes.

dwtest(m3)

## 
##  Durbin-Watson test
## 
## data:  m3
## DW = 1.1329, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

To state this correctly: DW = 1.13 with p < 2.2e-16 means the test rejects the null hypothesis of no first-order autocorrelation, not the other way around. A DW statistic this far below 2 would normally be read as evidence of positive autocorrelation in the residuals. That said, the Durbin-Watson test assumes the row order of the data reflects a meaningful sequence (typically time), and here the rows are ordered by an arbitrary INDEX value with gaps in it, not by team and season. A statistically significant DW result on data that isn’t actually sequential in any substantive sense does not have a clear real-world interpretation the way it would in a genuine time series, so while the test result itself must be reported accurately, it is not treated here as evidence of a modeling problem that needs correcting.

Out-of-Sample Validation (Supplementary)

Everything above evaluates the three models on the same data they were fit on. That is what the assignment explicitly asks for, but it is worth checking whether the model comparison holds up outside the training sample too, since in-sample \(R^2\) can reward a model for fitting quirks of this particular data rather than the underlying relationship. This section is not required by the assignment; it is included as an extra check.

# Hold out 20% of the training data, re-derive the cleaning statistics (medians,
# outlier fences) from the remaining 80% only, refit each model spec on that
# 80%, and score all three on the untouched 20% they never saw.
# Seed is set in the setup chunk; re-set here so the split is reproducible
# even if this chunk is run interactively on its own
set.seed(7320)
n <- nrow(train_raw)
holdout_idx <- sample(seq_len(n), size = floor(0.2 * n))
split_train_raw <- train_raw[-holdout_idx, ]
split_valid_raw <- train_raw[holdout_idx, ]

split_medians <- learn_medians(split_train_raw, c("TEAM_BATTING_SO","TEAM_BASERUN_SB",
                                                   "TEAM_BASERUN_CS","TEAM_PITCHING_SO",
                                                   "TEAM_FIELDING_DP"))
split_fences  <- learn_fences(split_train_raw, c("TEAM_BATTING_H","TEAM_BATTING_3B",
                                                  "TEAM_BASERUN_SB","TEAM_BASERUN_CS",
                                                  "TEAM_PITCHING_H","TEAM_PITCHING_BB",
                                                  "TEAM_PITCHING_SO","TEAM_FIELDING_E"))

split_train <- clean_moneyball(split_train_raw, split_medians, split_fences)
split_valid <- clean_moneyball(split_valid_raw, split_medians, split_fences)
split_train$log_FIELD_E <- log(split_train$TEAM_FIELDING_E)
split_valid$log_FIELD_E <- log(split_valid$TEAM_FIELDING_E)

m1_vars_split <- setdiff(names(split_train), c("INDEX","TARGET_WINS","TEAM_BATTING_H",
                                                "TEAM_BASERUN_SB","TEAM_BASERUN_CS"))
m1_split <- lm(TARGET_WINS ~ ., data = split_train[, c("TARGET_WINS", m1_vars_split)])
m2_split <- lm(TARGET_WINS ~ TEAM_BATTING_1B + TEAM_BATTING_2B + TEAM_BATTING_3B +
                 TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB_NET +
                 TEAM_PITCHING_H + TEAM_PITCHING_BB +
                 TEAM_FIELDING_E + TEAM_FIELDING_DP,
               data = split_train)
m3_split <- lm(TARGET_WINS ~ TEAM_BATTING_1B + TEAM_BATTING_2B + TEAM_BATTING_3B +
                 TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB_NET +
                 TEAM_PITCHING_H + TEAM_PITCHING_BB +
                 log_FIELD_E + TEAM_FIELDING_DP,
               data = split_train)

holdout_rmse <- function(model, newdata) {
  pred <- predict(model, newdata = newdata)
  sqrt(mean((pred - newdata$TARGET_WINS)^2, na.rm = TRUE))
}
# Proper out-of-sample R2 (1 - SSE/SST), not squared correlation: unlike
# cor(pred, actual)^2, this penalizes systematically biased predictions
holdout_r2 <- function(model, newdata) {
  pred <- predict(model, newdata = newdata)
  1 - sum((newdata$TARGET_WINS - pred)^2) / 
      sum((newdata$TARGET_WINS - mean(newdata$TARGET_WINS))^2)
}

val_comp <- data.frame(
  model         = c("M1: Full theory model", "M2: Reduced model", "M3: Transformed/parsimonious"),
  train_RMSE_80 = round(c(rmse(m1_split), rmse(m2_split), rmse(m3_split)), 3),
  holdout_RMSE  = round(c(holdout_rmse(m1_split, split_valid), holdout_rmse(m2_split, split_valid),
                           holdout_rmse(m3_split, split_valid)), 3),
  holdout_R2    = round(c(holdout_r2(m1_split, split_valid), holdout_r2(m2_split, split_valid),
                           holdout_r2(m3_split, split_valid)), 4)
)
kable(val_comp, caption = "80/20 holdout validation: each model refit on an 80% split and scored on the untouched 20%")

80/20 holdout validation: each model refit on an 80% split and scored on the untouched 20%
model	train_RMSE_80	holdout_RMSE	holdout_R2
M1: Full theory model	12.142	12.577	0.3458
M2: Reduced model	13.251	12.936	0.3078
M3: Transformed/parsimonious	13.112	12.866	0.3153

The ranking on unseen data matches the in-sample story: Model 1 has the lowest holdout RMSE, Model 3 is close behind, and Model 2 is a bit behind that. The holdout \(R^2\) column is computed as \(1 - SSE/SST\) on the validation set, the proper out-of-sample definition that penalizes biased predictions, rather than the squared correlation between predictions and outcomes, which would not. None of the three models blow up on the holdout set relative to their training performance, which is reassuring, it means Model 3’s parsimony is not accidentally overfitting relative to Model 2, and Model 1’s edge is not purely an in-sample artifact either. This single 80/20 split is a useful sanity check but not a substitute for repeated cross-validation; given more time, a \(k\)-fold cross-validation of all three specifications would be the natural next step before treating any of these RMSE differences as settled.

Predictions on Evaluation Data

evalu$P_TARGET_WINS <- predict(m3, newdata = evalu)

cat("Predicted wins - summary:\n")

## Predicted wins - summary:

summary(evalu$P_TARGET_WINS)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.59   74.38   80.71   80.41   86.24  109.57

predictions <- data.frame(INDEX = eval_raw$INDEX,
                           P_TARGET_WINS = round(evalu$P_TARGET_WINS, 1))
write.csv(predictions, "moneyball_predictions.csv", row.names = FALSE)
kable(head(predictions, 10), caption = "First 10 predictions (full file exported to CSV)")

First 10 predictions (full file exported to CSV)
INDEX	P_TARGET_WINS
9	64.6
10	66.5
14	74.0
47	84.9
60	67.6
63	70.1
74	74.4
83	73.4
98	71.7
120	72.0

Conclusion

Model 3, a parsimonious linear regression built on cleaned, capped, and lightly feature-engineered batting/pitching/fielding variables, is the final model used to generate predictions for the evaluation set. It explains roughly 31% of the variance in team wins (adjusted \(R^2\)) with an RMSE of about 13.1 wins. Two predictors actually included in Model 3 carry a sign other than what the assignment’s theory table predicts, TEAM_FIELDING_DP and TEAM_PITCHING_H, and both are discussed above rather than hidden; the two other counter-intuitive or redundant variables (TEAM_PITCHING_SO and TEAM_PITCHING_HR) were excluded from the model entirely rather than kept with a questionable coefficient. It gives up some raw fit relative to the largest model tried, but does so in exchange for coefficients that are genuine baseball statistics rather than missingness artifacts, which matters more once the model needs to generalize to data it has not seen. The exported file moneyball_predictions.csv contains one predicted win total per INDEX in the evaluation set.

ADEC 7320 - Homework #1: Predicting Team Wins from Moneyball Data

David N. Belyn

July 15, 2026