In-Class Activity #15 — Predicting the World Series Winner

##Read in and explore the data

baseball = read.csv("baseball.csv")
str(baseball)

'data.frame':   1232 obs. of  15 variables:
 $ Team        : chr  "ARI" "ATL" "BAL" "BOS" ...
 $ League      : chr  "NL" "NL" "AL" "AL" ...
 $ Year        : int  2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ...
 $ RS          : int  734 700 712 734 613 748 669 667 758 726 ...
 $ RA          : int  688 600 705 806 759 676 588 845 890 670 ...
 $ W           : int  81 94 93 69 61 85 97 68 64 88 ...
 $ OBP         : num  0.328 0.32 0.311 0.315 0.302 0.318 0.315 0.324 0.33 0.335 ...
 $ SLG         : num  0.418 0.389 0.417 0.415 0.378 0.422 0.411 0.381 0.436 0.422 ...
 $ BA          : num  0.259 0.247 0.247 0.26 0.24 0.255 0.251 0.251 0.274 0.268 ...
 $ Playoffs    : int  0 1 1 0 0 0 1 0 0 1 ...
 $ RankSeason  : int  NA 4 5 NA NA NA 2 NA NA 6 ...
 $ RankPlayoffs: int  NA 5 4 NA NA NA 4 NA NA 2 ...
 $ G           : int  162 162 162 162 162 162 162 162 162 162 ...
 $ OOBP        : num  0.317 0.306 0.315 0.331 0.335 0.319 0.305 0.336 0.357 0.314 ...
 $ OSLG        : num  0.415 0.378 0.403 0.428 0.424 0.405 0.39 0.43 0.47 0.402 ...

# How many team/year pairs in the whole dataset?
nrow(baseball)                    # 1232

[1] 1232

# How many distinct years?
length(table(baseball$Year))      # 47

[1] 47

Limit to teams that made the playoffs

baseball = subset(baseball, Playoffs == 1)
nrow(baseball)                    # 244

[1] 244

table(table(baseball$Year))

 2  4  8 10 
 7 23 16  1 

Add the NumCompetitors predictor

PlayoffTable = table(baseball$Year)
PlayoffTable[c("1990", "2001")]

1990 2001 
   4    8 

baseball$NumCompetitors = PlayoffTable[as.character(baseball$Year)]
table(baseball$NumCompetitors)    # 8-team years -> 128

  2   4   8  10 
 14  92 128  10 

Define the WorldSeries outcome

baseball$WorldSeries = as.numeric(baseball$RankPlayoffs == 1)
table(baseball$WorldSeries)       # 197 did NOT win, 47 won

  0   1 
197  47 

Twelve bivariate logistic regression models

model1  <- glm(WorldSeries ~ Year,           data=baseball, family="binomial")
model2  <- glm(WorldSeries ~ RS,             data=baseball, family="binomial")
model3  <- glm(WorldSeries ~ RA,             data=baseball, family="binomial")
model4  <- glm(WorldSeries ~ W,              data=baseball, family="binomial")
model5  <- glm(WorldSeries ~ OBP,            data=baseball, family="binomial")
model6  <- glm(WorldSeries ~ SLG,            data=baseball, family="binomial")
model7  <- glm(WorldSeries ~ BA,             data=baseball, family="binomial")
model8  <- glm(WorldSeries ~ RankSeason,     data=baseball, family="binomial")
model9  <- glm(WorldSeries ~ OOBP,           data=baseball, family="binomial")
model10 <- glm(WorldSeries ~ OSLG,           data=baseball, family="binomial")
model11 <- glm(WorldSeries ~ NumCompetitors, data=baseball, family="binomial")
model12 <- glm(WorldSeries ~ League,         data=baseball, family="binomial")

summary(model1)   # Year            -> significant

Call:
glm(formula = WorldSeries ~ Year, family = "binomial", data = baseball)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) 72.23602   22.64409    3.19  0.00142 **
Year        -0.03700    0.01138   -3.25  0.00115 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 239.12  on 243  degrees of freedom
Residual deviance: 228.35  on 242  degrees of freedom
AIC: 232.35

Number of Fisher Scoring iterations: 4

summary(model3)   # RA              -> significant

Call:
glm(formula = WorldSeries ~ RA, family = "binomial", data = baseball)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)  1.888174   1.483831   1.272   0.2032  
RA          -0.005053   0.002273  -2.223   0.0262 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 239.12  on 243  degrees of freedom
Residual deviance: 233.88  on 242  degrees of freedom
AIC: 237.88

Number of Fisher Scoring iterations: 4

summary(model8)   # RankSeason      -> significant

Call:
glm(formula = WorldSeries ~ RankSeason, family = "binomial", 
    data = baseball)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.8256     0.3268  -2.527   0.0115 *
RankSeason   -0.2069     0.1027  -2.016   0.0438 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 239.12  on 243  degrees of freedom
Residual deviance: 234.75  on 242  degrees of freedom
AIC: 238.75

Number of Fisher Scoring iterations: 4

summary(model11)  # NumCompetitors  -> significant (strongest)

Call:
glm(formula = WorldSeries ~ NumCompetitors, family = "binomial", 
    data = baseball)

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     0.03868    0.43750   0.088 0.929559    
NumCompetitors -0.25220    0.07422  -3.398 0.000678 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 239.12  on 243  degrees of freedom
Residual deviance: 226.96  on 242  degrees of freedom
AIC: 230.96

Number of Fisher Scoring iterations: 4

# The remaining eight (RS, W, OBP, SLG, BA, OOBP, OSLG, League) are NOT significant.

Significant bivariate predictors: Year, RA, RankSeason, NumCompetitors.

Multivariate model

LogModel = glm(WorldSeries ~ Year + RA + RankSeason + NumCompetitors,
               data=baseball, family=binomial)
summary(LogModel)   # None of the four is significant -> multicollinearity

Call:
glm(formula = WorldSeries ~ Year + RA + RankSeason + NumCompetitors, 
    family = binomial, data = baseball)

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)    12.5874376 53.6474210   0.235    0.814
Year           -0.0061425  0.0274665  -0.224    0.823
RA             -0.0008238  0.0027391  -0.301    0.764
RankSeason     -0.0685046  0.1203459  -0.569    0.569
NumCompetitors -0.1794264  0.1815933  -0.988    0.323

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 239.12  on 243  degrees of freedom
Residual deviance: 226.37  on 239  degrees of freedom
AIC: 236.37

Number of Fisher Scoring iterations: 4

Correlations among the four predictors

cor(baseball[c("Year", "RA", "RankSeason", "NumCompetitors")])

                    Year        RA RankSeason NumCompetitors
Year           1.0000000 0.4762422  0.3852191      0.9139548
RA             0.4762422 1.0000000  0.3991413      0.5136769
RankSeason     0.3852191 0.3991413  1.0000000      0.4247393
NumCompetitors 0.9139548 0.5136769  0.4247393      1.0000000

# Only Year / NumCompetitors are strongly correlated (0.914)

Two-variable models and AIC

model13 = glm(WorldSeries ~ Year + RA,                   data=baseball, family=binomial)
model14 = glm(WorldSeries ~ Year + RankSeason,           data=baseball, family=binomial)
model15 = glm(WorldSeries ~ Year + NumCompetitors,       data=baseball, family=binomial)
model16 = glm(WorldSeries ~ RA + RankSeason,             data=baseball, family=binomial)
model17 = glm(WorldSeries ~ RA + NumCompetitors,         data=baseball, family=binomial)
model18 = glm(WorldSeries ~ RankSeason + NumCompetitors, data=baseball, family=binomial)

# Compare AIC across all 10 models; lowest wins
AIC(model1, model3, model8, model11,
    model13, model14, model15, model16, model17, model18)

The lowest AIC belongs to model11 (NumCompetitors alone), AIC = 230.96.

Conclusion

None of the two-variable models had both variables significant, so none beats a simple bivariate model. The best model overall uses only NumCompetitors — the number of playoff teams — which has nothing to do with team quality. Once a team makes the playoffs, regular-season performance does not predict the World Series winner, supporting Billy Beane’s claim in Moneyball that winning the playoffs comes down to luck.