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 playerspecific 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 BradleyTerry 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 Bsplines basis functions, with athletespecific 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 Bsplines, while their shape is controlled by the associated control polygon (de Boor, 2001). In particular, to ensure the serve advantage is nonincreasing with rally length we constrain the spline function to be nonincreasing 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 outofsample 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 pointbypoint data for maindraw 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 Sackmann^{1}^{1}1https://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 
Short rallies  Long rallies  Total  

ATP  130577  14933  145510 
WTA  71592  10288  81880 
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 oddrally lengths compared to the evenrally lengths. This pattern can be explained by the fact that evennumbered 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 evenshots 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 oddnumbered rallies are the winners or errors made by the receiver. In particular, the oddshots in Figure 1 represent the errors done by the receivers.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.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 athletespecific 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 nonlinear 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 Bsplines basis functions of order on , and the ’s are athletespecific basis function coefficients. Specifically, let and be a nondecreasing sequence of knots such that for all , and and . Function is a linear combination of the Bsplines , 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) logodds 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 Bsplines 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 subinterval of its domain , for example if is monotone decreasing in . To simplify the notation, we will omit index that denotes individualspecific 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 nonincreasing in . While the nonincreasing 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 nonincreasing 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 nonincreasing 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 nondecreasing 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.To ensure the serve advantage is nonincreasing 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 nonincreasing on by imposing that the associated restricted control polygon is nonincreasing 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 nonincreasing on its support , then the restricted spline function is also nonincreasing 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 nonincreasing constraint on , the broken line with vertexes is nonincreasing 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 nonincreasing 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 athletespecific 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 Bspline 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 nonincreasing 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
:(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 nonincreasing, 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. Nonidentifiability 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 sumtozero 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 noninformative 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 BradlyTerry 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 subjectspecific 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 noninformative priors on the hyperparameters:
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 paircomparison 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 Bsplines with support in , thus the resulting spline function is nonincreasing in (partially monotone); and 3) spline function constrained to be nonincreasing 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 crossvalidation 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 Bspline 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 subinterval of the spline function domain where monotonicity is to be imposed. For setting 2) above, this subinterval 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 Bspline 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 datarich 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 nonincreasing (bottom right) display a nonincreasing behaviour for large values of rally length.For a quantitative evaluation of the performance of the four approaches, we compute the goodnessoffit measures for these models, reported in Table 3.
GoodnessoffitMeasures  

LPML  WAIC  DIC  RMSE  
Exponential model  52813.7  105821.3  105876  22.52 
Spline models  
No constraints  52760.4  105853.9  105818  20.71 
Nonincreasing in  52739.1  105747.6  105828  19.32 
Nonincreasing in  52744.9  105790.4  105875  20.37 
Although no dramatic difference in performance emerge, the BILR model under setting 2) above (spline function nonincreasing 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 pointbypoint data for maindraw 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 holdout 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 burnin 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.
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 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.
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) 
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.180.35), whereas the estimates on clay, grass, and hard courts are, respectively, 0.17 (0.050.30), 0.17 (0.110.21) and 0.27 (0.170.37).
In Figure 8 we observe the outofsample 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 nonsubject 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 holdout 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 Bspline basis function decomposition, thus achieving more flexible results. Constraints on the basis function coefficients guarantee that the serve advantage is nonincreasing 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 nonzero 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 tradeoff 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 paircomparison 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 Bsplines. 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.230.46)  0.26 (0.170.36)  0.25 (0.140.36)  0.30 (0.220.40) 
Rafael Nadal  0.31 (0.190.40)  0.49 (0.390.60)  0.15 (0.030.28)  0.25 (0.160.32) 
Roger Federer  0.16 (0.070.25)  0.09 (0.010.21)  0.27 (0.180.35)  0.11 (0.030.20) 
Caroline Wozniacki  0.17 (0.030.29)  0.03 (0.000.11)  0.11 (0.040.26)  0.21 (0.110.31) 
Simona Halep  0.11 (0.030.21)  0.20 (0.090.30)  0.11 (0.010.22)  0.11 (0.020.19) 
Angelique Kerber  0.03 (0.010.10)  0.09 (0.030.13)  0.16 (0.060.32)  0.15 (0.060.24) 
References
 Abraham and Khadraoui (2015) Abraham, C. and Khadraoui, K. (2015). Bayesian regression with bsplines 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. SpringerVerlag.
 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 matchplay 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 highdimensional dynamic model. Tinbergen Institute Discussion Papers 18009/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 highspeed 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., CapelDavis, 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.
Comments
There are no comments yet.