Spatial modeling of shot conversion in soccer to single out goalscoring ability

4.1 Proposed model

We start with a probit regression model equipped with a spatially correlated error process, in light of the discussion in the previous section. For shots i = 1, …, N, suppose Y _i is the binary outcome (which takes the value 1 if the ith shot was converted to a goal), s _i denotes the location of the shot and X _i denotes the column vector with the values of the covariates.

Our proposed model starts with the assumption that we have a latent variable r _i , such that Y _i  ∼ ber (p _i ), p _i  = P (r _i  > 0). r _i depends on the values of the covariates following a simple linear regression model, and z_m(i) is an additive random effect that quantifies the player’s ability. Here, m (i) corresponds to the player who took the ith shot. In order to avoid the issue of over-fitting, we separately consider the players whose number of shots is more than a minimum cutoff, say s _m . Throughout this study, we use s _m  = 10. For all other players, m (i) is considered to be same. If we use Pl₁, …, Pl_M-1 to denote the players whose number of shots is greater than s _m , and we associate z _j to Pl _j , then

m ( i ) = { j if i th shot is taken by Pl j M if i th shot is not taken by Pl 1 , … , Pl M - 1 .(4.1)

We further assume that, given the values of the covariates, the outcomes will be dependent of each other through the shot locations and this dependence would decrease with distance. These assumptions support our findings while exploring the data in the previous section.

We are going to consider a hierarchical structure. First of all, Y i = I ( r i > 0 ) and r i = X i ′ θ + z m ( i ) ε j ,(4.2) where θ is a parameter vector of appropriate order corresponding to the covariates, z_m(i) is the player-effect on the shot, and _epsi is the error process. The term X _i ′θ is considered to be an additive combination of the effects of some continuous and some discrete covariates. z_m(i)’s are considered to be independent and identically distributed (iid) random effects, coming from N ( 0 , σ p 2 ) distribution. Note that there will be M different random effects, if there are M - 1 different players whose number of shots is greater than s _m . We further assume that the error process has the following structure: ε j = w i + e i .(4.3)

Here w _i denotes a zero-mean spatially correlated process and e _i stands for a zero-mean independent and identically distributed (iid) white noise process. We assume that e _i ’s are iid N (0, σ²) while for w _i ’s, we take a correlation structure that decays exponentially. In particular, the covariance between w _i and w _j (i.e. for locations s _i and s _j ) is Cov ( w i , w j ) = σ w 2 exp ( - φ ∥ s i - s j ∥ ) ,(4.4) where the distance function norm ∥s _i  - s _j ∥ is taken as the Euclidean distance of the two locations.

Now, let us use Y to denote the vector of all outcomes, r and p to denote the vector of r _i ’s and p _i ’s, respectively. X is the design matrix such that X′ = [X₁ : … : X _N ]. We will denote the covariance matrix of w by σ w 2 Σ w . Finally, let z be a vector of length M, corresponding to the player-effects. Then, the vector of z_m(i)’s can be written in the form Az, where each row of A (a N × M matrix) denotes the index of the player effect. Thus, we can write the full model in the following form: Y i = I ( r i > 0 ) ,(4.5) r = X θ + Az + w + e , such that z ∼ N ( 0 , σ p 2 I M × M ) , w ∼ N ( 0 , σ w 2 Σ w ) , and e ∼ N ( 0 , σ 2 I N × N ) .

We are going to use a Bayesian framework to get estimates for the parameters in the model. For each of the components in the parameter vector θ, we will consider improper Jeffrey’s prior distribution. On the other hand, we will assume the two variance components σ 2 , σ w 2 to be equal. While we acknowledge the simplicity associated with this assumption, it is worth mentioning that we have performed some simulation studies and found that the performance is approximately the same with the assumption. The computational complexity becomes cumbersome already with complicated priors in Gibbs sampling and the computational burden decreases manifold if we consider the variances equal. Also, we achieve satisfactory prediction power with this parsimonious assumption and hence, don’t want to increase our number of parameters. Further discussion on this assumption is provided at the end of Section 4.2.

The priors we associate with σ² (same as σ w 2 ) and σ p 2 are independent inverse gamma (IG) distribution with parameters a, b. We would take a > 1 to avoid improper posterior distributions.

And finally, the parameter φ (refer to Equation (4.4)) will be kept fixed throughout the analysis. In order to find out the best possible value for φ, we would consider a cross-validation scheme. The possible choices for φ used in the study were between 0.05 and 1. Based on the relationship e^-φd ≈ 0.05, we chose the above range, which corresponds to actual distances in the field ranging from 3 yards to 60 yards. The validation scheme would consider prediction for 20% of all the shots after estimating the parameters from the other 80% and would find out the mean squared error for those predictions. For shots considered for the validation purpose (let us denote them by i ∈ V = {v₁, …, v _m }), we would predict the probability of them being converted to goals and let the predicted probabilities be denoted by p ˆ i while Y _i ’s are the actual outcomes. Then, for each φ, we calculate the mean squared error of the predictions using the following formula: MSE φ = 1 m ∑ i ∈ V ( Y i - p ˆ i ) 2 .(4.6)

The value of φ with the least mean squared validation error will be our optimal choice for the analysis. The prediction procedure is discussed in detail in Section 4.3.

Throughout the discussion below, K will always indicate a constant term, that may vary from time to time. Using the full model described in Section 4.1 and the prior distributions described in the paragraph following that equation, we can write the joint posterior distribution as:

log π ( r , θ , σ 2 , w , z Y ) = K + ∑ i = 1 N { Y i log ( P ( r i > 0 ) ) + ( 1 - Y i ) log ( P ( r i ≤ 0 ) ) } - w ′ Σ w - 1 w 2 σ 2 - ∥ r - X θ - Az - w 2 ∥ 2 σ 2 - ( a + N + 1 ) log σ 2 - b σ 2 - ∥ z 2 ∥ 2 σ p 2 - ( a + M 2 + 1 ) log σ p 2 - b σ p 2 .(4.7)

Since simulating directly from the joint posterior distribution is difficult, we would use the principles of Gibbs sampling, which sequentially updates every parameter in an iterative way until convergence. For σ², we can get the full conditional distribution as: log π ( σ 2 Y , r , θ , w , z ) = K - ( a + N + 1 ) log σ 2 - 1 σ 2 [ b + 1 2 ∥ r - X θ - Az - w ∥ 2 + 1 2 w ′ Σ w - 1 w ] and thus,

( σ 2 Y , r , θ , w , z ) ∼ IG ( a + N , b + 1 2 ∥ r - X θ - Az - w ∥ 2 + 1 2 w ′ Σ w - 1 w ) .(4.8)

Similarly, we get that ( σ p 2 ∣ Y , r , θ , w , z ) ∼ IG ( a + M 2 , b + 1 2 ∑ k = 1 m z k 2 ) .(4.9)

For θ, we get that log π ( θ ∣ Y , r , σ 2 , w , z ) = K - 1 2 θ ′ ( X ′ X σ 2 ) θ + 1 2 · 2 θ ′ · X ′ ( r - Az - w ) σ 2 . Straightforward calculations then tell us that the conditional posterior of θ is

( θ ∣ Y , r , σ 2 , w , z ) ∼ N ( ( X ′ X ) - 1 X ′ ( r - Az - w ) , σ 2 ( X ′ X ) - 1 ) .(4.10)

On the other hand, using similar techniques, we can get the posterior distribution for w as: ( w ∣ Y , r , θ , σ 2 , z ) ∼ N ( ( I + Σ w - 1 ) - 1 ( r - X θ - Az ) , σ 2 ( I + Σ w - 1 ) - 1 ) .(4.11)

Further, if we define ∑ i = 1 N I { m ( i ) = k } = n k , we can get the posterior distribution for z _k as:

( z k ∣ Y , r , θ , σ 2 , w ) ∼ N ( ( n k + σ 2 σ p 2 ) - 1 ∑ i : m ( i ) = k ( r i - X i ′ θ - w i ) , ( n k σ 2 + 1 σ p 2 ) - 1 ) .(4.12)

Note that the means of the conditional distributions of both θ and w require the computation of the products of a matrix and a vector, but the dispersion matrices of the same are updated at every step only through the value of σ and so, one needs to compute the inverses of three big matrices, namely X′X, Σ _w and ( I + Σ w - 1 ) , only once throughout the MCMC analysis. Hence, directly using the above conditional distributions at every step will not suffer from huge computational burden. However, one can always use component-wise conditional posterior distributions to further increase the efficiency.

Finally, we find the conditional distribution of r. The joint posterior distributions above clearly show that conditional on data and other parameters, each component of r is independent of others and simple calculations would reveal that the conditional distributions are truncated normal as follows. ( r i ∣ Y , θ , w , σ 2 , z ) ∼ { TN ( α i , σ 2 ; 0 , ∞ ) if Y i = 1 TN ( α i , σ 2 ; - ∞ , 0 ) if Y i = 0(4.13) In the above, α _i is ith component of α = Xθ + w + Az and TN (μ, τ² ; a, b) denotes a normal distribution with mean μ and variance τ², truncated in the interval [a, b].

Now to note that, if we assume σ w 2 ≠ σ 2 , then the posterior distribution of w (see Equation (4.11)) will become

( w ∣ Y , r , θ , σ 2 , σ w 2 , z ) ∼ N ( σ 2 ( I + σ 2 σ w 2 Σ w - 1 ) - 1 ( r - X θ - Az ) , σ 2 ( I + σ 2 σ w 2 Σ w - 1 ) - 1 ) .(4.14)

This implies that at each step of Gibbs sampling, we will have to invert a n × n matrix ( I + σ 2 σ w 2 Σ w - 1 ) , which requires O (n³) time. Naturally, it is computationally not feasible to do it for thousands of iterations when n is large (in our case n = 3957). Whereas, if we follow Equation (4.11), we will have to invert the matrix ( I + Σ w - 1 ) only once and we can use that in subsequent iterations. Hence, this parsimonious assumption leads to a huge decrease in computational complexity and we do not lose much in prediction.

In order to make a new prediction for a shot from location s′, we will take resort to the posterior estimates obtained from the Gibbs sampler. Let us denote the outcome of the shot from s′ by Y (s′) (X (s′), r (s′), z (s′) and w (s′) are defined accordingly). Note that the posterior predictive distribution f (Y (s′) ∣Y) can be written using an integral for the product of the density functions, where the integral is taken with respect to all the parameter vectors. However, instead of solving that integral, a better and more convenient idea is to use the Gibbs sampler estimates to draw observations from the posterior predictive distribution and one can do it sequentially.

At first, we draw samples for r, w, θ, σ² using the conditional posterior distributions described in Section 4.2. Also, if it comes from kth player used in our training data, we define z (s′) as the posterior mean of z _k , else, we simulate z (s′) from N ( 0 , σ p 2 ) . After that, we draw samples for w (s′) using the conditional distribution of (w (s′) ∣w, σ²) as discussed below. Using the generated observations, an estimate for r (s′) can be obtained by taking a realization from a normal distribution with mean X (s′) ∣θ + w (s′) + z (s′) and variance σ². Finally, we set Y (s′) = ind (r (s′) >0). Alternately, one can also use Y (s′) ∼ ber (p (s′)), where p (s′) = P (r (s′) >0) can be computed using the sampled values of θ, σ², z (s′) and w (s′).

To compute the conditional distribution of (w (s′) ∣w, σ²), we start with the fact that I ( w ( s ′ ) w ) ∼ N ( 0 , σ 2 [ 1 Σ ′ s - s ′ Σ s - s ′ Σ w ] ) .

Here, Σ_s-s′ is the column vector denoting the covariance of w (s′) with the elements of w. Following Equation (4.4), the jth element of Σ_s-s′ is exp(- φ∥s _j  - s′∥). Using the principle of conditional distribution for multivariate normal distribution, we can then say that w ( s ′ ) w ∼ N ( μ new , σ new 2 ) ,(4.15) where σ new 2 = σ 2 ( 1 - ( Σ s − s ' ' Σ w - 1 Σ s - s ′ ) ) , μ new = Σ s − s ' ' Σ w - 1 w .

Based on our model, we propose two measures to quantify a player’s ability on the pitch.

To start with, we would like to introduce, what we call, the Shooting Prowess (SP) of a player. We define this measure SP _k as the posterior mean of z _k i.e. the player-specific random effect. Recall that we use a cutoff of s _m  = 10 matches to evaluate this random effect, and hence, discrete cases of players taking few shots to score will be considered as general effect in our model. This would effectively help us in identifying the players who are in general good shooters.

Next, we would like to quantify a player’s Positioning sense (PS), which is a pivotal thing for a striker on the field in terms of scoring goals. One has to be at the right place at the right time to get a goal in his name. From our training model, we can estimate the shot conversion probability of an average player from a particular distance, angle and type of opportunity. It is a different issue whether the player has grabbed the opportunity or not, but getting into a perfect position, he could ensure that a shot was fired.

In order to do that for the kth player, suppose that in a total of g _k matches, he has taken n _k shots (calculated over a time period, possibly over a whole season or for the entire career) in total. The expected probability of the ith shot, i = 1, …, n _k , is denoted by p ˆ i . Then, summing up the predicted probability of shot conversion for that player can give us a perfect idea about his positioning sense. Necessarily, the PS _k (Positioning Sense) of the kth player is defined by PS k = 1 g k ∑ i = 1 n k p ˆ i .(4.16)

The fact that we calculate the average over the number of matches and not over the set of all shots captures the information of whether a player takes more shots in general.

5.1 Fitting the model

Following the discussions in the previous sections, we consider two different models for ‘other shots’ and ‘headers’, and they can be written as follows.

For other shots, let Y _i  = I   (r _i  > 0) denote whether ith shot was a goal. Then, we set r _i  = μ _i  + z_m(i) + w _i  + ϵ _i , where z_m(i), w _i , ϵ _i are as described in the previous section and μ _i is the mean function defined by

μ i = β 0 + β 1 log ( D i ) + β 2 cos ( A i ) + β 3 K i + β 4 P i + β 5 I H i = 1 } + β 6 I { F i = 1 } + ∑ j = 0 1 β 7 , j I { GD i = j } + β 8 I { S i = 1 } + ∑ k = 2 3 β 9 , j I { BP i = j } .(5.1)

Here, D _i , A _i , K _i and P _i denote the continuous covariates, namely distance, angle, keepreach and proportion of shots converted against the same opponents. Indicators for H _i , F _i and S _i equal to 1 are used to denote if it was a home game, if it was in the first half and whether it happened during stoppage time. GD _i denotes whether the goal difference at the time of the shot was zero (j = 0) or positive (j = 1). And finally, BP _i shows which body part (j = 1, 2, 3 corresponding to other, left foot and right foot) was used in taking the shot. Note that we take β_7,-1 = β_9,1 = 0 for identifiability purpose.

On the other hand, for ‘headers’, we do not consider the covariate BP, but the other model remains unchanged.

Now, using the mean function (Section 5.1), we fit the model described in Section 4. First of all, using the procedure described towards the end of Section 4.1, we evaluate the best choice of φ to be used in the spatial covariance function (Equation (4.4)) It is obtained to be 0.75 (and 0.45) for ‘headers’ (and for ‘other shots). Following the approximate relation e^-φd ≈ 0.05, we can say that the spatial correlation in converting headers (and other shots) is insignificant for locations approximately 4 yards (and 6.7 yards) or more apart from each other. Then, we use the Gibbs sampler to estimate the parameters. We have used 10000 runs in the burn-in period and have used MCMC samples of size 1000 to evaluate the estimates, standard errors and the confidence intervals of the parameters. Summary of the results for the ‘headers’ are shown in Table 3, while the same for ‘other shots’ are displayed in Table 4.

From the results, we can see that the effect of logarithm of distance is significantly negative for both headers and other shots. This matches with the established theory that the closer you shoot from, the greater chances of the shot being converted to goal. However, the conclusion about the ‘angle’ of shots is interesting. In the simple logistic regression model used by current researchers in this field, it was found that a smaller angle (higher value of the cosine of the angles) significantly increases the goal-scoring chance of any shot. But, in our model, we have found that the angle of the shot is not a significant covariate for headers. This happens probably because the effect of the angle is already being explained by the spatial process. On the contrary, even with the spatial dependence in the model, the angle has significant positive effect on the conversion rate for other shots.

So far as the other variables are concerned, only keeper’s reach and stoppage time are significant for headers. The chance of converting from a header increases when the keeper’s reach is less. Also, headers taken during the stoppage time has a significantly higher chance of being converted. On the other hand, we also find that the two variance components are very different. While the spatially correlated process shows a variance of 0.999, the player effect explains more variability (estimated value of σ p 2 is 1.656) in the data.

Further, for non-headed shots, barring the distance and angle of the shot, no other situational or match-specific covariate is significant. One interesting observation is the estimates and confidence intervals for the effect of left or right footed shots. They are very close to each other, thereby supporting the nature of symmetry in the game of soccer. The variability in the data, similar to the previous case, is explained more by the player-effect ( σ p 2 is estimated to be 1.675) while the spatially correlated process has a variance of 0.995. These values are similar to the ‘headers’ data.

We next delve into the diagnostics of the model, and compare the performance of our model to three different models. First, we consider the simple logistic regression model (SLRM) which is the most common one to use for such data. Second, we consider the spatial probit model proposed by Klier and McMillen (2008). It is in essence similar to our proposed model. This model, henceforth denoted as KMM, uses a spatially related weight matrix, where the weights are corresponding to the number of observations contiguous to every observation. The parameters, including the spatial interaction coefficient, are estimated by a linearized version of generalized method of moments estimators for probit models. We include this model in this comparison study to understand if our specification of the spatial correlation can be deemed appropriate. Finally, we include a machine-learning algorithm. A type of neural network (later denoted as NN) that performs well on such a binary classification problem from a set of vectors is a simple stack of fully connected (dense) layers with rectified linear unit activation function. We use two dense layers with 16 and 8 hidden units and then pass it onto a layer with sigmoid function as output to train classification probabilities.

In order to compare the performance of the different techniques, at first we calculate the Brier score for all the models. This is recognized as one of the most important score function to measure the accuracy of model fit for binary type data. It is defined as ∑ t = 1 N ( f t - o t ) 2 / N where f _t is the expected probability of forecast and o _t is the actual outcome. For the non-headed shots, we compute it to be 0.067, while for the headers it is approximately 0.061. Note that for both cases, this score is around 30% less than what we observed in the other three candidate models, namely SLRM, KMM or NN.

Another standard method to evaluate probabilistic prediction is the log score, which is defined by ∑ t = 1 N { o t log f t + ( 1 - o t ) log ( 1 - f t ) } . A higher value of the log score will indicate a better predictive ability of the model. We have observed improvement in our model according to this measure as well. For headers, SLRM, KMM and NN yield log score of approximately -301.7, -338.9 and -299.4, respectively, while the same for our model is nearly -184.3. Thus, our model performs more than 38% better in this aspect. Similar phenomena has been observed for non-headed shots as well.

Next, we evaluate the proportion of times the models could classify the outcomes correctly. Here, the improvement with our model is about 20%, as compared to the other three methods.

The above results are summarized in Table 5.

As a final piece, we want to see how accurately we can predict out-of-sample outcomes based on the model we developed. At this point, one should keep in mind that this is a case of rare events binary regression and one of the main issues with these scenarios is that the prediction problem is immensely difficult and some calibration method is always needed, which further requires more information on the whole population, as pointed out by King and Zeng (2001). Since we do not have that information, we have to rely on usual prediction methods (refer to Section 4.3). We use cross-validation techniques to evaluate how good the predictions are, by taking a subset of the full data (approximately 80%) to learn the model and applying it on the rest (approximately 20%). Afterwards, in order to evaluate the prediction efficiency, we follow the general scoring rules outlined by Merkle and Steyvers (2013).

In the aforementioned study, the authors discussed the beta family of proper scoring rules, proposed by Buja et al. (2005) for binary regression problems, in detail. Using different combinations of the parameters (α, β) for the beta distribution, one can choose the most appropriate scoring rule. For example, α = β = 0 corresponds to log score while α = β = 1 corresponds to Brier score. On the other hand, Merkle and Steyvers (2013) presented another way of thinking about the parameters, where one can assume that the value of α/(α + β) (or 1 - α/(α + β)) is the cost of a false positive (or false negative). In light of this, one should note that scoring rules with α < β emphasize low-probability forecasts, in the sense that such a scoring rule would heavily penalize a prediction that attach a low probability to a successful outcome (goals). Similarly, scoring rules with α > β emphasize high-probability forecasts. Using the R package scoring (Merkle (2014)), we have computed the scores for different values of α/(α + β). The logarithm of the scores for the three models are plotted against different cost functions in Fig. 5.

OPEN IN VIEWER

Fig.5

Logarithm of the score function, corresponding to different cost values, for simple logistic regression model (SLRM), Klier-McMillen model (KMM), Neural network (NN) and for our proposed model.

We can see that the performances of the three methods are very similar for other shots. But, our model shows consistently better (higher score function) performance when the cost for a false positive is less than 0.7. From a practical point of view, it is probably more important to predict a goal successfully, and naturally, the cost for a false positive should be less. And we can see that in such cases, our model performs better than the other two. Even for α = β, when the cost is equal for goal and miss, calculations showed that for our model, the score function is approximately 20% higher than the other candidate models.

For headers, though, the results are a bit different. While similar conclusions remain true for our model as compared to SLRM or NN, KMM performs the best. Hence, one can expect that for headers, the specification for the spatial correlation is probably not being best explained by the exponential decay function.

To conclude this section, we once again emphasize that goals are example of rare events (Anderson and Sally (2013)), as only about 10% of the shots get converted to goals on an average. It is nearly impossible to get a very high prediction accuracy without calibration. Since we do not have the necessary information to calibrate the predictions, we should primarily focus on validating the model. And from the above discussions, we can say that in our proposed model performs better than the three different types of methods we considered.

We now proceed to compute the two measures described in Section 4.4. In line with the way we analyze the data, SP and PS are evaluated separately for headers (denoted as SPH and PSH henceforth) and other shots (denoted as SPS and PSS henceforth), in order to identify the players with best heading or shooting ability. Histograms of the measures for all these players are shown in Fig. 6.

OPEN IN VIEWER

Fig.6

Histogram of the shooting prowess (left) and positioning senses (right) for 381 players, calculated based on fitted probability and actual outcomes.

For the SP values in both cases (figures on the left), one can see that there is a huge peak around 0 and that the histogram is slightly skewed to the right, suggesting that with more shots per game, players become more adept in converting them to goals. Similarly for the figures on the right, we can see that most players showed standard positioning senses, while there were some outliers as well.

Next, to examine if our measures can successfully capture the shooting abilities of the players, we identified the top five shooters in terms of SPS and PSS. However, at this point, to avoid potential small-sample issues, we considered players with more than 5 headers or more than 20 shots. The list of the best players in this regard, along with their real life data from the MLS 2016/17 season are presented in table Table 6. In that table, we also present the data for the top players according to SPH and PSH.

We note that Ignacio Piatti, Giovani Dos Santos, Fanendo Adi, Ola Kamara, Bradley Wright-Phillips and Dominic Dwyer were in the list of top ten highest scorers in 2016/17 season and our measures show that the reason behind this is that all have great shooting or heading abilities. Michael Barrios, a winger for FC Dallas that season, is very good with his positioning senses and shooting abilities. Another winger, Kekuta Manneh also showed good shooting prowess and found a spot in the list of top five.

The results for the heading abilities were something intriguing. At first, we like to point out the special cases of Drew Moor, Sebastian Hines and David Horst. All of them are defenders and their heading abilities are supposed to be good. In soccer, it is often observed that the defenders move ahead during a corner or a long free-kick. Naturally, it would be beneficial for the teams to identify which defenders can take good positions during the set-piece situations and convert the chances that come in the way. It is evident that our model can serve that purpose well. Among forwards, Bradley Wright-Phillips, Cyle Larin and Dominic Dwyer showed really good positioning sense and conversion ability when it comes to heading. Moreover, it was found that positioning sense is significantly positively correlated with heading prowess.

It is worth mention here that the training data we use in our model contained information only from a small part of the season and the above measures are always calculated based on that data. Thus, it is interesting that the list of best players according to these measures consists mostly of the players who were actually at the top of the goal-scoring chart after the season ended. This clearly reflects that we have been able to single out the players’ goal-scoring abilities and positioning senses perfectly. In real life, it has huge potential to identify most suitable players in different scenarios.