Bayesian isotonic logistic regression via constrained splines: an application to estimating the serve advantage in professional tennis

08/29/2019 ∙ by Silvia Montagna, et al. ∙ Università di Torino 0

In professional tennis, it is often acknowledged that the server has an initial advantage. Indeed, the majority of points are won by the server, making the serve one of the most important elements in this sport. In this paper, we focus on the role of the serve advantage in winning a point as a function of the rally length. We propose a Bayesian isotonic logistic regression model for the probability of winning a point on serve. In particular, we decompose the logit of the probability of winning via a linear combination of B-splines basis functions, with athlete-specific basis function coefficients. Further, we ensure the serve advantage decreases with rally length by imposing constraints on the spline coefficients. We also consider the rally ability of each player, and study how the different types of court may impact on the player's rally ability. We apply our methodology to a Grand Slam singles matches dataset.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Predicting the outcome of tennis matches has attracted much attention within sport analytics over the years for a number of applications. For example, prediction models can provide coaches useful feedback about how players are improving over time and who they should be able to beat. Further, prediction models could help assess fan engagement and determine who is the favourite player, by how much, and who is currently the best player. See, for example, Glickman (1999); Klaassen and Magnus (2003); Barnett and Clarke (2005); Newton and Keller (2005); Gilsdorf and Sukhatme (2008); Gomes et al. (2011); Smith (2013); Irons et al. (2014); Kovalchik (2016) and references therein.

It is nowadays generally acknowledged that the service is one of the most important elements in tennis. Indeed, it has been observed that the serving player wins more points than the receiving player in elite tennis (Lees, 2003). With the advances in racquet technologies, most top male players can hit service speeds of over 200 Kph. Kotze J. and Rothberg (2000) point out that if the serving speed reaches the receiver’s reacting threshold, it becomes virtually impossible for the receiving player to return the ball. In the extreme, a strong serve strategy that gets rarely broken reduces the competitiveness of the game, and this may result in a loss of spectator interest. For this reason, the International Tennis Federation (ITF) monitors the importance of the serve and can undertake measures, such as slowing surface speeds, to ensure the game’s combativeness is not endangered.

While it is reasonable to assume that the serve advantage gets lost as the rally length increases, there are only a few contributions in the literature attempting to quantify the serve advantage and relate it to rally length via a statistical model. An early contribution is given by O’Donoghue and Brown (2008), where the authors describe the advantage of serving in elite tennis by comparing points won by both the server and the receiver for a given rally length. They conclude that the serve advantage is lost after the 4th rally shot on men’s first serve. Subsequently, Kovalchik (2018b) proposes a Bayesian hierarchical model to estimate player-specific serve curves that also adjust for the opponent rally abilities. In particular, the author uses an exponential decay function to model the decline in serve advantage plus a random effect, representing the difference between the rally ability of the opponents.

In this paper, we focus on the role of the serve advantage in winning a point as a function of the rally length. Our approach falls into the Bradley-Terry class of models (Bradley and Terry, 1952) and is built upon of Kovalchik (2018b). We propose a Bayesian isotonic logistic regression model by representing the logit of the probability of winning a point on serve, , as a linear combination of B-splines basis functions, with athlete-specific basis function coefficients. We point out that while the term isotonic is used to denote regression models where monotonicity is imposed everywhere, in our application we may also want to accommodate for monotonicity only in a subinterval of the function domain. The smoothness of is controlled by the order of the B-splines, while their shape is controlled by the associated control polygon (de Boor, 2001). In particular, to ensure the serve advantage is non-increasing with rally length we constrain the spline function to be non-increasing by controlling its control polygon. This essentially results in imposing a constraint on the coefficients of the spline function. Further, we allow for the probability to win on serve to also depend on the rally abilities of the opponents. We note that the rally advantage component of the model draws on Kovalchik (2018b), but we extend it further to study how the different types of court (e.g., clay, hard) may impact on the player’s rally ability. It is indeed well known that some players favour and perform better on particular surfaces (e.g., Nadal holds 11 French Open (clay) titles and 2 Wimbledon (grass) titles). Each surface material presents its own unique characteristics and provides different challenges to the players, with certain playing styles working better on some types of court and less effectively on others. For example, a grass court is the fastest type of court because of its low bounce capacity. Players must get to the ball more quickly than with clay or hard courts, thus players with stronger serve will generally perform better on grass. The rally advantage component of our model reflects how a player is likely to perform on a particular surface, and this in turns affects the win probability.

Our contribution is twofold: first, the basis function decomposition allows for a more flexible modelling of the longitudinal curve for serve advantage than that attainable via an exponential decay function. Our hierarchical Bayesian framework further accommodates for the borrowing of information across the trajectories of the different athletes, and allows for out-of-sample prediction; second, our construction allows for the inclusion of covariates (e.g., the type of terrain) in the modelling of the rally abilities of the opponents. Therefore, it becomes possible to examine how covariates impact on the rally abilities, and ultimately on the distribution of the serve advantage curves.

The remainder of the paper is organised as follows. In Section 2 we provide an overview of the data used for our analysis. In Section 3 we present our hierarchical Bayesian isotonic logistic model. In Section 4 we compare our model with Kovalchik (2018b). Section 5 presents the results of our real data analysis, and conclusions are outlined in Section 6.

2 Grand slam data

We consider point-by-point data for main-draw singles Grand Slam matches from 2012 forward. Organisations such as the ITF and Grand Slam tournaments record some data on professional tennis matches, but rarely make it available to the public. In this paper, we use data scraped from the four Grand Slam websites shortly after each event by Jeff Sackmann111https://github.com/JeffSackmann. The data is also available with the R package deuce (Kovalchik, 2018a).

There are four Grand Slam tournaments, namely, the Australian Open, French Open, Wimbledon, and US Open. These tournaments are subdivided in two types of associations: the Association of Tennis Professionals (ATP), containing all the matches played by male athletes, and the Women’s Tennis Association (WTA), containing all matches played by female athletes. We consider the male and female tournaments separately, thus obtaining two datasets. For both datasets, we include players with three matches or more in the training data. For a more robust inference, we consider only rally lengths between zero and thirty, counting as zero the first shot played by the server. For both datasets, we extract the following variables: rally length, the series of return hits of the ball from a player to the opponent (an integer between zero and thirty); the ID (name and surname) of the players serving and receiving, respectively; an indicator variable denoting if the server wins the point; the tournament name, used to derive the type of court in the different tournaments. Indeed, the Australian Open and US Open tournaments are played over a hard court, the French Open is played on clay while Wimbledon is played on grass. Unfortunately, other information of potential interest, e.g. the serve’s speed and direction, is not available for every rally. Table 1 reports some summary statistics about the dataset. In Table 2 we report the total number of rallies in the ATP and WTA tournaments, respectively. Short rallies, i.e. rally lengths smaller than or equal to , constitute of the rallies played during the Grand Slam tournaments.

Tournament
US Open Aus Open French Open Wimbledon
Matches
ATP 147 106 192 214
WTA 158 106 186 118
Players
ATP 108 107 132 107
WTA 108 127 119 104
Table 1: Number of matches and players in the four Grand Slam tournaments by association from 2012 forward. The data was collected by Jeff Sackmann, and is also available with the R package deuce (Kovalchik, 2018a).
Short rallies Long rallies Total
ATP 130577 14933 145510
WTA 71592 10288 81880
Table 2: Number of rallies by tournament. We define as short a rally whose length is less than or equal to 4.

In Figure 1

we report the observed relative frequency of rallies won by the server given the number of shots. It appears that the server has a higher chance of winning the point on odd-rally lengths compared to the even-rally lengths. This pattern can be explained by the fact that even-numbered rallies end on the server’s racquet, so he/she can win or make a mistake. Since the

-axis report the observed relative frequency of winning for the server, all the even-shots in Figure 1 represent the case in which the server wins a point with a winner. A winner is a shot that is not reached by the opponent and wins the point. Occasionally, the term is also used to denote a serve that is reached but not returned into the court. On the other hand, the odd-numbered rallies are the winners or errors made by the receiver. In particular, the odd-shots in Figure 1 represent the errors done by the receivers.



Figure 1: Observed frequency of rallies won by the server given the number of shots for the ATP tournament (left) and the for the WTA (right). The blue points represent the odd-shots, while pink points are the even-numbered rallies. The size of a point is proportional to the number of the server’s victories with rally length equal to divided by the total number of points won by the server.

Because errors are more common than winners, we aggregate odd and even rally lengths (see also Kovalchik (2018b)

). We obtain a vector of integers, where 1 corresponds to rally lengths equal to zero or one, 2 corresponds to values 2 or 3 of rally length and so on. This ensures that the same set of outcome types for the server and receiver are represented within each group. The resulting frequencies are showed in Figure

2. We observe that after the first shot the server’s chance of winning the point drastically decreases. This is clear for both men and women. As conventional wisdom suggests, the server has the highest chance of winning a point at the beginning of the rally, owing to the strength of the serve. As the rally progresses, the serve advantage is expected to get increasingly small and have increasingly less influence on the outcome of the rally with each additional shot taken.



Figure 2: Conditional percentage of winning a point given the number of shots for the ATP tournament (left) and the for the WTA (right). Since we aggregated the odd and even results, rally length is between 1 and 15. The size of a point is proportional to the number of the server’s victories with rally length equal to divided by the total number of points won by the server.

3 Hierarchical Bayesian isotonic logistic regression model

In this section, we present our Bayesian hierarchical isotonic regression approach to model the serve advantage. Let be the discrete variable representing rally length, where is the set of all integers between and . Let

be a binary random variable which is equal to one if server

wins the point against receiver , and zero otherwise. We assume

(1)

where is the probability that server wins a point against receiver at rally length , e.g . We consider two components to model , the first describing the serve advantage and the second representing the rally ability of the players. Specifically, we model the logit of as follows:

(2)

with

(3)

for and , where is the total number of servers and the total number of receivers. The ’s are athlete-specific parameters representing the rally ability of the player. We observe that the function ( index omitted for simplicity) is defined for each in the continuous interval , however only its values , for in the discrete set , enter the sampling model in Equation (2). The continuous structure of takes into account the overall trend of the serving advantage (i.e., it estimates the drop of the serving advantage as the rally length increases), accommodating for the longitudinal structure of the data.

The conditional log odds for the probability of the -th server winning a point against the -th receiver is a non-linear function of the rally length, . Function in Eq. (3) represents the decay part of the model via a linear combination of basis functions , where is the dimension of the spline basis. In particular, we opted for B-splines basis functions of order on , and the ’s are athlete-specific basis function coefficients. Specifically, let and be a non-decreasing sequence of knots such that for all , and and . Function is a linear combination of the B-splines , and is called a spline function of order and knot sequence (de Boor, 2001). In other words, each is a piecewise polynomial of degree with breakpoints , and the polynomials are times continuously differentiable at . Here denotes the cardinality of . Moreover, we recall that spline has support on the interval and here we are going to assume and . A more extensive presentation of spline functions is given in de Boor (2001).

In our model, the spline function is defined on the whole interval and does not go to zero for high values of rally length. Indeed, by looking at the last value of rally length in our application, e.g. , it is clear that the logit of the conditional probability of winning the point against reduces to

We consider this as an asymptote, describing the server’s (-th player) log-odds of winning a point against opponent when the serve advantage has vanished. We can interpret parameter as the rally ability of the - player. When the serve advantage vanishes, the probability of winning the point depends on the discrepancy between the rally abilities of the two players plus a constant, obtained from the B-splines basis.

In the next sections, we provide more details regarding the modelling of the serve advantage and the rally ability, respectively.

3.1 Modelling the serve advantage

Hereafter we will denote a spline function as partially monotone if is monotone only in a sub-interval of its domain , for example if is monotone decreasing in . To simplify the notation, we will omit index that denotes individual-specific objects, so we will let and . Further, the spline coefficients will be also referred to as control points in the following.

Given that the serve advantage is expected to decrease as rally length increases, the spline function should be non-increasing in . While the non-increasing behaviour can be directly learnt from the data for small values of rally length, this could be harder to achieve for large values of rally length due to data sparsity in this part of the function domain. In other words, we may have to impose that the spline function is non-increasing for large values of rally length. Given a threshold , we may allow to be free to vary for small values of rally length (i.e., for ), while it is crucial to ensure that is non-increasing as goes above . Thus, we would like to be partially monotone, according to our definition. Then, we need to investigate which condition the spline function must verify to guarantee the partial monotonicity constraint. To this end we will first provide the following definition.

Definition 1 (Control Polygon).

Let be a non-decreasing sequence of knots and let be a spline function of order and knot sequence . The control polygon of is defined as the piecewise linear function with vertices at , where

is called the th knot average.

We note that because it is assumed that for all . The left panel of Figure 3 provides an illustrative example of a spline function with its associated control polygon. The spline function has order and knot vector

The control points

are randomly drawn from standard Normal distribution. The control polygon approximates the spline function

, and the approximation becomes more accurate as the number of control points increases.

Figure 3: Examples of spline functions (black solid lines) and associated control polygons (piecewise linear lines). In both panels, the spline functions were generated assuming B-splines functions of order on and with knot sequence , and is the dimension of the spline basis. Left panel: control points are generated as , for . Right panel: the control polygon is restricted to be non-increasing in , and the resulting spline function is such from the smallest knot greater than . Black rug bars indicate the knot averages, while orange rug bars denote the interior knots.

To ensure the serve advantage is non-increasing with rally length, we need to control the shape of the spline function on an interval . To do so, we will follow the notation and the construction of Abraham and Khadraoui (2015) hereafter. In particular, for all we denote by the smallest knot greater than , with the proviso that if belongs to . Further, we assume that for all so that for all . We identify splines as those whose support intersects with . In other words, for and for , and the spline function restricted to the interval reduces to

(4)

We define the restricted control polygon on as the control polygon associated to the spline (4

), that is, the piecewise linear function that interpolates the vertexes

for . We denote by this restricted control polygon, which is defined for . We also observe that is defined on a interval that contains . Indeed and .

The spline function can be restricted to be non-increasing on by imposing that the associated restricted control polygon is non-increasing as the following Proposition states.

Proposition 1.

Let be a spline function with , and let be the associated control polygon. Consider a real such that . If the restricted is non-increasing on its support , then the restricted spline function is also non-increasing on .

The proof of Proposition 1 is given in Appendix A. Abraham and Khadraoui (2015) remark that if the control polygon is unimodal in , then is unimodal or monotone on , but it is possible to force to be unimodal by increasing the number of knots.

Controlling the shape of the control polygon reduces to controlling the magnitude of the sequence of control points . For example, for a non-increasing constraint on , the broken line with vertexes is non-increasing if . In particular, we impose that the spline coefficients of the restricted spline satisfy

(5)

Thus we have that the spline coefficients are free, while the spline coefficients must be chosen such that condition in Equation (5) is satisfied. The right panel of Figure 3 shows an example of spline function whose shape is constrained to be non-increasing on . Given knots , it is simple to realise that . Thus, , for , while the remaining coefficients are such that . The resulting control polygon is decreasing from , and the spline function decreases from .

3.2 A new prior for the isotonic model

Now, let us return to the athlete-specific notation, that is, let’s denote with the spline function and with the spline coefficients for athlete . Hereafter, we will construct the spline function as (an intercept plus) a combination of B-spline functions defined on the closed set , with order and knot vector . For rally lengths , the decreasing behaviour of the spline function is learnt from the data (Figure 2), whereas the non-increasing trend for rally length larger is assured by controlling the trend of to the restricted spline .

In order to specify a Bayesian model which takes into account the constraints on , we have to specify a prior distribution on the spline coefficients such that the conditions described in section 3.2 are satisfied. With this goal in mind, the free spline coefficients are given a Normal prior distribution with mean

and variance

:

(6)

The prior mean and precision

are given conditionally conjugate prior distributions:

(7)

where is the mean and is the variance of . Further:

With regard to the constrained coefficients, we need to ensure the condition in Equation (5) is verified. Therefore, we define these parameters recursively by letting

(8)

where are random decrements with

(9)

The last equation shows how a spline basis can be easily constrained to be non-increasing, while still retaining its essential flexibility.

For our application, we choose the same setup discussed in the example of Section 3.1. The constrain satisfies the empirical conclusion of O’Donoghue and Brown (2008), namely, that the serve advantage is lost after the 4th rally shot. We assume independent Normal priors as in Equation 6 for the first three spline () coefficients and we adopt the recursive construction (8), with the prior in (9), for the remaining () coefficients.

3.3 Modelling the rally ability

As outlined above, parameter in Equation (2) can be interpreted as the rally ability of server . It is clear that parameters and in Equation (2) are not identifiable, that is, adding and subtracting a constant to these parameters leaves , thus inference, unchanged. Non-identifiability of the ’s is not a concern if one is solely interested in learning the logit of . However, we are also interested in direct inference of the rally ability parameters, thus we need to include an identifiability constraint (Gelfand and Sahu, 1999; Baio and Blangiardo, 2010). In particular, we adopt a sum-to-zero constraint by setting

where is the total number of players in the dataset. We specify a Gaussian prior distribution for the rally ability parameter as in Kovalchik (2018b):

(10)

Finally, we specify the following conditionally conjugate non-informative priors on the hyperparameters:

3.4 Estimating a court effect

In the Grand Slam tournaments, players play over three different types of court: clay, grass, and hard, respectively. The tennis season begins with hard courts, then moves to clay, grass, and back again to hard courts. While very popular in the past, nowadays only Wimbledon is played on grass. Each surface elicits different ball speed, bounce height, and sliding characteristics. For example, grass courts produce little friction with the ball, which will typically bounce low and at high speed on this court. Conversely, clay slows down the ball a little and allows players more time to return it, resulting in longer rallies. Players have to adapt their technique effectively to the surface. However, adapting training and playing schedules is extremely physically demanding on the player. As a result, it is very difficult for one player to dominate across all the courts, and thus all the slams (Starbuck et al., 2016).

It is therefore reasonable to state that the surface type can impact on a player’s performance. Gorgi et al. (2018) study the effect of the different courts for ATP players using a Bradly-Terry model and conclude that taking this information into account leads to improved rankings of the players. In our model, it is straightforward to include court as a covariate within a regression model for the rally ability of each player, and observe the best player for each court. In particular, we define the probability of winning a point on serve given both the rally length and the type of court:

(11)

for and , where is the total number of servers, the total number of receivers, and is the dimension of the splines basis. Here index denotes the type of court, with , where 1 means clay, 2 stands for grass and 3 means hard.

Adding this covariate to our model does not affect the serve advantage, which is modeled as in Section 3.1. Conversely, we now have a subject-specific vector of rally abilities , where refers to court type . To ensure the ’s are identifiable for all players, we impose

where is the total number of players in the dataset. We specify a Gaussian prior distribution on the rally ability parameters:

(12)

for , , and for and . Finally, we specify the following conditionally conjugate non-informative priors on the hyper-parameters:

4 Model comparison

In this section, we compare different models in order to identify the best on the tennis data. In particular, we consider the pair-comparison exponential decay model, proposed by Kovalchik (2018b), and three versions of our Bayesian isotonic logistic regression (BILR) model: 1) no constraints on the spline coefficients, thus the spline function is free of monotonicity constraints; 2) set () and impose an order constraint on the coefficients of the B-splines with support in , thus the resulting spline function is non-increasing in (partially monotone); and 3) spline function constrained to be non-increasing in , with an order constraint on all basis function coefficients, . Setting 2) draws on O’Donoghue and Brown (2008), who observe the serve advantage is lost after the 4th rally shot on men’s first serve.

To compare the performance of the different methods, we compute four goodness of fit indices broadly used in the Bayesian framework. In particular, we consider the Log Pseudo Marginal Likelihood (LPML) (Geyser and Eddy, 1979)

, which derives from predictive considerations and leads to pseudo Bayes factors for choosing among models. Further, we compute the Deviance Information Criterion (DIC)

(Spiegelhalter and der Linde A., 2002), which penalizes a model for its number of parameters, and the Watanabe Akaike information criterion (WAIC) (Watanabe, 2010). The latter can be interpreted as a computationally convenient approximation to cross-validation and it is not effected by the dimension of the parameter vector. Finally, we also compute the root mean squared error (RMSE).

Prior to implementing the BILR model, one has to choose the order of the B-spline bases , the number of knots and their location, which together determine the dimension of the spline basis. We recall that in Eq. (3) is determined as number of interior knots. Further, one has to choose the sub-interval of the spline function domain where monotonicity is to be imposed. For setting 2) above, this sub-interval is chosen to be . We performed some preliminary sensitivity analysis to investigate changes in performance of the BILR model due to different choices for the number of bases, their degree, the knots location, and range . Results of our sensitivity analysis are reported in Orani (2019). For all three versions of our BILR model 1)-3) above, the spline functions , with , are constructed from a B-spline basis of dimension as in Eq. (3), defined on the closed interval . The spline functions have order and knot vector . This model maximises the LPLM criterion and minimises the DIC and the WAIC, as reported in Orani (2019).

In Figure 4 we compare the fit obtained for Andy Murray with the exponential model and our three versions of our BILR model. The exponential decay model (top left) has a decreasing behaviour until , and after that the chance of winning a point is a constant given by the difference between Murray’s rally ability and the average rally ability of his opponents, i.e. to plot the figure we substituted in Eq. (2) by obtained averaging all the for . Our BILR model with no constraints (top right) allows for an increasing behaviour in the probability of winning the point for some intermediate values of rally length, and this behaviour is unlikely to be justified in practice. While for small values of rally length the decreasing trend in serve advantage can be learned from the data, for large values of rally length this behaviour must be imposed through the model given data sparsity. We recall that short rallies, i.e. , constitute

of the rallies in the dataset. In this data-rich part of the domain, no constraint is needed to adequately describe the data. Conversely, for long rallies, the decreasing behaviour imposed via the prior on the coefficients leads to a model which is not influenced by outliers. Both the model with monotonicity constraint in

(bottom left) and the model with all spline coefficients constrained to be non-increasing (bottom right) display a non-increasing behaviour for large values of rally length.





Figure 4: Probability of winning a point as a function of rally length for Andy Murray estimated with the exponential decay model (top left), BILR with no constrained splines (top right), BILR with order constraint on the coefficients of the splines with support in (bottom left), and BILR with order constraint on all spline coefficients (bottom right). The points represents the real data, while the black lines are the posterior mean estimate of the probability of winning a point as a function of rally length obtained with these models. The blue dashed lines are the credible intervals.

For a quantitative evaluation of the performance of the four approaches, we compute the goodness-of-fit measures for these models, reported in Table 3.

Goodness-of-fit-Measures
LPML WAIC DIC RMSE
Exponential model -52813.7 105821.3 105876 22.52
Spline models
No constraints -52760.4 105853.9 105818 20.71
Non-increasing in -52739.1 105747.6 105828 19.32
Non-increasing in -52744.9 105790.4 105875 20.37
Table 3: Predictive goodness-of-fit measures for different model settings.

Although no dramatic difference in performance emerge, the BILR model under setting 2) above (spline function non-increasing in ) simultaneously maximises the LPML criterion and minimises the WAIC, DIC, and RMSE, respectively. According to the results in Table 3, we select the model with six constrained splines. Thus, our final model has a spline function for server :

(13)

with six constrained splines, that is, for all servers .

5 Results

In this Section, we report results of the model fitted to point-by-point data for main-draw singles Grand Slam matches from 2012 forward, which were described in Section 2. We divide both the male and female datasets into training and test sets. In both training sets we have randomly chosen servers, while the receivers are in the male training set and in the female training set. Conversely, in the male test set we have servers and receivers, whereas in the female test set we have servers and receivers. We fit model (1)-(10) separately on both male and female training sets, and perform predictions on the hold-out test sets. Our aim is to predict the conditional probability of winning a point for servers in the test sets by borrowing information from the training set results.

The posterior update of the model parameters was performed via Gibbs sampling, implemented by the rjags package (Plummer et al., 2016) in the R programming language (R Core Team, 2013). Posterior summaries were based on

draws from the posteriors, with a burn-in of 1000 iterations and thinning every 20 iterations to reduce the autocorrelation in the posterior samples. Convergence of the Markov Chain ha been assessed by visual inspection ad using the

coda package (Plummer et al., 2006). The sampler appeared to converge rapidly and mix efficiently.

Summaries of the serve advantage model parameters suggest a strong serve advantage. The probability of the server winning conditional on the point ending on serve, , is with credible interval (C.I.) for men and for women with C.I. , respectively. When the serve advantage is lost, e.g. , the probability of winning a point is mainly given by the rally ability of the server against the rally ability of the opponent. In this case, has credible intervals for men and for women, respectively.

Figure 5: Probability of winning a point as a function of rally length for Rafael Nadal and Roger Federer. The points represents the real data, while the line is the posterior mean estimate of the probability of winning a point as a function of rally length obtained with the model. The blue dashed lines are the credible intervals.

In Figure 5 we show the estimated probability of winning a point as a function of rally length for Rafael Nadal and Roger Federer. The mean posterior curves are very similar: both players have the highest chance of winning the point on serve, and then this probability decreases. It is evident that the posterior mean estimate of the probability of winning a point on serve does not undergo an exponential decay, and a similar pattern for this estimate is observed on other players as well. We remark again that the curve estimates a global trend, namely it describes how the serve advantage drops with rally length. Nevertheless, the value , for , is the estimate of the serve advantage for athlete at the (discrete) value of rally length, . In our Figure we decide to plot the posterior estimate of as a continuous trajectory to underline the longitudinal structure of the data.

Figure 6: The posterior median rally ability estimated under the baseline model (Section 3.3), on the -axis, against the total serve advantage (Eq. (3)), on the -axis, for male players in the ATP tournaments. The red lines represent the median rally ability of male athletes, parallel to the -axis, and the median of the serve advantage. The red points indicate the top four players of the ATP tournaments. We only display those athletes whose 95% CI for serve advantage and rally ability do not include zero.

Figure 6 displays the estimated posterior median serve advantage versus the estimated posterior median rally ability. Specifically, the -axis displays the posterior median of estimated under the baseline model (Section 3.3), whereas the -axis displays the total serve advantage , where is defined in Eq. (3). Let us observe the three top players according to the ATP singles ranking as of January 2019, namely, Roger Federer, Rafael Nadal and Novak Djokovic. We notice that Djokovic excels in terms of rally ability, confirming that he is better in defence than in attack. Conversely, Federer wins more at the first shot than on the long play. Finally, Nadal stands out in both terms of serve advantage and rally ability. Figure 7 displays the same plot for the female dataset. We observe that Serena Williams excels on the long play, while Angelique Kerber and Simona Halep are better on serve. Caroline Wozniacki displays a good balance between serve and rally abilities.

Figure 7: The posterior median rally ability estimated under the baseline model (Section 3.3), on the -axis, against the total serve advantage (Eq. (3)), on the -axis, for female players in the WTA tournaments. The red lines represent the median rally ability of female athletes, parallel to the -axis, and the median of the serve advantage. The red points indicate the top four players of the WTA tournaments. We only display those athletes whose 95% CI for serve advantage and rally ability do not include zero.

Further, we want to investigate the effect of the surface on the rally ability. To this end we fit the extension of our model described in Section 3.4. Since the court is likely to have an effect on the player’s rally abilities, we study how the court affects the players’ skills. We report here the posterior median estimate for the rally ability along with 95% credible intervals for the three different courts for the best players in the ATP and WTA tournaments, respectively. We also compute the posterior median estimate, with 95% credible intervals, for , obtained with the model which does not take the court effect into account (Equation (2)). These estimates, reported in Appendix B, are used to rank the athletes and understand how the different courts impact to the game of these top players.

Type of court
Ranking Baseline Clay Grass Hard
ATP
1 Djokovic (0.35) Nadal (0.52) Federer (0.28) Djokovic (0.34)
2 Nadal (0.31) Djokovic (0.29) Djokovic (0.28) Nadal (0.28)
3 Federer (0.16) Federer (0.09) Nadal (0.17) Federer (0.14)
WTA
3 Wozniacki (0.17) Halep (0.20) Kerber (0.16) Wozniacki (0.21)
1 Halep (0.11) Wozniacki (0.09) Halep (0.11) Kerber (0.15)
2 Kerber (0.03) Kerber (0.01) Wozniacki (0.11) Halep (0.11)
Table 4: Ranking of the players on different court surfaces. The second column lists the official ATP and WTA year-end final rankings (by points) for singles for the 2018 championships season.

The results (Table 4) confirm common knowledge about these athletes. Djokovic and Nadal are both great at rallying. Djokovic is good on all courts, while Nadal is very good on clay and hard courts, but less favorite on grass. Federer appears to be weaker in rallying compared to the other two athletes, though he is the strongest on grass courts. Regarding the WTA tournament, Angelique Kerber is good at rallying on both hard and grass court, but underperforming on clay courts. Caroline Wozniacki is good on all types of court, and in fact she is the player with the highest estimated rally ability among the three female athletes. Simona Halep is good on clay, but does not outperform other players either on grass and hard courts. In general, however, the female athlete with the highest estimated in the WTA dataset is Serena Williams (Figure 7). The median estimate of her rally ability in the baseline model is 0.26 (0.18-0.35), whereas the estimates on clay, grass, and hard courts are, respectively, 0.17 (0.05-0.30), 0.17 (0.11-0.21) and 0.27 (0.17-0.37).

Figure 8: Probability of winning a point as a function of rally length for Gilles Simon, on the left, and for Eugenie Bouchard, on the right. The points represent the real data, while the line is the estimated posterior mean probability of winning a point as a function of rally length obtained with the model. The blue dashed lines are the credible intervals.

In Figure 8 we observe the out-of-sample prediction for two players belonging to the male and female test sets, Gilles Simon and Eugenie Bouchard. The estimated probability of winning a point for a server in the test set (e.g., the black solid curve in Figure 8), is obtained by drawing the basis functions coefficients and the positive random decrements as per Equations (6), (9), (8), using the posterior estimates of the non-subject specific parameters, that is, , , and . The rally ability is just computed in the training phase. The strength of the hierarchical model is the ability to infer the conditional probability of winning for a hold-out subject by borrowing strength from athletes in the training dataset. The estimated trajectory for these players is in line with the observed realisations given by the points in Figure 8.

6 Conclusions

In this paper, we presented a framework to modelling the serve advantage in elite tennis. Our approach extends Kovalchik (2018b) by replacing a simple decay exponential function for the serve advantage with a B-spline basis function decomposition, thus achieving more flexible results. Constraints on the basis function coefficients guarantee that the serve advantage is non-increasing with rally length. As in Kovalchik (2018b), we allow the conditional probability of winning on serve to also depend on the rally ability of the two players, and investigate how the different types of court may impact on such rally ability. When the exponential decay function in Kovalchik (2018b) goes to zero, the conditional probability of winning a point is only given by the difference between two rally abilities, thus a constant. Conversely, our spline function is defined on by construction, and therefore non-zero everywhere in the spline domain. This results in higher uncertainty as represented by wider credible intervals for large values of rally length. This should be considered as a positive feature of our model, that is able to reflect larger uncertainty in presence of sparser data.

Our results show a sort of trade-off between serve advantage and rally ability. The most successful tennis players in the dataset show higher rally ability (rally ability above training median value) relative to their serve advantage. Indeed, if two players have the same chance of winning the point on the first shot, the match will be won by the player with the higher rally ability. We can conclude that although the service is important, what makes a tennis player great is his/her rally ability.

Although motivated by the analysis of tennis data, our methodology can be applied to pair-comparison data in general, with applications ranging from experimental psychology to the analysis of sports tournaments to genetics.

Appendix A Proof of Proposition 1

Consider the restricted spline function

where is the order of the B-splines. Following Formula (12) on page 116 and Formula (13) of de Boor (2001), we can compute its derivative as

where we used the fact that that . The latter follows from Formula (37) on page 96 of de Boor (2001), and from the fact that for and . It is straightforward to observe that the constant is smaller or equal than zero if and only if

The latter property is equivalent to requiring that the restricted control polygon is not increasing on its support, i.e. for .

Appendix B Rally abilities on different types of courts

Type of court
Players Baseline Clay Grass Hard
Novak Djokovic 0.35 (0.23-0.46) 0.26 (0.17-0.36) 0.25 (0.14-0.36) 0.30 (0.22-0.40)
Rafael Nadal 0.31 (0.19-0.40) 0.49 (0.39-0.60) 0.15 (0.03-0.28) 0.25 (0.16-0.32)
Roger Federer 0.16 (0.07-0.25) 0.09 (0.01-0.21) 0.27 (0.18-0.35) 0.11 (0.03-0.20)
Caroline Wozniacki 0.17 (0.03-0.29) 0.03 (0.00-0.11) 0.11 (0.04-0.26) 0.21 (0.11-0.31)
Simona Halep 0.11 (0.03-0.21) 0.20 (0.09-0.30) 0.11 (0.01-0.22) 0.11 (0.02-0.19)
Angelique Kerber 0.03 (0.01-0.10) 0.09 (0.03-0.13) 0.16 (0.06-0.32) 0.15 (0.06-0.24)
Table 5: Credible intervals for the rally abilities of the top players for the ATP and the WTA tournaments.

References

  • Abraham and Khadraoui (2015) Abraham, C. and Khadraoui, K. (2015). Bayesian regression with b-splines under combinations of shape constraints and smoothness properties. Statistica Neerlandica, 69(2):150–170.
  • Baio and Blangiardo (2010) Baio, G. and Blangiardo, M. (2010). Bayesian hierarchical model for the prediction of football results. Journal of Applied Statistics, 69:150–170.
  • Barnett and Clarke (2005) Barnett, T. and Clarke, S. R. (2005). Combining player statistics to predict outcomes of tennis matches. IMA Journal of Management Mathematics, 16:113–120.
  • Bradley and Terry (1952) Bradley, R. A. and Terry, M. E. (1952). Rank analysis of incomplete block designs: the method of paired comparisons. Biometrika, 39:324–345.
  • de Boor (2001) de Boor, C. (2001). A practical guide to splines. Springer-Verlag.
  • Gelfand and Sahu (1999) Gelfand, A. E. and Sahu, S. K. (1999). Identifiability, improper priors, and gibbs sampling for generalized linear models. Journal of the American Statistical Association, 94(445):247–253.
  • Geyser and Eddy (1979) Geyser, S. and Eddy, F. (1979). A predictive approach to model selection. Journal of the American Statistical Association, 3:153–160.
  • Gilsdorf and Sukhatme (2008) Gilsdorf, K. F. and Sukhatme, V. A. (2008). Testing rosen’s sequential elimination tournament model: Incentives and player performance in professional tennis. Journal of Sport Economics, 9:287—303.
  • Glickman (1999) Glickman, M. E. (1999). Parameter estimation in large dynamic paired comparison experiments. Journal of the Royal Statistical Society: Series C (Applied Statistics), 48:377—394.
  • Gomes et al. (2011) Gomes, R. V., Coutts, A. J., Viveiros, L., and Aoki, M. S. (2011). Physiological demands of match-play in elite tennis: A case study. European Journal of Sport Science, 11(2):105–109.
  • Gorgi et al. (2018) Gorgi, P., Koopman, S. J. S., and Lit, R. (2018). The analysis and forecasting of ATP tennis matches using a high-dimensional dynamic model. Tinbergen Institute Discussion Papers 18-009/III, Tinbergen Institute.
  • Irons et al. (2014) Irons, D. J., Buckley, S., and Paulden, T. (2014). Developing an improved tennis ranking system. Journal of Quantitative Analysis in Sports, 10(2).
  • Klaassen and Magnus (2003) Klaassen, F. J. and Magnus, J. R. (2003). Forecasting the winner of a tennis match. European Journal of Operational Research, 148:257 – 267.
  • Kotze J. and Rothberg (2000) Kotze J., Mitchell, S. J. and Rothberg, S. (2000). The role of the racket in high-speed tennis serves. Sports Engineering, 3:67–84.
  • Kovalchik (2016) Kovalchik, S. (2016). Searching for the goat of tennis win prediction. Journal of Quantitative Analysis in Sports, 12:127–138.
  • Kovalchik (2018a) Kovalchik, S. (2018a). deuce: Resources for Analysis of Professional Tennis Data. R package version 1.2.
  • Kovalchik (2018b) Kovalchik, S. (2018b). The serve advantage in professional tennis. Under Review.
  • Lees (2003) Lees, A. (2003). Science and the major racket sports: a review. Journal of Sports Sciences, 2:707–732.
  • Newton and Keller (2005) Newton, P. K. and Keller, J. B. (2005). Probability of winning at tennis i. theory and data. Studies in Applied Mathematics, 114(3):241–269.
  • O’Donoghue and Brown (2008) O’Donoghue, G. P. and Brown, E. (2008). The importance of service in grand slam singles tennis. International Journal of Performance Analysis in Sport, 8(3):70–78.
  • Orani (2019) Orani, V. (2019). Bayesian isotonic logistic regression via constrained splines: an application to estimating the serve advantage in professional tennis. Master’s thesis, University of Torino.
  • Plummer et al. (2006) Plummer, M., Best, N., Cowles, K., and Vines, K. (2006). Coda: Convergence diagnosis and output analysis for mcmc. R News, 6(1):7–11.
  • Plummer et al. (2016) Plummer, M., Stukalov, A., and Denwood, M. J. (2016). rjags: Bayesian Graphical Models Using MCMC.
  • R Core Team (2013) R Core Team (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Smith (2013) Smith, S. (2013). Effects of playing surface and gender on rally durations in singles grand slam tennis. PhD thesis, Cardiff Metropolitan University.
  • Spiegelhalter and der Linde A. (2002) Spiegelhalter, D. J. Best, N. G. and der Linde A., V. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, 64:583–639.
  • Starbuck et al. (2016) Starbuck, C., Damm, L., Clarke, J., Carré, M., Capel-Davis, J., Miller, S., Stiles, V., and Dixon, S. (2016). The influence of tennis court surfaces on player perceptions and biomechanical response. Journal of Sports Sciences, 34(17):1627–1636. PMID: 26699792.
  • Watanabe (2010) Watanabe, S. (2010). Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.

    Journal of Machine Learning Research

    , 11:3571–3594.