Nuclear penalized multinomial regression with an application to predicting at bat outcomes in baseball

A baseball game comprises a sequence of matchups between one batter and one pitcher. Each matchup, or plate appearance (PA), results in one of several outcomes. Disregarding some obscure possibilities, we consider nine categories for PA outcomes: flyout (F), groundout (G), strikeout (K), base on balls (BB), hit by pitch (HBP), single (1B), double (2B), triple (3B) and home run (HR).

A problem which has received a prodigious amount of attention in sabermetric (the study of baseball statistics) literature is determining the value of each of the above outcomes, as it leads to scoring runs and winning games. But that is only half the battle. Much less work in this field focuses on an equally important problem: optimally estimating the probabilities with which each batter and pitcher will produce each PA outcome. Even for “advanced metrics” this second task is usually done by taking simple empirical proportions, perhaps shrinking them toward a population mean using a Bayesian prior.

In statistics literature, on the other hand, many have developed shrinkage estimators for a set of probabilities with application to batting averages, starting with Stein’s estimator (Efron and Morris,, 1975). Since then, Bayesian approaches to this problem have been popular. Morris, (1983) and Brown, (2008) used empirical Bayes for estimating batting averages, which are binomial probabilities. We are interested in estimating multinomial probabilities, like the nested Dirichlet model of Null, (2009) and the hierarchical Bayesian model of Albert, (2016). What all of the above works have in common is that they do not account for the “strength of schedule” faced by each player: How skilled were his opponents?

The state-of-the-art approach, Deserved Run Average (Judge and BP Stats Team,, 2015, DRA), is similar to the adjusted plus-minus model from basketball and the Rasch model used in psychometrics. The latter models the probability (on the logistic scale) that a student correctly answers an exam question as the difference between the student’s skill and the difficulty of the question. DRA models players’ skills as random effects and also includes fixed effects like the identity of the ballpark where the PA took place. Each category of the response has its own binomial regression model. Taking HR as an example, each batter $B$ has a propensity $\beta_{B}^{\mbox{\scriptsize HR}}$ for hitting home runs, and each pitcher $P$ has a propensity $\gamma_{P}^{\mbox{\scriptsize HR}}$ for allowing home runs. Distilling the model to its elemental form, if $Y$ denotes the outcome of a PA between batter $B$ and pitcher $P$,

(Actually, in detail DRA uses the probit rather than the logit link function.)

One bothersome aspect of DRA is that the probability estimates do not sum to one; a natural solution is to use a single multinomial regression model instead of several independent binomial regression models. Making this adjustment would result in a model very similar to ridge multinomial regression (described in Section 3.3), and we will compare the results of our model with the results of ridge regression as a proxy for DRA. Ridge multinomial regression applied to this problem has about 8,000 parameters to estimate (one for each outcome for each player) on the basis of about 160,000 PAs in a season, bound together only by the restriction that probability estimates sum to one. One may seek to exploit the structure of the problem to obtain better estimates, as in ordinal regression. The categories have an ordering, from least to most valuable to the batting team:

K $<$ G $<$ F $<$ BB $<$ HBP $<$ 1B $<$ 2B $<$ 3B $<$ HR,

with the ordering of the first three categories depending on the game situation. But the proportional odds model used for ordinal regression assumes that when one outcome is more likely to occur, the outcomes close to it in the ordering are also more likely to occur. That assumption is woefully off-base in this setting because as we show below, players who hit a lot of home runs tend to strike out often, and they tend not to hit many triples. The proportional odds model is better suited for response variables on the Likert scale (Likert,, 1932), for example.

The actual relationships among the outcome categories are more similar to the hierarchical structure illustrated by Figure 1. The outcomes fall into two categories: balls in play (BIP) and the “three true outcomes” (TTO). The three true outcomes, as they have become traditionally known in sabermetric literature, include home runs, strikeouts and walks (which itself includes BB and HBP). The distinction between BIP and TTO is important because the former category involves all eight position players in the field on defense whereas the latter category involves only the batter and the pitcher. Figure 1 has been designed (roughly) by baseball experts so that terminal nodes close to each other (by the number of edges separating them) are likely to co-occur. Players who hit a lot of home runs tend to strike out a lot, and the outcomes K and HR are only two edges away from each other. Hence the graph reveals something of the underlying structure in the outcome categories.

This structure is further evidenced by principal component analysis of the player-outcome matrix, illustrated in Figures 2 and 3. For batters, the principal component (PC) which describes most of the variance in observed outcome probabilities has negative loadings on all of the BIP outcomes and positive loadings on all of the TTO outcomes. For both batters and pitchers, the percentage of variance explained after two PCs drops off precipitously.

Principal component analysis is useful for illustrating the relationships between the outcome categories. For example Figure 3(a) suggests that batters who tend to hit singles (1B) more than average also tend to ground out (G) more than average. So an estimator of a batter’s groundout rate could benefit from taking into account the batter’s single rate, and vice versa. This is an example of the type of structure in outcome categories that motivates our proposal, which aims to leverage this structure to obtain better regression coefficient estimates in multinomial regression.

In Section 2 we review reduced-rank multinomial regression, a first attempt at leveraging this structure. We improve on this in Section 3 by proposing nuclear penalized multinomial regression, a convex relaxation of the reduced rank problem. We compare our method with ridge regression in a simulation study in Section 4. In Section 5 we apply our method and intepret the results on the baseball data, as well as another application. The manuscript concludes with a discussion in Section 6.

Suppose that we observe data $\mathbf{x}_{i}\in\mathbb{R}^{p}$ and $Y_{i}\in\{1,...,K\}$ for $i=1,...,n$. We use $\mathbf{X}$ to denote the matrix with rows $\mathbf{x}_{i}$, specifically $\mathbf{X}=(\mathbf{x}_{1},...,\mathbf{x}_{n})^{T}$. The multinomial logistic regression model assumes that the $Y_{i}$ are, conditional on $\mathbf{X}$, independent, and that for $k=1,...,K:$

were $\alpha_{k}\in\mathbb{R}$ and $\beta_{k}\in\mathbb{R}^{p}$ are fixed, unknown parameters. The model (1) is not identifiable because an equal increase in the same element of each of the $\beta_{k}$ (or in each of the $\alpha_{k}$) does not the change the estimated probabilities. That is, for each choice of parameter values there is an infinite set of choices which have the same likelihood as the original choice, for any observed data. This problem may readily be resolved by adopting the restriction for some $k_{0}\in\{1,...,K\}$ that $\alpha_{k_{0}}=0$ and $\beta_{k_{0}}=\vec{0}_{p}$. However, depending on the method used to fit the model, this identifiability issue may not interfere with the existence of a unique solution; in such a case, we do not adopt this restriction. See the appendix for a detailed discussion.

In contrast with logistic regression, multinomial regression involves estimating not a vector but a matrix of regression coefficients: one for each independent variable, for each class. We denote this matrix by $\mathbf{B}=(\beta_{1},...,\beta_{K})$. Motivated by the principal component analysis from Section 1, instead of learning a coefficient vector for each class, we might do better by learning a coefficient vector for each of a smaller number $r$ of latent variables, each having a loading on the classes. For $r=1$, this is the stereotype model originally proposed by Anderson, (1984), who observed its applicability to multinomial regression problems with structure between the response categories, including ordinal structure. Greenland, (1994) argued for the stereotype model as an alternative in medical applications to the standard techniques for ordinal categorical regression: the cumulative-odds and continuation-ratio models.

Yee and Hastie, (2003) generalized the model to reduced-rank vector generalized linear models. In detail, the reduced-rank multinomial logistic model (RR-MLM) fits (1) by solving, for some $r\in\{1,...,K-1\}$, the optimization problem:

If rank$(\mathbf{B})<r$, then there exist $\mathbf{A}\in\mathbb{R}^{p\times r}$, $\mathbf{C}\in\mathbb{R}^{K\times r}$ such that $\mathbf{B}=\mathbf{A}\mathbf{C}^{T}$. Under this factorization, the $r$ columns of $\mathbf{C}$ can be interpreted as defining latent outcome variables, each with a loading on each of the $K$ outcome classes. The $r$ columns of $\mathbf{A}$ give regression coefficient vectors for these latent outcome variables, rather than the outcome classes.

The optimization problem (2) is not convex because rank$(\cdot)$ is not a convex function. Yee, (2010) implemented an alternating algorithm to solve it in the R (R Core Team,, 2016) package VGAM. However this algorithm is too slow for feasible application to datasets as large as the one motivating Section 1 ($n=176559$, $p=796$, $K=9$).

Because of the computational difficulty of solving (2), we propose replacing the rank restriction with a restriction on the nuclear norm $||\cdot||_{*}$ (defined below) of the regression coefficient matrix. For some $\rho>0$, this convex optimization problem is:

We prefer to frame the problem in its equivalent Lagrangian form: For some $\lambda>0$,

This optimization problem (4) is what we call nuclear penalized multinomial regression (NPMR). We use $\ell(\alpha,\mathbf{B};\mathbf{X},Y)$ to denote the log-likelihood of the regression coefficients $\alpha$ and $\mathbf{B}$ given the data $\mathbf{X}$ and $Y$. The nuclear norm of a matrix is defined as the sum of its singular values, that is, the $\ell_{1}$-norm of its vector of singular values. Explicitly, consider the singular value decomposition of $\mathbf{B}$ given by $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$, with $\mathbf{U}\in\mathbb{R}^{p\times p}$ and $\mathbf{V}\in\mathbb{R}^{K\times K}$ orthogonal and $\mathbf{\Sigma}\in\mathbb{R}^{p\times K}$ having values $\sigma_{1},...,\sigma_{\min\{p,K\}}$ along the main diagonal and zeros elsewhere. Then

In the same way that the lasso induces sparsity of the estimated coefficients in a regression, solving (4) drives some of the singular values to exactly zero. Because the number of nonzero singular values is the rank of a matrix, the result is that the estimated coefficient matrix $\mathbf{B}^{*}$ tends to have less than full rank. Thus (4) is a convex relaxation of the reduced-rank multiomial logistic regression problem, in much the same way as the lasso is a convex relaxation of best subset regression (Tibshirani,, 1996). The convexity of (4) makes it easier to solve than (2), and we discuss algorithms for solving it in Sections 3.1 and 3.2. In practice, we recommend using standard cross-validation techniques for selecting the regularization parameter $\lambda$, which controls the rank of the solution.

Consider the singular value decomposition $\mathbf{U}^{*}\mathbf{\Sigma}^{*}\mathbf{V}^{*T}$ of the $p\times K$ estimated coefficient matrix $\mathbf{B}^{*}$. Each column of the $K\times K$ orthogonal matrix $\mathbf{V}^{*}$ represents a latent variable as a linear combination of the $K$ outcome categories. Meanwhile, each row of $\mathbf{U}^{*}\mathbf{\Sigma}^{*}$ specifies for each predictor variable a coefficient for each latent variable, rather than for each outcome category. By estimating some of the singular values of $\mathbf{B}^{*}$ (the entries of the diagonal $p\times K$ matrix $\mathbf{\Sigma}^{*}$) to be zero, we reduce the number of coefficients to be estimated for each predictor variable from (a) one for each of $K$ outcome categories; to (b) one for each of some smaller number of latent variables. These latent variables learned by the model express relationships between the outcomes because two categories for which a latent variable has both large positive coefficients are both likely to occur for large values of that latent variable. Similarly, if a latent variable has a large positive coefficient for one category and a large negative coefficient for another, those two categories oppose each other diametrically with respect to that latent variable.

In this section we present the results of two different simulations, one using a full-rank coefficient matrix and the other using a low-rank coefficient matrix. In both settings we vary the training sample size $n$ from 600 to 2000, and we fix the number of predictor variables to be 12 and the number of levels of the response variable to be 8. Given design matrix $\mathbf{X}\in\mathbb{R}^{n\times 12}$ and coefficient matrix $\mathbf{B}\in\mathbb{R}^{12\times 8}$, we simulate the response according to the multinomial regression model. Explicity, for $i=1,...,n$, and $k=1,...,8$,

For both simulations the entries of $\mathbf{X}$ are i.i.d. standard normal:

for $i=1,...,n$. However the simulations differ in the generation of the coefficient matrix $\mathbf{B}$. In the full rank setting, the entries of $\mathbf{B}$ follow an i.i.d. standard normal distribution: For $k=1,...,8$,

In the low rank setting we first simulate two intermediary matrices $\mathbf{A}\in\mathbb{R}^{12\times 2}$ and $\mathbf{C}\in\mathbb{R}^{8\times 2}$ with i.i.d. standard normal entries, and we then define $\mathbf{B}\equiv\mathbf{A}\mathbf{C}^{T}$ so that the rank of $\mathbf{B}$ is 2. In each simulation we fit ridge regression and NPMR to the training sample of size $n$ and estimate the out-of-sample error by simulating 10,000 test observations, comparing the model’s predictions on those test observations with the simulated response. The results of 3,500 simulations in each setting, for each training sample size $n$, are presented in Figure 4.

In the full rank setting we expect ridge regression to out-perform NPMR because ridge regression shrinks all coefficient estimates toward zero, which is the mean of the generating distribution for the coefficients in the simulation. If this were a Gaussian regression problem instead of a multinomial regression problem, then the ridge regression coefficient estimates would correspond (Hastie et al.,, 2009) to the posterior mean estimate under a Bayesian prior of (10). In fact ridge regression does beat NPMR in this simulation (for all training sample sizes $n$), but NPMR’s performance is surprisingly close to that of ridge regression.

The low rank setting is one in which NPMR should lead to a lower test error than does ridge regression. NPMR bets on sparsity in the singular values of the coefficient matrix, and in this setting the bet pays off. The simulation results verify that this intuition is correct. NPMR beats ridge regression for all training sample sizes $n$ but especially for smaller sample sizes. By betting (correctly in this case) on the coefficient matrix having less than full rank, NPMR learns more accurate estimates of the coefficient matrix. As the training sample size increases, learning the coefficient matrix becomes easier, and the performance gap between the two methods shrinks but remains evident.

In summary, this simulation demonstrates that each of NPMR and ridge regression is superior in a simulation tailored to its strengths, confirming our intuition. Furthermore, in a simulation constructed in favor of ridge regression, NPMR performs nearly as well. Meanwhile NPMR leads to more significant gains over ridge regression in the low rank setting.

The potential for reduced-rank multinomial regression to leverage the underlying structure among response categories has been recognized in the past. But the computational cost for the state-of-the-art algorithm for fitting such a model is so great as to make it infeasible to apply to a dataset as large as the baseball play-by-play data in the present work. Using a convex relaxation of the problem, by penalizing the nuclear norm of the coefficient matrix instead of its rank, leads to better results.

The interpretation of the results on the baseball data is promising in how it coalesces with modern baseball understanding. Specifically, the NPMR model has quantitative implications on leveraging the structure in PA outcomes to better jointly estimate outcome probabilities. Additional application to vowel recognition in speech shows improved out-of-sample predictive performance, relative to ridge regression. This matches the intuition that NPMR is well-suited to multinomial regression in the presence of a generic structure among the response categories. We recommend the use of NPMR for any multinomial regression problem for which there is some nonordinal structure among the outcome categories.

The authors would like to thank Hristo Paskov, Reza Takapoui and Lucas Janson for helpful discussions, as well as Balasubramanian Narasimhan for computational assistance.