World Series

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 ...

Basic checks

nrow(baseball)

[1] 1232

length(table(baseball$Year))

[1] 47

There are 1232 team/year rows and 47 years.

Playoff teams

baseball <- subset(baseball, Playoffs == 1)
nrow(baseball)

[1] 244

table(baseball$Year)

1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1973 1974 1975 1976 1977 1978 1979 1980 1982 1983 
   2    2    2    2    2    2    2    4    4    4    4    4    4    4    4    4    4    4    4    4 
1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 
   4    4    4    4    4    4    4    4    4    4    8    8    8    8    8    8    8    8    8    8 
2006 2007 2008 2009 2010 2011 2012 
   8    8    8    8    8    8   10 

There are 244 playoff team/year rows.

playoff_table <- table(baseball$Year)
baseball$NumCompetitors <- as.numeric(playoff_table[as.character(baseball$Year)])
table(baseball$NumCompetitors)

  2   4   8  10 
 14  92 128  10 

There are 128 rows from seasons with 8 playoff teams.

World Series variable

baseball$WorldSeries <- as.numeric(baseball$RankPlayoffs == 1)
table(baseball$WorldSeries)

  0   1 
197  47 

197 teams did not win and 47 teams won.

Bivariate logistic models

predictors <- c("Year", "RS", "RA", "W", "OBP", "SLG", "BA",
                "RankSeason", "OOBP", "OSLG", "NumCompetitors", "League")

models <- lapply(predictors, function(x) {
  glm(reformulate(x, "WorldSeries"), data = baseball, family = binomial)
})
names(models) <- predictors

p_values <- sapply(models, function(m) coef(summary(m))[2, 4])
p_values

          Year             RS             RA              W            OBP            SLG 
  0.0011544695   0.2013364471   0.0262137758   0.0577260348   0.2959201177   0.0504256593 
            BA     RankSeason           OOBP           OSLG NumCompetitors         League 
  0.8385900817   0.0438313556   0.9015582281   0.7571369452   0.0006784323   0.6264465507 

names(p_values[p_values < 0.05])

[1] "Year"           "RA"             "RankSeason"     "NumCompetitors"

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

Combined model

full_model <- glm(
  WorldSeries ~ Year + RA + RankSeason + NumCompetitors,
  data = baseball,
  family = binomial
)
summary(full_model)

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

None of the predictors are significant in the combined model.

Correlation

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

Year and NumCompetitors have the strongest correlation, about 0.914.

Compare models by AIC

formulas <- list(
  Year = WorldSeries ~ Year,
  RA = WorldSeries ~ RA,
  RankSeason = WorldSeries ~ RankSeason,
  NumCompetitors = WorldSeries ~ NumCompetitors,
  Year_RA = WorldSeries ~ Year + RA,
  Year_RankSeason = WorldSeries ~ Year + RankSeason,
  Year_NumCompetitors = WorldSeries ~ Year + NumCompetitors,
  RA_RankSeason = WorldSeries ~ RA + RankSeason,
  RA_NumCompetitors = WorldSeries ~ RA + NumCompetitors,
  RankSeason_NumCompetitors = WorldSeries ~ RankSeason + NumCompetitors
)

comparison <- sapply(formulas, function(f) {
  AIC(glm(f, data = baseball, family = binomial))
})
sort(comparison)

           NumCompetitors                      Year RankSeason_NumCompetitors 
                 230.9592                  232.3508                  232.5240 
        RA_NumCompetitors       Year_NumCompetitors           Year_RankSeason 
                 232.7429                  232.8972                  233.5524 
                  Year_RA                        RA             RA_RankSeason 
                 233.8766                  237.8839                  238.2180 
               RankSeason 
                 238.7546 

The NumCompetitors-only model has the lowest AIC.