Beyond runs expectancy

1.1George Lindsey and the runs expectancy matrix

One the pioneers in sabermetics was George Lindsey, a defense consultant in Canada who had a great love for baseball. Lindsey wrote two remarkable papers in the 1960’s that had a great influence on the quantitative analysis of baseball. In particular, Lindsey (1963) focused on several questions of baseball strategy such as stealing a base, sacrificing to sacrifice a runner to second base, and issuing an intentional walk. Lindsay observed that these questions could be answered by the collection of appropriate data.

Lindsay noted that the probability of winning at a point during the game depends on the current score and inning and presented tables of W (I, H_i), the probability a team with a lead of I runs at the home half of the ith inning H_i will win the game.

In this paper, Lindsey also considered the runs scoring distribution of a team during an inning. He focused on the run potential, the number of runs a team will score in the remainder of the inning. He noted that the run potential depends on the current number of outs and the runners on base. Since there are three possible number of outs and eight possible runner configurations, there are 3 × 8 = 24 possible out/runner situations, and Lindsey focused on computing

To learn about run potential, Lindsey’s father collected run-scoring data for over 6000 half-innings in games during the 1959 and 1960 seasons, and Lindsey was able to find empirical run-scoring distributions for each of the 24 situations. By computing the mean runs for each situation, he produced the runs expectancy matrix as displayed in Table 1.

At the beginning of an inning with no outs and bases empty, a team will score, on average, 0.46 runs in the remainder of the inning. In contrast, when the bases are loaded with one out, a team will score 1.64 runs in the rest of the inning.

The runs expectancy table has been illustrated in a number of sabermetrics books such as Thorn et al. (1984), Albert and Bennett (2003), Keri et al. (2007) and Tango et al. (2007) for deciding on proper baseball strategy, for learning about the value of plays, and for evaluating players. For example, Chapter 7 of Albert and Bennett (2003) shows the use of run expectancies in determining the values of singles, doubles, triples, and home runs that leads to use of linear weights to measure batting performances. Chapter 9 of Albert and Bennett (2003) uses the run expectancy table to assess the value of common baseball strategies such as a sacrifice bunt, intentional walk, and stealing a base.

Lindsay (1963) cautioned that the run expectancy table represented the situation where all players were “average”.

Also Lindsey (1963) mentioned the possibility of extending this analysis to run-scoring of teams with non-average players.

1.2Modeling runs scored in an inning

A related line of research is the use of probability models to represent the number of runs scored in an inning or a game. Fits of basic discrete probability models do not accurately represent actual runs scored in an inning. To demonstrate this claim, the Poisson (λ) and negative binomial (n, p) models were separately fit to all inning runs scored in the 2013 season and the fitted probabilities for these two models and the actual run-scoring distribution are displayed in Table 2. Looking at the table, the Poisson fit underestimates the probability of a scoreless inning. Generally, the Poisson fit understates the actual variability in runs scored. The negative binomial fit is better than the Poisson, but this model overestimates the probability of scoring one run and underestimates the probability of scoring two runs.

Since basic probability models appear inaccurate, more sophisticated models for run scoring have been proposed in the literature. For example, Rosner et al. (1996) represents the probability of scoring x runs in an inning by the sum

1.3Multilevel modeling

Another research theme in modeling of sports data is the use of multilevel or hierarchical models to estimate parameters from several groups. Efron and Morris (1975) illustrate the use of an exchangeable model to estimate true batting averages of 18 players from the 1970 season. Albert (2004, 2006) provide further illustrations of the use of an exchangeable model to estimate a set of true batting rates or pitching rates. These models can be extended to the regression framework. Morris (1983) uses an multilevel regression model to assess if Ty Cobb was ever a true.400 hitter. Albert (2009) uses an exchangeable regression model to simultaneously estimate the true pitching trajectories for a number of pitchers.

1.4Plan of the paper

This paper provides an multilevel modeling framework for understanding run scoring of teams in Major League Baseball, and understanding how team run scoring is affected by several factors, such as ballpark, pitcher quality, and runners in scoring position. The basic multilevel model, described in Section 2.1, simultaneously estimates means from several populations, say the average number of runs scored for 30 MLB teams, when the means are believed to follow a normal curve model with unknown mean and standard deviation.

As in Lindsay’s work, we focus on the number of runs scored in a half-inning. A run scoring distribution for a particular team is represented by a multinomial distribution with underlying probabilities over the classes 0 runs, 1 run, 2 runs, 3 runs, and 4 or more runs. Section 2.2 describes a multilevel model for simultaneously estimating a set of multinomial probability vectors. By using this model to estimate the run scoring distributions for the 30 MLB teams, one obtains “improved” estimates at the teams’ abilities to score runs.

After estimating teams’ scoring distributions, we next explore team situational effects. An ordinal regression model, described in Section 2.3, is a useful method for describing how run-scoring probabilities change as one changes a predictor such the quality of a pitcher. This model can be described in terms of the odds of scoring at least a particular number of runs. It is called a proportional odds model as a change in the predictor will result in the odds of scoring runs being increased by a constant factor. Section 3.2 introduces the use of odds in comparing run scoring distributions of the National and American Leagues. Section 3.3 extends this approach to compare run scoring of the 30 teams and illustrates the value of multilevel modeling.

Section 4 focuses on the effect of the following covariates on team run scoring.

2.1Estimating a group of means

Consider the problem of estimating a collection of means from several populations. One observes sample means y₁, . . . , y_N, where the sample mean from the jth population, y_j is distributed normal with mean θ_j and known standard error σ_j. We describe this situation in a baseball setting, where N = 30, y₁, . . . , y₃₀ correspond to the sample mean numbers of runs scored in a half-inning for the thirty teams for a season, and θ₁, . . . , θ_N are the corresponding population means of runs scored.

One can obtain obtained improved estimates at the population means by means of the following multilevel model. We assume that the means θ₁, . . . , θ₃₀ represent a sample from a normal curve with mean μ and standard deviation τ. The locations of the normal curve parameters μ and τ are unknown, and so (under a Bayesian framework) vague or imprecise prior distributions are assigned to these parameters.

One can fit this multilevel model from the observed sample means of runs scored for the 30 teams. In particular, we estimate the population mean and standard deviation μ and τ by μˆ and τˆ, respectively. The estimate μˆ represents the overall or combined estimate of a population mean of the runs scored from the 30 teams and the estimate τˆ represents the spread of these population means.

Three sets of estimates

There are three types of estimates of the population means. If we estimate each population mean only by using data from the corresponding sample, we obtain individual estimates

By use of the multilevel model, we obtain improved estimates of the population means that compromise between the individual and combined estimates. The multilevel (ML) estimate of the population mean for team j, θˆjML, is given by

2.2Estimating collections of multinomial data

A related problem is estimating a multinomial population. Suppose we classify data into k bins and observe the vector of frequencies W = (w₁, . . . , w_k). The vector W is assumed to have a multinomial distribution with sample size N = w₁ + . . . + w_k and probability vector p = (p₁, . . . , p_k), where p_j represents the probability that an observation falls in the jth bin. In our setting, we classify the number of runs scored in an inning in the five bins 0, 1, 2, 3, and 4 or more runs, and the frequencies w₁, . . . , w₅ represent the number of innings where a particular team scores the different number of runs.

This type of multinomial data on runs scored is collected for each of the 30 baseball teams. Let W ¹, . . . , W ³⁰ denote the vectors of frequencies of runs scored for the 30 teams. We assume the vector for the jth team W ^j is multinomial with probability vector p ^j . We are interested in estimating the probability vectors p ¹, . . . , p ³⁰. As in the population means case, we can consider individual, combined, and multilevel estimates for the probability vectors.

Individual estimates

Suppose we use the data from only the jth team to learn about its run scoring tendencies. Then the probabilities of falling in the different run groups are estimated by the individual estimates:

Combined estimates

Instead, suppose we assume that the probabilities of falling in the different runs group are the same for all teams; that is, p ¹ = . . . = p ³⁰. Then we can estimate each team’s probability vector by the combined estimate

Multilevel estimates

We wish to get improved estimates at the team probability vectors p ¹, . . . , p ³⁰ that compromise between the individual estimates and the combined estimates. This is accomplished by means of the following multilevel model. We assume that the unknown probability vectors p ¹, . . . , p ³⁰ are a random sample from a Dirichlet distribution with mean vector η and precision K with density proportional to

One fits this multilevel model to the observed count data and one obtains estimates at K and η from the posterior distribution – call these estimates Kˆ and ηˆ. (Details of this computation are provided in the appendix.) The estimate ηˆ is approximately the combined estimate pˆC. Then the multilevel estimate at the probability vector of team j is given by

2.3Groups of ordinal multinomial data with regression

Ordinal logistic regression

Suppose one observes the vector of run frequencies w = (w₁, . . . , w_k) which is multinomial with probability vector p = (p₁, . . . , p_k). The categories 1,..., k are ordinal – in our setting, these categories will be “0 runs”, “1 run”, “2 runs”, etc. Define the probability

Suppose one observes a variable x that influences the run-scoring probabilities p. The ordinal logistic model (see McCullagh (1980) and Johnson and Albert (1999)) says that the log of the odds of scoring at least c runs is a linear function of x. This model is written as

Ordinal logistic regression over groups

As a more general set-up, suppose we observe run-scoring frequencies for N = 30 teams, where the frequencies for the jth team w ^j is multinomial with probability vector pj=(p1j,...,pkj). With the same covariate x, we fit the ordinal logistic model to w ^j of the form

One can now use the estimation of separate means methodology of Section 2.1 to get improved (multilevel) estimates at the regression effects. The true regression effects β₁, . . . , β_N are given a normal distribution with mean μ and standard deviation τ, we place vague priors on μ and τ. Improved estimates are provided by fitting this multilevel model – the estimates have the general form

3.1Run scoring of all teams in 2013 season

We begin by collecting in Table 3 the runs scored for all complete innings (three outs) in the 2013 season. Since scoring five or more runs in an inning is relatively rare, Table 4 collapses the data from Table 3 by combining the counts of runs four or greater into the class “4+”.

3.2Logits and comparing scoring of the NL and AL

To motivate our general approach, a useful reexpression of a proportion p is the logit L=log(p1-p). For example, using data from Table 4, the logit of the proportion of innings where at least one run is scored is

Logits are useful in comparing groups of ordinal data such as runs scored. To illustrate, suppose we want to compare the runs scored per inning of National and American League teams in the 2013 season. First we obtain the counts of runs scored by the two leagues as displayed in Table 6. Each row of counts of the table is converted to proportions, and then logits are computed using each of the four breakpoints. The logits for each league are displayed in Table 7 and the “Difference” row gives the difference in the logits for the two leagues.

An ordinal regression model is a useful way to summarize the relationship seen in Table 7. If θ_c is the probability of scoring at least c runs, then the model is

To gain some historical perspective on the run-scoring advantage of the American League, this ordinal regression model was fit for each of the 16 seasons 1998 through 2013. Figure 1 displays a graph of the AL advantage (on the logit scale) against season. It is interesting to note that there is substantial variability in the AL advantage – for example the advantage was 0.121 in 1998 and decreased sharply to 0.019 in 2007. But there is no clear pattern in the changes in AL advantage over this time period. Generally, the probability an AL team scores a large number of runs is about 6% larger than the probability a NL team scores a large number of runs.

3.3Comparing scoring of teams

How does run scoring vary across teams? If we break down this runs table by the batting team in Table 8, we see some interesting variation in runs scored. If one looks at the percentage of big innings (four runs or more), Boston’s percentage of big innings is 3.2, contrasted with Miami’s percentage of 1.4. Boston’s percentage of scoreless innings is 68.9, while the New York Mets were scoreless in 75.8 percent of their innings.

The logit approach described in Section 3 is used to facilitate comparisons in team scoring. Logits were computed for each team using the four breakpoints. For example, for Philadelphia, we computed the four logits

3.4Multilevel modeling of team run distributions

Since there appears to be chance variability in the estimates of run scoring for the individual teams, there may be some advantage to estimating the run scoring distributions simultaneously using a multilevel model. The exchangeable model of Section 2.2 is fit to the run distributions for the 30 teams. The estimate of the smoothing parameter K is 1973. Most teams play close to the median number of innings 1457. So the size of the shrinkage of the individual estimates towards the combined estimate is

4.1General method

Section 3 introduced the use of logits in the comparison of run-scoring distributions. To explore the effect of specific covariates like home/away, pitcher quality, and clutch hitting on scoring, we fit the general ordinal logistic model of the form

4.2Home versus away effects

To consider a team’s scoring advantage of playing at home, write the logit of the probability of scoring at least c runs by

This ordinal regression model is initially fit to run-scoring data for all teams in the 2013 season and one obtains the covariate estimate βˆ=0.032. The interpretation is that the ratio of the odds of scoring at least c runs at home and away is equal to exp(0.032) =1.032. Specifically, the probability that a team scores at least one run in an inning is increased by 3% at home, the probability a team scores at least two runs in an inning is increased by 3% , and so on.

Since ballparks are known to have a significant impact on scoring, it is natural to fit this home/away model for each of the 30 teams. Figure 4 displays the individual estimates of the parameters β₁, . . . , β₃₀ as black dots where the endpoints of the bars correspond to the estimates plus and minus the standard errors. As expected, the individual home/away estimates show substantial variability from the New York Mets (βˆj=-0.31) to the Colorado Rockies (βˆj=0.52).

Next we apply the multilevel model of Section 2.1 to simultaneously estimate the 30 home effects β₁, . . . , β₃₀. One obtains the multilevel estimates μˆ=0.071 and τˆ=0.114. The estimate μˆ represents an average ballpark effect and τˆ is an estimate at the spread of the true ballpark effects. Here the multilevel model estimates of {β_j} shrink the individual estimates 52% of the way towards the combined estimate μˆ. The red bars in Fig. 4 show the posterior means plus and minus the posterior standard deviations. Since the size of the shrinkage is relatively small, this indicates that there is sizeable variation between the true home park effects.

4.3Pitcher effects

The same ordinal regression approach can to be used to learn about the pitcher effect on team run-scoring. The logit of the probability of scoring at least c runs in a half-inning is given by

The challenge in this approach is to find a good measure of pitcher quality. Many pitchers are of the relief type with a small number of plate appearances, and the measurements of pitcher quality typically show high variability. We use a multilevel model approach to obtain improved estimates of pitcher quality and these improved estimates are used as covariates in the ordinal model for run scoring.

To begin, we use a run value approach illustrated in Albert and Bennett (2003), Keri et al. (2007) and Tango et al. (2007) to measure the quality of all pitchers who played in the 2013 season. For each plate appearance (PA), we measure the run value as

Figure 5 displays a histogram of the run value measures of the 679 pitchers in the 2013 season. These multilevel estimates shrink the individual mean run value estimates towards the overall value, the degree of shrinkage depending on the number of batters faced. For pitchers who faced fewer than 100 batters, the shrinkage exceeds 80%; for a starter facing 900 batters, the shrinkage was only 35%. In this particular season, Clayton Kershaw stands out with a run value measure of −0.036 – since he faced 934 batters, he saved 934 × -0.036 = -33.6 runs during the 2013 season.

If the ordinal regression model with pitcher quality covariate is fit to data for all teams, one obtains the estimate βˆ=19.51. The standard deviation of all pitcher abilities is 0.012. So if the pitcher’s ability measure is increased by one standard deviation, one predicts a team’s log odds of scoring runs to increase by 19.51 × 0.012 = 0.23.

By fitting the multilevel ordinal regression model, we learn about the pitcher effects for all 30 teams. Figure 6 displays individual (team) effects by black bars and the multilevel model estimates are displayed by red bars. The individual estimates show substantial variation. For example, Seattle’s regression coefficient is βˆj=9.84 and Philadelphia’s coefficient is βˆj=28.24, indicating that the Phillies were much more able to take advantage of poor pitching than Seattle in the 2013 season. The multilevel estimates in this case shrink the individual estimates about 76% towards the average value. Seattle and Philadelphia’s pitcher effects, under the multilevel model, are corrected to 16.85 and 21.05, respectively. These multilevel estimates are better reflections of the true team pitcher effects than the individual team estimates, and would lead to better predictions of run-scoring against pitchers of different abilities.

4.4Advancing baserunners home (clutch effects)

Run scoring generally is viewed as a two-step process. A team places runners on bases and then these runners are advanced to home. There is a special focus in the media on advancing runners who are in scoring position (on second and third base). One commonly hears about the number of runners left on base and the proportion of runners in scoring position who score. Do teams really differ in their abilities to advance runners from scoring position?

To address this clutch hitting question, an ordinal regression model is constructed. For each half-inning, one records the total number of (unique) runners SP who are in scoring position (either on second or third base) and the number of these runners R who eventually score. (Note that R can be smaller than the number of runs that score in the half-inning since runners who score don’t need to be in scoring position.) As before, we classify R into the categories 0, 1, 2, 3, and “4 or more”, and consider the ordinal regression model

To see how this clutch hitting statistic varies between teams, this ordinal regression model was fit separately for all teams. Anaheim and the New York Mets had the smallest and largest regression estimates of βˆj=2.19 and βˆj=2.60, respectively. This indicates that the Mets were the best team in baseball in 2013 from a clutch-hitting perspective and Anaheim was the worst. However, when we use the multilevel model to simultaneously estimate the true clutch regression coefficients for all teams, we learn that this observed variability in clutch measures is primarily due to chance. Figure 7 displays the individual and multilevel clutch estimates for all teams. Here the shrinkage is about 93% and we see that the individual clutch estimates are shrunk almost entirely towards the common value in the multilevel model fit. The interpretation is that although teams obviously have different abilities to score runs, teams appear to have similar abilities to advance runners in scoring position.

5.1Introduction

In our exploration of the effects of covariates on run-scoring, we focused on the application of an ordinal regression model with a single covariate. This was done since it is simpler to describe the covariate effects and the shrinkage behavior of the multilevel model. A natural extension of this approach is to model run-scoring of all teams by a large ordinal regression model where one includes all covariates that one believes have an impact on run scoring. Using our general notation, one can write the ordinal logistic model as

where θcj represents the probability of scoring at least c runs for the jth team, x₁, . . , x_k represent k possible covariates (such as league, ballpark, and pitcher quality), and β_j1, . . . , β_jk represent the regression effects for the jth team. In the multilevel modeling framework, one could assign γc1,...,γcN a common multivariate normal distribution, and likewise assume each of the sets of team covariate effects {β_jk, j = 1, . . . , N} come from a common normal distribution. This model is more complicated to fit due to the large number of unknown parameters, but it would accomplish the same smoothing effect as demonstrated in this paper. Team scoring distributions would be shrunk towards an overall scoring distribution – this is accomplished by the common multivariate normal distribution played on the bin cutpoint parameters γcj. In a similar fashion, team effects for a particular covariate such as pitcher quality would be shrunk or moved towards a common covariate effect. This particular approach is promising and can help us better understand relevant covariates that affect run scoring.

5.2Example

To illustrate this multiple regression approach, suppose we are interested in learning about the home and inning effects on team run scoring. For the jth team, we write the ordinal regression model as

This model was fit on the 2013 season data. One can represent an estimate of the pair (β_j1, β_j2) by an ellipse that contains the unknown pair of parameters with 95% probability. Figure 8 displays the estimates of the home and inning effects where one uses individual regression models for the 30 teams, and Fig. 9 displays estimates of the same effects using the multilevel regression model. Note that there are general tendencies for teams to score more at home and in early innings (the ellipses tend to be located in the upper-right section of the graph), but the locations of the ellipses for the multilevel model estimates tend to be shrunk towards a common estimate.