3.1.1 Bayesian regression (BR)

Bayesian linear regression is the most common model of choice when assessing the association between a response variable y and p predictors x=(x1,x2,⋯,xp). In the simplest formulation,

where i = 1, 2, ⋯ , n is the observation number and ε is the residual that is assumed to be normally distributed with zero mean and constant variance σ². Priors are then placed on the parameter vector β and σ². See Gelman et al. (2014) for more details. Some examples of the use of this modeling paradigm are given below.

A Bayesian regression approach for small treatment effects, which are commonly encountered in sports performance studies, was proposed by Mengersen et al. (2016) as an alternative to the traditional magnitude-based inference suggested by Batterham and Hopkins (2006). These authors addressed the effect of altitude training regimens on the running performance and blood parameters (hemoglobin mass, maximum blood lactate concentration) in triathlon. They considered G =3 treatments: live high-train low (LHTL), intermittent hypoxic exposure (IHE) and placebo, and 8 participants per group. Another predictor in the model (X) was the change (in %) in training load before and after for each participant.

Letting I₁ and I₂ be indicator values for treatments 1 and 2, respectively, this model can be described in a similar manner to Eq. (2) or in an equivalent hierarchical manner as:

where the i-th observation reside within groups (j), with each group having a potentially different variance.

In another example, Deshpande and Jensen (2016) estimated basketball players’ contributions to their team winning probabilities using a high dimensional regression. Let y_i be the win probability of the home team in the i^th shift, where the shifts are the periods between substitutions. They used the following regression equation:

where μ is the home court advantage. The subscripts h and a represent the home and away teams, respectively. The θ’s are the player’s effect and hence θhi1 and θai1 are the effects from player 1 in the home team and the away team, respectively. For each of the 488 players on the league, they obtained θ estimates. The parameters τHi and τAi are the partial effect associated with the home and away teams. σ denotes a measure of the variability. The marginal posterior densities for θ provide a good picture of the player contribution to winning.

Often, the response y_i is a binary variable following a Bernoulli distribution, for instance whether the player scored or not. This case is frequently approached using a logistic regression model, which can be defined as follows:

Here, ε_i represents extra-Bernoulli variation; see Besag et al. (1995) for details.

Deshpande and Wyner (2017) approached the issue of baseball pitch “framing” from a Bayesian hierarchical perspective using logistic regression. The term “framing” refers to an action often carried out by catchers making a pitch look more like a strike. In order to assess the impact of the catcher in the decision, they estimated the probability that a pitch is called a strike by a given umpire and other covariates such as the pitch location (x, z), the pitcher, etc.

They employed the following logistic model:

where b, ca, co, p and u are the partial effects of the batter, catcher, count, pitcher, and the umpire. The pitch location is given by fu(x,z). Other factors such as the type of pitch (fastball, curveball, etc.) and the speed could also have been included in the model in a straightforward manner. In other examples, Miskin, Fellingham, and Florence (2010) used logistic regression to assess the importance of several skills in volleyball and Cafarelli, Rigdon, and Rigdon (2012) to obtain the probability of converting a third down in National Football League (NFL) based on the number yards to go.

Another useful model for binary response variables in the literature is the probit regression, which is based on the probit link function:

where Φ is the standard normal cumulative distribution function.

Probit regression was used by Jensen et al. (2009b) to construct a baseball defensive model and predict the catching probability given the location of the defender in the field, the velocity of the ball and the direction. They defined their model as:

This model gives the probability of the ball j being caught by the player i. Here, D_ij represents the distance travelled by the player and V_ij is the velocity. The variable F_ij is 1 when moving forward and 0 otherwise. Note that this model considers the interaction between predictors and includes a categorical variable. It also allows for the computation of the player’s contribution to the defence in terms of runs saved, compared to the rest of the players in the same position.

The multinomial logistic regression in McFadden (1973) is a generalized method for cases where a categorical response variable takes more than two values, for example, the possible outcomes of a shot in basketball y={0,1,2,3}. Reich et al. (2006) used this model to assess the relationship between predictors such as the defensive strength of the opposition team and playing home or away, first or second half, and the response variables: (1) location when shooting, (2) the shooting frequency and (3) the efficiency. They let y_i be a region in the court for a given shot i that follows a multinomial distribution with parameter θ(η), and they defined a predictor

where A is a vector of the area of the section j, since some of the sections have different areas. The model is then defined as:

See also Glickman and Hennessy (2015) for another example of the use of a multinomial logit model for competitors’ rank ordering in Alpine skiing competitions.

Bayesian linear mixed models are also on the rise. Revie et al. (2017) considered mixed models to show the viability of players’ perceptions (via surveys) to predict fitness levels for cases when a direct fitness measurement is inconvenient or not possible.

Bayesian non-parametric and semi-parametric regression models have also been employed, albeit less commonly, in sports analysis. For example, a semi-parametric latent variable approach was used by Wimmer et al. (2011) for modeling points performance in decathlon events assuming four latent abilities (sprint, jumping, throwing, and endurance). Another Bayesian nonparametric model which is based on a Dirichlet process mixture model was suggested by Pradier, Ruiz, and Perez-Cruz (2016) for modeling the effect of covariates such as the age, gender and environment in marathon runners’ performance.

Other popular regression techniques are log-linear models. See for instance, Boys and Philipson (2018), who employed an additive log-linear model for ranking cricketers, accounting for factors such as the year and player’s age. Several other log-linear models will be discussed in the football subsection 3.2.2.

3.1.2 Accounting for time

Time is a critical factor in modelling and analysis of many sports (Kovalchik and Albert 2017). The performance of athletes and teams are known to change during the season and even during the course of a game due to factors such as fatigue. Often, it is of interest to analyse factors such as fatigue or momentum that are difficult or impractical to measure. A common approach for such analysis over time is to employ a state space model (SSM) or hidden Markov model (HMM). Both of these models assume that there is an underlying latent variable, z say, that governs the value of the observed variable, y. The form of z determines the type of model employed: if z is categorical or ordinal then a hidden Markov model (HMM) is appropriate, whereas a continuous z leads to a state space model (SSM).

HMMs have been used by several authors (e.g. Albert 1993, Jensen, McShane, and Wyner 2009a, Dadashi et al. 2013, Koulis, Muthukumarana, and Briercliffe 2014). Dadashi et al. (2013) proposed a HMM to estimate timing coordination between hands and feet via estimation of temporal phases of breaststroke swimming. The model uses three axis information from wearable inertial measurement units (IMU) worn on arms and legs, and predicts three hidden states [Q=(q1,q2,q3)] corresponding to glide, propulsion and recovery in leg and arm movements.

Typically, a HMM model is defined as λ = (A, B, π), where A is the state transition probability matrix, B is the emission probability matrix that relates hidden states to observations from the wearable sensor and π is the initial state probability. The HMM model was trained using supervised learning from expert annotated video. The authors reported that the model detected correctly the phases 93.5% of the time in arm strokes and 94.4% in leg strokes. See Dadashi et al. (2013) for a detailed description.

A Poisson HMM model was used by Koulis et al. (2014) for modelling batting performance in cricket. They used a Bayesian approach with multiple states related to the batsman’s performance, where the observed variable is the number of runs per game produced. A Bayesian HMM was also used by Franks et al. (2015) to model basketball defensive placements, where the hidden states are the offensive player being guarded by each defender.

Glickman and Stern (1998) suggested a Bayesian state-space approach to predict American football teams strengths using a first-order auto-regressive process. They found that this model was able to predict the outcomes of games outcomes slightly better than the Las Vegas Betting Line oddsmaker. A nonlinear version of this model was suggested by Glickman (2001) to evaluate paired comparisons in NFL football and chess.

Some other approaches accounting for time were suggested by Stephenson and Tawn (2013) and Kovalchik and Albert (2017). Stephenson and Tawn (2013) applied concepts of extreme value theory to model annual best racing times in athletics considering an exponentially decreasing trend. This facilitates the comparison of athletes who performed in different decades. In tennis, Kovalchik and Albert (2017) fitted temporal data (time-to-serve) using a Bayesian hierarchical model and the covariates point importance and the length of the previous rally.

3.1.3 Accounting for space and time

As discussed above, modern tracking technology is providing the location of players and the ball at regular small intervals of time. This is opening the door to spatio-temporal analysis, in which the court (basketball), the course (golf) or the field (baseball) is often discretized using a grid. This grid is often a square in the Cartesian systems, e.g. one-square-foot quadrats (Miller et al. 2014) or a slice in the polar coordinates system, e.g. Reich et al. (2006); Yousefi and Swartz (2013).

In the case of basketball, a common task is to compute the probability of scoring as a function of the location. This generally yields a heat map of scoring probabilities. For instance, Reich et al. (2006) used logit multinomial Bayesian regression to assess the relationship between the shot location in the court and some covariates such as the presence of key players from the same team in the court, defensive strength, playing home or away, etc. They also assessed the significance of these predictors on shooting frequency and efficiency in different regions measured using polar coordinates (the distance to the basket and the angle). They used conditionally autoregressive (CAR) and two neighbor relation CAR priors to achieve a smoother surface borrowing information from neighbors.

Shortridge, Goldsberry, and Adams (2014), for example, extended this idea by computing the spatial variability in scoring within an empirical Bayesian framework. They used a shrinkage approach to obtain a smoother scoring probability surface. Spatial shooting patterns have been also modelled using a log-Gaussian Cox process (LGCP) (Miller et al. 2014; Franks et al. 2015). Also in basketball, Cervone et al. (2016) suggested the use of a conditional autoregressive model to compute the expected score as a function of factors such as the player in possession of the ball, defensive stance, etc. Here the CAR prior accounts for spatial autocorrelation by adding a random effect for the player. See Cervone et al. (2016) for more details.

In baseball, Jensen et al. (2009b) used a hierarchical model to estimate the probability of a defensive player catching a ball considering among other parameters the location of the player. Pitch horizontal and vertical coordinates around the strike zone were considered for estimating the probability of a strike Deshpande and Wyner (2017). Yousefi and Swartz (2013) introduced a metric for golf putting performance considering the distance to the pin and the angle. This approach splits the “green” area into eight slices centred at the pin, where the within slice probability of scoring is dependent on the distance.

3.1.4 Other methods

Variations and extensions of the above general classes of models have also been used in sports analytics. For instance, Swartz, Gill, and Muthukumarana (2009) developed a simulator for predicting one-day cricket games outcomes using a latent variable approach. Ofoghi et al. (2013) considered the selection of racing athletes in the multi-event cycling omnium using Bayesian Networks.

See also for example Stenling et al. (2015) for Bayesian structural equation model (SEM) with application to sport psychology settings. The authors estimated several latent factors associated with the athlete’s behavioral regulation measured with the Sport Motivation Scale. They found a better fit to the data using a Bayesian approach compared to the traditional method based on the maximum likelihood.

Empirical Bayes is a very popular statistical technique in which the prior distribution in Eq. (1) is obtained from observed data e.g. obtained from previous games or from players with similar characteristics or the same position. This makes EBA to be considered as pseudo or not fully Bayesian.

EBA models are generally hierarchical where the parameter of interest is assumed to come from a common pooled distribution. EBA is particularly useful in part because it leads to fast computation. For decades EBA has enjoyed consistent use in baseball for modeling batting averages Efron and Morris (1973); Brown (2008); Neal et al. (2010); Jiang et al. (2010). For instance, estimates of baseball averages of players with a few at-bats can be obtained using the league average as a prior distribution. See an extended discussion in Robinson (2017). In basketball, spatial models of shooting effectiveness are commonly built using EBA because a smooth scoring intensity surface can be obtained (e.g. Shortridge et al., 2014). Another example can be found in Baker and McHale (2017) who recently suggested an approach for estimating player strengths based on empirical Bayes.

Finally, another area that is of interest in sports analytics is experimental design. See, for example, Glickman (2008) who developed a Bayesian locally optimal design approach for knockout-based competitions.

3.2.1 Basketball

Basketball is one of the most popular, dynamic and competitive sports worldwide. In this game, two teams of players interact in different locations of the court according to a given set of rules with the aim of scoring in the opponent team’s basket.

An early paper by Reich et al. (2006) used logit multinomial Bayesian regression to assess the relationship between the shot location in the court and some covariates such as the presence of key players from the same team in the court, defensive strength, playing home or away, etc. The inference in this work was limited to a one NBA player (Sam Cassell) during the season 2003–2004.

The adoption of the SportVU player tracking technology after 2010 in the NBA marks a milestone in basketball analytics. It enhanced the individual statistics levels by capturing (at 25 frames per second) the coordinates of each player (x, y) and the ball (x, y, z). Detailed statistics such as the players’ distance covered in a game and the speed developed when approaching the basket, became available as result.

The success of the paper by Goldsberry (2012) on spatial modeling of shooting effectiveness motivated a large number of publications in this area. See e.g. Shortridge et al. (2014) who suggested metrics like the expected number of points per shot for a given location within the offensive court and points above league average for each player.

Other spatial models that explicitly incorporate spatial information have also been proposed. For example, Miller et al. (2014) employed a log-Gaussian Cox process as a spatial prior and combined this with dimension reduction to obtain the players’ shooting intensities and the identification of shooting habits.

Cervone et al. (2016) introduced a statistic called expected possession value (EPV) plotted as the expected number of points in a given offensive play that a team might score versus time (from 0 to 25 seconds). This metric depends on the player in possession of the ball, its location, the defense placement, etc. The contribution of players to the team’s win probability was assessed by Deshpande and Jensen (2016) using Bayesian linear regression model. See also Lam (2018) who used a Bayesian regression to predict the outcome of NBA basketball games.

3.2.2 Football

Numerous studies have attempted to model the outcome of football matches (Rue and Salvesen 2000; Karlis and Ntzoufras 2008; Baio and Blangiardo 2010). For instance, Karlis and Ntzoufras (2008) used a Poisson difference distribution to model the difference of goals in football games using data from the English Premier League. Let X_i and Y_i be the score of the home and away teams in the i^th game. They defined a statistic Z_i as follows

where PD is the Poisson difference distribution with rates λ_1i and λ_2i, that are obtained using the following log-linear link functions:

where μ is a constant parameter. H is the home team coefficient. A and D are the parameters for the team attack and defense.

Baio and Blangiardo (2010) suggested some improvements on Karlis and Ntzoufras (2008) model to predict football results in the Italian Serie A championship. They obtained the number of goals from each team using a Poisson distribution rather than modeling the difference as suggested by Karlis and Ntzoufras (2008). The model produces estimates of the posterior distributions of attack and defense.

Suzuki et al. (2010) suggested a Bayesian model for forecasting the results of the 2006 World Cup taking into consideration the FIFA World Ranking. In this approach, the number of goals on each team is fit using Poisson distributions.

The values R_A and R_B are the ratings of the team A and B, respectively. The prior distributions for λ_A and λ_B were set as Gamma distributions, and expert knowledge was incorporated via elicitation. This approach does not consider other relevant variables such as the offensive/defensive skills.

Other papers used Bayesian methods for assessing the importance of player skills (Thomas, Fellingham, and Vehrs 2009), functional performance (Carvalho et al. 2017), optimum substitution times (Silva and Swartz 2016).

3.2.3 Baseball

A wide range of Bayesian techniques has been used in baseball. Predicting batsmen averages has fascinated researchers and statisticians for a long time and “several authors recently have taken a swing at the subject” (Neal et al. 2010). For instance, Efron and Morris (1973); Brown (2008); Neal et al. (2010); Jiang et al. (2010) predicted baseball averages within an empirical Bayes approach. In another examples (Albert 1993, 2008) used hidden Markov models for assessing streakiness among batsmen.

Jensen et al. (2009a) used a log-linear hierarchical model to predict player’s number of home runs per season. Age, player position and home ballpark are predictors included in the model along with previous seasons performance. A mixture model on the intercept term is used to create groups of home run hitters (elite and non-elite). The probability that a hitter is a member of the elite group is determined using a hidden Markov model. This model was shown to have better predictive accuracy than other competing methods. McShane et al. (2011) used also a hierarchical Bayesian model for the selection of performance variables that better describes offensive abilities.

A good defense is critical for winning games. However, this aspect is difficult to quantify since the traditional assessment is rather subjective, making it hard to compare the contribution of the players. Jensen et al. (2009b) addressed this issue using a Bayesian probit regression model for assessing the fielder’s effectiveness. They computed the player’s contribution to the defense in terms of runs saved, compared to the rest of the players playing the same position. The model predicts the catching probability given the location of the defender in the field, the velocity of the ball and direction that the fielder has to move to (forward or backward). It would be interesting to see an extension of this analysis considering baseball park’s constraints.

Healey (2017) suggested new statistics for players performance based on batted-ball parameters (speeds, vertical and horizontal angles) within the Bayesian philosophy. Probability density estimates are obtained using a non-parametric approach. In another example, Bendtsen (2017) suggested Bayesian networks for modeling career regimes.

3.2.4 Other team sports

In ice hockey, Thomas (2006) modeled the scoring probability as a continuous-time Markov process, and Gramacy, Jensen, and Taddy (2013) introduced a Bayesian logistic regression model for evaluating the impact of hockey players on team’s scoring. This last model is an alternative to the plus-minus approach and is based on a Laplace prior for the regression coefficients to facilitate variable selection, which is required in these high dimensional regression problems characterized by a large number of players. See also Thomas et al. (2013) who introduced a method for determining players abilities by modelling the team’s scoring rate as a semi-Markov process using hazard functions.

Giles et al. (2017) used a Bayesian regression model to assess the association between mental toughness and behavioral perseverance accounting for physical fitness in Australian rules footballers. They found association between these variables except in presence of fatigue. Multiple examples can also be found in cricket. See e.g. Damodaran (2006); Brewer (2008); Swartz et al. (2009); Boys and Philipson (2018).

3.2.5 Other related issues

Streakiness: Another cluster of publications addressed the issue of streakiness (also known as the hot hand phenomenon). Say, for example, when a baseball player shows a pattern indicating a substantially larger (than average) proportion of hits (successes) in a period of time. This apocryphal phenomenon is supposed to be experienced by athletes during the season and it has been largely studied among others by Gilovich, Vallone, and Tversky (1985); Albright (1993); Bar-Eli, Avugos, and Raab (2006).

Albert (1993), for example, inspired by the thesis of Albright (1993), used two-state hidden Markov chains while Albert (2008) employed the Bayes factor for detecting non-random changes in the batting performance. A similar approach was followed by Wetzels et al. (2016) who analyzed the streakiness rates in basketball. Yang (2004) suggested a Bayesian binary segmentation method for analyzing consecutive successes or failures. This method relies on the Bayes factor to assess the change in the success rates. The author analyzed several popular events considered to be the result of a streaky performance in basketball, baseball and golf. Whether the player missed the previous shot has been considered by Reich et al. (2006) as a predictor of the shooting frequency and location in the court in basketball. However, they found no relationship between them.

Doping: Antidoping studies on biological markers are generally addressed using a longitudinal approach within a Bayesian framework to account for within and between athlete variation. Elements of Bayesian inference are particularly useful, ranging from the established population-based reference antidoping approach to an individual passport system.

Sottas et al. (2006), for example, suggested a method for the detection of abnormal T/E (testosterone glucuronide/epitestosterone glucuronide) ratio values. This approach compares the test results against a cutoff threshold obtained using Bayesian inference and the estimated population and intraindividual mean and coefficient of variation. Robinson et al. (2007) extended this approach for the detection of another illegal drug (recombinant human erythropoietin) and Schulze et al. (2009) added genotype (UGT2B17) as a predictor into a Bayesian framework suggested by Sottas et al. (2006) to achieve an increased test sensitivity. Bayesian inference is also used by Van Renterghem et al. (2011) for the detection of testosterone based on new biomarkers.

Relative age effect: The Relative Age Effect (RAE) establishes that children/athletes who were born in the first months after the school year cutoff have more chances of success. Ishigami (2016) used a Poisson Bayesian regression model to investigate the impact of the RAE and birthplace on the chances of becoming a professional athlete in Japan. The author reported that those who were born in the first month after the cutoff were three times more likely to become a professional athlete.