Spatial modeling of shot conversion in soccer to single out goalscoring ability
Soudeep Deb \thanks{Email: sdeb@uchicago.edu
Debangan Dey \thanks{Email: ddey1@jhu.edu
\noindent Goals are results of pin-point shots and it is a pivotal decision in soccer when, how and where to shoot. The main contribution of this study is two-fold. At first, after showing that there exists high spatial correlation in the data of shots across games, we introduce a spatial process in the error structure to model the probability of conversion from a shot depending on positional and situational covariates. The model is developed using a full Bayesian framework. Secondly, based on the proposed model, we define two new measures that can appropriately quantify the impact of an individual in soccer, by evaluating the positioning senses and shooting abilities of the players. As a practical application, the method is implemented on Major League Soccer data from 2016/17 season. \\ \noindent {Keywords: Soccer analytics, Probit regression, Bayesian modeling, MCMC, Spatial correlation.
main.tex
\documentclass[letterpaper,12pt]{article}

\usepackage[margin=1in]{geometry}
\usepackage[affil-it]{authblk}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{parskip}
\usepackage{graphicx}
\usepackage[round]{natbib} 

\usepackage{color,hyperref,xcolor}
\hypersetup{colorlinks=true,urlcolor=blue, citecolor=purple}
\usepackage{cleveref}

\newtheoremstyle{general}
{3mm} % Space above
{3mm} % Space below
{\it} % Body font
{} % Indent amount
{\bfseries} % Theorem head font
{.} % Punctuation after theorem head
{.5em} % Space after theorem head
{} % Theorem head spec (can be left empty, meaning `normal')


\def\N{\mathbb{N}}
\def\Z{\mathbb{Z}}
\def\Q{\mathbb{Q}}
\def\R{\mathbb{R}}
\def\C{\mathbb{C}}
\def\p{{\mbox{\open P}}}
\def\E{\mathbb{E}}
\def\I{\mathrm{I}}

\newcommand{\abs}[1]{\left\lvert #1 \right\rvert}
\newcommand{\norm}[1]{\left\| #1 \right\|}
\newcommand{\cov}{\mathrm{Cov}}
\newcommand{\rss}{\mathrm{RSS}}
\newcommand{\bic}{\mathrm{BIC}}
\newcommand{\F}{\mathcal{F}}
\renewcommand{\P}{\mathcal{P}}
\newcommand{\eps}{\varepsilon}
\renewcommand{\L}{\mathbb{L}}
\newcommand{\pdmy}{\delta_{i,\psi}}
\newcommand{\pdmp}{\delta_{i,p}}
\newcommand{\given}{\hspace{2pt}\vert\hspace{2pt}}
\newcommand{\x}{\mathbf{X}}
\newcommand{\Xm}{\tilde{X}}
\newcommand{\dash}{^{\prime}}
\newcommand{\ddash}{^{\prime\prime}}
\newcommand{\gap}{\hspace{3pt}}
\newcommand{\bx}{\bar{X}}
\newcommand{\ind}{\mathbb{I}}
\newcommand{\poi}{\mathrm{Poisson} \ }
\newcommand{\ber}{\mathrm{Bernoulli} \ }
\newcommand{\bm}[1]{$\mathbf{ #1}$}
\newcommand{\argmin}{\operatornamewithlimits{argmin}}
\newcommand{\argmax}{\operatornamewithlimits{argmax}}
\newcommand{\one}{\mathbf{I}}
\newcommand{\logit}{\mathrm{logit}}

\theoremstyle{general}
\newtheorem{theorem}{Theorem}
\newtheorem{assumption}{Assumption}

\numberwithin{equation}{section}
\numberwithin{theorem}{section}


\begin{document}

\title{Spatial modeling of shot conversion in soccer to single out goalscoring ability}

\author{Soudeep Deb
  \thanks{Email: \texttt{sdeb@uchicago.edu}; Corresponding author}}
\affil{Department of Statistics,\\ The University of Chicago,\\ 5747 S Ellis Ave,\\ Chicago, IL, USA. 60615 \\ Phone: +1(312)709-0673}

\author{Debangan Dey
  \thanks{Email: \texttt{ddey1@jhu.edu}}}
\affil{Department of Biostatistics, \\ Johns Hopkins Bloomberg School of Public Health,\\ 615 N Wolfe Street,\\ Baltimore, MD, USA. 2120\\ Phone: +1 (917)292-7447}

\date{}

\maketitle

\begin{abstract}
\noindent Goals are results of pin-point shots and it is a pivotal decision in soccer when, how and where to shoot. The main contribution of this study is two-fold. At first, after showing that there exists high spatial correlation in the data of shots across games, we introduce a spatial process in the error structure to model the probability of conversion from a shot depending on positional and situational covariates. The model is developed using a full Bayesian framework. Secondly, based on the proposed model, we define two new measures that can appropriately quantify the impact of an individual in soccer, by evaluating the positioning senses and shooting abilities of the players. As a practical application, the method is implemented on Major League Soccer data from 2016/17 season.  \\

\noindent {\bf Keywords:} Soccer analytics, Probit regression, Bayesian modeling, MCMC, Spatial correlation.
\end{abstract}




\newpage

\section{Introduction}
\label{sec:introduction}

Association football or Soccer is arguably the most popular sport in the world. Along with the excitement involved with the international fixtures, millions of viewers across the world also watch the club football on a regular basis. And with the ever-growing competitive nature of the sport, the field of soccer analytics is becoming more important with every passing day. The most interesting aspect about soccer is that it is a dynamic game and a team's success relies both on team strategies and individual player contributions. Hence, it is generally more difficult to develop sophisticated statistical methods for soccer as compared to, say baseball, where events are more discrete. 

% French mathematician \citet{laplace12} hypothesized that if we know the location and momentum of an atom precisely, we can determine its past and future values. 
In soccer, the soccer-ball is certainly like the principal atom in the game, around which the players on the pitch trace their path. In soccer, a shot gives momentum to the atom of the game, read soccer ball, to reach its intended destination, the back of the net. In this aspect, shooting ability of a player has always been a key factor in building team strategies, as the outcome of a game is decided based on the number of goals scored by the two teams. So, it is immensely important to identify players who have great positioning sense and can successfully convert the chances to score more frequently than others. And there comes the need to define a measure which can capture the information about a player's positioning sense and how good he is in converting shots to goals. 

A simple way of looking at this problem is to calculate the proportion of shots one converts successfully, which works around the assumption that probability of scoring a goal is the same for all shots taken by a certain player. While it gives some idea about the efficiency of a striker, it has its shortcomings and would fail to serve the purpose from many aspects. For example, this crude measure only gives the percentage of shots converted successfully, but it does not tell anything about the players' positioning senses. Besides, this measure fails to identify difficulty of chances. For example, a player scoring 2 goal from 4 shots, all from several yards out of the box is considered equally efficient as another who scores 1 from 2 shots from inside the box, and that is definitely not the case. Naturally, one should take into account many other factors - location, time, type of play, how the shot was taken and opponents being some examples - to investigate the conversion rates of soccer players.

While the problem of predicting goals has not been overwhelmingly explored in soccer, there are some studies that have tried to predict the outcome of an average shot from an average player from a particular position on the pitch in a particular situation. Of course, there are limitations as we do not know the shot speed, the position of defenders and the exact goalkeeper position but it is worthy to check how much variance we can explain from the existing covariates which, at the first look, seem to be really important. Current researchers in this field have taken various routes to assign a goal-scoring probability to a particular shot and finally, worked on to develop an Expected Goals (ExpG) metric for a team or a player.

In general, there are different types of discrete regression models that people have used so far. For example, \citet{goddard05} used bivariate Poisson and ordered probit regression models to forecast number of goals scored and conceded in soccer matches based on different covariates. The author here analyzed data from English football and learned that the number of goals in a match depends mostly on the recent home performances of the home side and the recent away performances of the away side. In another match-level study, \citet{williams11} studied the effects of altitude on the result of soccer matches in South America. Taking a regression-based approach for change in altitude of the away side, they showed that the win probability decreases significantly for an away team traveling up, while traveling down does not have a significant impact on the results.

Player-related or shot-related approaches have been studied by a few authors as well. One of the most popular methods in this regard is a study in \citet{asa16}. Analyzing data from the Major League Soccer (MLS), the authors showed that the shot conversion rates significantly depend on many factors, including the distance, angle and the body part used when a shot is taken. Simple logistic regression techniques were the main components of this study. In an earlier paper, \citet{moura07} dealt with data from Brazilian first division championship matches and analyzed shots to goal strategies based on the field position from where the ball was collected and the length of the passing sequence leading to a goal. Contrary to the popular belief, they found that the ball possessions resulting in shots to goal usually start in the defensive field of the scoring teams. So far as the passing sequences are concerned, this study showed that either short (2 to 5) or very long (10 or more) sequences significantly improve the conversion rates. This second finding corroborated what \citet{hughes05} found before. 

Finally, we cite two works that used the situation in a small time-window just preceding the event of interest. First, \citet{jordet13} used English Premier League data to investigate how visual exploratory behaviors before receiving the ball affect the probability of completing the next pass. They found that the players who show extensive visual exploratory behaviors right before receiving a ball are more successful in completing a pass. In another recent study, \citet{lucey14} analyzed the spatiotemporal patterns of the ten-second window of play before a shot was fired. They extracted some strategic features from the data and used them to present a method that estimates the likelihood of a team scoring a goal. This study showed that not only is the game phase (corner, free-kick, open-play, counter attack etc.) important, but the features such as defender proximity, interaction of surrounding players and speed of play coupled with the shot location also play an impact on determining the likelihood of a team scoring a goal. A few more relevant discussions can be found in \citet{leitner10}, \citet{clark13}, \citet{lasek13} and \citet{brooks16}.

In this study, we want to look at the problem of modeling shot conversion rates from a slightly different point of view. While each of the above models has its advantages, there are some limitations too. In order to perform the last two methods, one would need the player tracking data from some reputable source and that requires a considerable amount of money. On the other hand, usual binary regression based approaches often show a lack of fit and fail to have good predictive abilities. This has been discussed more in \Cref{sec:motivation}.

Our primary goal here is to build a model that fits the data well, identifies the factors significantly affecting the shot conversion rate and provides an estimate of the probability of scoring a goal if the location and other details are given. In this regard, we use a spatially correlated error process in the binary regression model. Another key contribution of this paper is to quantify the impact of an individual in soccer. In particular, we define two new measures that evaluate the positioning senses and the goal-scoring abilities of the players. As a practical application, we would implement our model, in conjunction with the measures we devise, to Major League Soccer (MLS) data from 2016/17 season.

The paper is organized as follows. A brief description of the data and some exploratory analysis are provided in the next section. In \Cref{sec:motivation}, we provide motivations behind our idea of including a spatial correlation term in a probit regression model. The proposed method and the measures to evaluate shooters' efficiency are presented in \Cref{sec:methods}. Detailed analysis of the data using the proposed model is carried out in \Cref{sec:results} while the final section includes a discussion, along with some important concluding remarks.


\section{Data}
\label{sec:data}

Throughout this study, we use the data for the first half of the MLS 2016/17 season. Since the focus is on the shots to goal conversions, relevant information for all shots (excluding penalties, own goals and cases with some missing data) and their outcomes - encoded as {\it goal} or {\it miss} - are recorded. But, we choose to exclude cases when a shot was taken from beyond the half line, for often these shots are simply the results of some unusual mess-up of the opponent side or the result of a moment of some striker's presence of mind. In other words, a player will not take a shot from beyond the center line unless he is handed a golden opportunity and naturally, such a shot has a higher probability of being converted. So, including these instances are likely to affect our conclusions. After removing these events, the data comprised of 3957 attempts at goal, among which 482 were converted successfully (approximately 12.18\% of the total number of shots). A heat map is shown in \Cref{fig:heatmap} to show the proportions of goal scored from different locations in the field.

\begin{figure}[!thb]
\begin{center}
\includegraphics[width=11cm,keepaspectratio]{R-scatter-plot.eps}
\end{center}
\caption{Heat map of the full data: Total number of attempts were 3957, among which 482 were goals and 3475 were misses. Each pixel represents the proportion of goals scored from that position, color coded on a scale from 0\% (white) to 100\% (black).}
\label{fig:heatmap}
\end{figure}

Then, in order to fit an appropriate model, it is crucial to identify possible factors affecting the conversion of a shot. For that, the first thing we considered are the distance and angle of the location from where a shot is taken. The way we evaluated these two covariates is displayed on the left panel of \Cref{fig:shot-location}. It should be noted that this definition of the angle and the distance of a shot helps us consider them as uncorrelated covariates in this study. The right panel of the same figure shows another important covariate in our study, the keeper's reach. It is evaluated as the shortest distance the goalkeeper has to cover if he stands at the best position corresponding to the location of the shot.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=11cm,keepaspectratio]{R-distance-angle.eps}
\end{center}
\caption{For location X, the distance is evaluated as the Euclidean distance between X and the center (O) of the goalline, while the angle (with appropriate sign) is calculated as the angle XO makes with OY, the line that bisects the field horizontally.}
\label{fig:shot-location}
\end{figure}

Further, exploratory analysis suggested that a transformation of variable is needed for the distance and the angle of a shot. Similar to the study by \citet{asa16}, we found that it is more appropriate to model the conversion rates as linear function of the logarithm of distance and cosine of angles instead of the usual values. The plots in \Cref{fig:rate-distance-angle} support this.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=11cm,keepaspectratio]{R-rate-distance-angle.eps}
\end{center}
\caption{Conversion rates corresponding to different distances and angles.}
\label{fig:rate-distance-angle}
\end{figure}
%elapsed time (scaled down to an interval between 0 and 1) when the shot is taken, 
In addition to the distance, angle and keeper's reach, there are information on several other covariates and five of them are considered in this study - if it is a home game, the half of play when the shot is taken, body part used in the shot (header, left foot, right foot or other), goal difference at the time of the shot i.e. whether the shooter's team is leading, level or trailing, if the shot takes place during the stoppage time and the proportion of shots converted against the same opponent so far. 

Note that the set of covariates we are using includes several situational covariates (home or away, half of the game, goal difference, and the stoppage time) while the last covariate takes into account the strength of the opposing teams and helps us to identify if it significantly affects the rate of conversion. 

A quick summary of the total number and the proportions of shots converted for different levels of the aforementioned five factored covariates in the study are presented in \Cref{tab:covariates}. For example, the first row says that 2232 (56.4\% of the total of 3957) shots were taken by the home side and 286 of them (which amounts to 12.2\% of the shots taken by the home side) were converted successfully. 

\begin{table}[!hbt]
\caption{Conversion rates for different levels of the factored covariates: For goal or miss, proportions given in brackets are conversion rates for that level while the proportion in the last column is the proportion of shots taken for that level}
\centering
\label{tab:covariates}
\vspace{0.1in}
\begin{tabular}{|c|cc|c|}
\hline
Factor & Goal  & Miss & Total  \\
\hline
Home & 286 (12.2\%) &  1946 (87.8\%) & 2232 (56.4\%) \\
Away & 196 (11.4\%) & 1529 (88.6\%) & 1725 (43.6\%) \\
\hline
First half & 209 (12.0\%) &  1536 (88.0\%) & 1745 (44.1\%) \\
Second half & 273 (12.3\%) & 1939 (87.7\%) & 2212 (55.9\%) \\
\hline
Header & 106 (11.0\%) & 855 (89.0\%) & 961 (24.3\%) \\
Left Foot & 130 (12.5\%) & 911 (87.5\%) & 1041 (26.3\%) \\
Right Foot & 242 (12.5\%) & 1698 (87.5\%) & 1940 (49.0\%) \\
Other & 4 (26.7\%) & 11 (73.3\%) & 15 (0.4\%) \\
\hline
Leading & 127 (13.9\%) &  784 (86.1\%) & 911 (23.0\%) \\
Scores level & 218 (11.4\%) & 1687 (88.6\%) & 1905 (48.1\%) \\
Trailing & 137 (12.0\%) & 1004 (88.0\%) & 1141 (28.8\%) \\
\hline
Regulation time & 447 (11.9\%) & 3290 (88.1\%) & 3737 (94.4\%) \\
Stoppage time & 35 (15.9\%) & 185 (84.1\%) & 220 (5.6\%) \\
\hline
\end{tabular}
\end{table}

One can clearly see that the home factor shows only a little change from one level to another. Same is true for first half versus second half. Very few shots were taken by other means than left or right foot or headers and it should not show significant effect in the analysis. For the three different types of body parts, the header effect is slightly less than the other two. As expected, the conversion rate of both the right-footed and left-footed shots are almost same, at least empirically, which accounts for the symmetry of the play. Finally, the conversion rate seems to be slightly higher when the shooter's team is leading while for the stoppage time shots, the conversion rate is 4\% more than the same for regulation time scenarios.

Further taking a deeper look into the covariate corresponding to body-parts, we noticed that the headers were usually taken from a close range, mostly from inside the box. In fact, the median distance from where the headers were taken was 11.2, starkly different from the overall median of 17.6. \Cref{tab:header-other} shows a summary of the distances and angles (in radian) for headers and other shots. The angles are similar for both cases, but the headers are always taken from a closer range. Naturally, the dependence on the covariates might be totally different for headers and other shots, and so, we suggest separate analysis for the two cases. It is worth mention that this has been overlooked in earlier studies. 

\begin{table}[!hbt]
\caption{Comparison of shot locations (mean, median and standard deviation) for headers and other shots}
\centering
\label{tab:header-other}
\vspace{0.1in}
\begin{tabular}{|l|ccc|ccc|}
\hline
& \multicolumn{3}{|c|}{Summary of distance} & \multicolumn{3}{|c|}{Summary of angle} \\
\hline
        &   Mean  & Median & St. Dev. & Median & Mean & St. Dev. \\
\hline
Header      & 11.3 & 11.2 & 3.92 & $-0.007$ & 0.005 & 0.53 \\ 
\hline
Other shots & 22.0 & 20.8 & 8.77 & $0.004$ & 0.019 & 0.62 \\
\hline
\end{tabular}
\end{table}



\section{Motivation of spatially correlated error process}
\label{sec:motivation}

In most of the works discussed before, the usual practice is to fit a logistic regression with appropriate covariates. However, for these models, what we have commonly observed is a lack of goodness-of-fit checks, a moderately high average error rate and the use of wrong resorts to validate a model (e.g. usual R-squared in case of logistic regression). Even for the study by \citet{lucey14} - one of the most sophisticated techniques in our opinion - where conditional random field was used to take into account defender position and attacking context to model the probability of scoring a goal, we found the error rate to be 14.3\%, only 2\% less than what would have been for the usual logistic regression model without even considering defender attributes. And thus, there is a dire need to investigate these models in detail, so as to identify why they fail in practice. In this regard, the reader is instructed to read the \emph{Deadspin} article by \citet{bertin15}, who discussed the inefficiencies of common statistical models in practice in this regard.

We start our analysis with typical binary regression model with probit links for the MLS data we have. For compactness of the paper, we will discuss only the summary of what we observed with the typical model. First of all, for other shots (not headers), only the distance and angle, transformed to logarithmic and cosine scales, respectively, appeared to be significant. It was found that a smaller distance and a bigger angle significantly (both $p$-values very close to 0) improve the conversion rates. The estimate for the {\it keepreach} variable was estimated to be negative with a $p$-value of 0.067 while the same was positive with a $p$-value of 0.097 for the parameter corresponding to the factor when the shooter's team is leading. This shows that the conversion rate somewhat significantly increases when {\it keepreach} is less and when the shooter's team has a positive goal difference. The other covariates in the study showed no significant effect at all. 

On the other hand, for headers, it was found that a smaller distance, smaller keepreach and a higher angle increase the conversion chance significantly (all $p$-values less than 0.05), while the other covariates are not significant.

However, standard model diagnostics techniques revealed that the above models suffer from high deviance and high classification errors, when exposed to cross-validation. And so, we felt that there is a need to improve upon these models to build a feasible one. In this regard, our initial guess was that two shots may not at all be spatially independent and hence, we cannot take independence assumptions for model fitting. 

Recall that we can consider the occurrence of goals to be a spatial point process with the field being the location set. Plotting Ripley's K-function (\citet{dixon02}) is a standard technique to check whether there exists any spatial clustering or not. It is defined as $K(a)=\lambda^{-1}E_a$ where $\lambda$ is the density (number per unit area) of events (here goals), and $E_a$ is the expected number of extra events within distance $a$ of a randomly chosen event. In \Cref{fig:ripley}, $K_{theo}(a)$, $\hat K_{high}(a)$ and $\hat K_{low}(a)$ give us the expected values and confidence intervals of $K(a)$ under the assumption of spatial homogeneity and independence. We can clearly observe that the $K(a)$ values diverge from the expected values indicating possible spatial correlation and more spatial clustering as $a$ increases.

\begin{figure}[!thb]
%% Use the graphics package to insert figures
%% \includegraphics[width=0.8\textwidth]{figure.eps}
% Use \centering to center the table
\begin{center}
\includegraphics[width=0.8\textwidth]{Ripley_HNH.eps}
\end{center}
\caption{Ripley's K-function $K(a)$ is plotted against $a$ for both headed and non-headed goals considering goals as spatial point processes. Under spatial homogeneity and independence, $K(a)=\pi a^2$. Expected values under this hypothesis and confidence intervals are constructed for reference.}
\label{fig:ripley}
\end{figure}

Now that we have a hint that there might be a spatial dependency involved, we test whether there exists spatial autocorrelation or not. We have a binary variable with $1$ denoting a goal and $0$ otherwise. We create a graph of the shot locations based on $k$-nearest neighbor criteria using Euclidean distance (taking $k=63$ as it has been accepted as a thumb rule to take $k=\sqrt{n}$ where $n$ is the number of observations, here $n=3957, \sqrt{n} \approx 63 $). We join edges between two points if a point belongs to the set of $63$ nearest neighbours of another point and vice-versa. Thus, joining the neighbours, we can get the join-count statistics (\citet{cliff81}) to test for spatial autocorrelation. We perform join-count test against the alternative that the number of like joins are more than expected from random under 5\% significance level and found clear evidence for positive spatial auto-correlation as the number of $0-1$ join-counts were significantly lower than the expected values and both $1-1$ and $0-0$ counts were significantly higher than the expected ($p$-value nearly 0 for all the cases).

Note that positive autocorrelation indicates clustering of similar values and hence we conclude that there is sufficient spatial correlation to account for with the correlation between probabilities of shot conversions from two locations on the field decreasing  with increasing distance between the locations. Thus, one cannot perform the analysis considering spatial independence of the shots. Driven by the exploratory findings above, our assumption is that the shots are dependent even if they are from different matches and different players. This is like seeing it as a big shooting experiment in soccer which has been done on the same hypothetical field under different circumstances. We further emphasize that, to the best of our knowledge, this interesting and very important phenomena of spatial dependence has never been accounted for before in soccer. To tackle this, we introduce a spatially correlated error process in the standard probit model with the spatial covariance following an exponential decay law and use complete Bayesian framework to estimate the parameters and variances correctly. This idea, in form of spatial autoregressive models, is close to what has been used in spatial econometrics before (\citet{elhorst2017transitions} and \citet{baltagi2016spatial}). 


A detailed description of the model is provided in the next section.





\section{Methods}
\label{sec:methods}



\subsection{Proposed model}
\label{subsec:proposed-model}

We will start with a binary probit regression model equipped with a spatially correlated error process, in light of the discussion in the previous section. For shots $i=1,\hdots,N$, suppose $Y_i$ is the binary outcome (which takes the value 1 if the $i$th 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 \sim \ber(p_i)$, $p_i=P(r_i>0)$, and that $r_i$ depends on the values of the covariates following a simple linear regression model. We would 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 be less for locations far apart. These assumptions support our natural findings while exploring the data in the previous section. Now, let us use $Y$ to denote the vector of all outcomes, $r$ to denote the vector of $r_i$'s, $p$ to denote the vector of $p_i$'s and let $X$ be the design matrix such that $X\dash =[X_1:\hdots:X_N]$. 

We are going to consider a hierarchical structure. First of all, $Y_i=\mathbb{I}(r_i>0)$ and
\begin{equation}
\label{eqn:hierarchical-model}
r_i = X_i\dash \theta + \eps_i,
\end{equation}
where $\theta$ is a parameter vector of appropriate order, corresponding to the covariates and $\eps_i$ is the error process. Note that the mean function $\mu_i = X_i\dash \theta$ is considered to be an additive combination of the effects of some continuous and some discrete covariates. We would further assume that the error process has the following structure:
\begin{equation}
\eps_i = w_i + e_i.
\end{equation}

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,\sigma^2)$ 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 
\begin{equation}
\label{eqn:spatial-covariance}
\cov(w_i,w_j) = \sigma_w^2 \ \exp(-\phi\norm{s_i-s_j}), 
\end{equation}
where the distance function $\norm{s_i-s_j}$ is taken as the Euclidean distance of the two locations. Writing $w=(w_1,\hdots,w_N)$, we will denote the covariance matrix of $w$ by $\sigma_w^2\Sigma_w$. Thus, we can now write the full model in the following form:
\begin{eqnarray}
\label{eqn:full-model}
Y_i &=& \mathbb{I}(r_i>0), \\ 
%p_i = 1-\Phi(Z_i) \nonumber \\
r &=& X \theta + w + e, \nonumber \\
\text{such that } \ w &\sim& N(0,\sigma_w^2\Sigma_w), \nonumber \\
\text{and } \ e &\sim& N(0,\sigma^2I). \nonumber
\end{eqnarray}

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 $\theta$, we will consider improper Jeffrey's prior distribution. On the other hand, we will assume the two variance components $\sigma^2,\sigma_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 prior we associate with $\sigma^2$ is the inverse gamma (IG) distribution with parameters $a,b$. We would take $a>1$ to avoid improper posterior distributions.

And finally, the parameter $\phi$ (refer to \Cref{eqn:spatial-covariance}) will be kept fixed throughout the Bayesian analysis. Ideally, one should consider appropriate prior distribution and estimate it within the Bayesian model as well. But, simulation studies show that the computational burden reduces to a great extent and the performance does not deteriorate if we keep it fixed and that motivated us to make this adjustment. In order to find out the best possible value for $\phi$, we would consider a cross-validation scheme. The possible choices for $\phi$ used in the study were between 0.05 and 1. Based on the relationship $e^{-\phi d}\approx 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 $30\%$ of all the shots after estimating the parameters from the other $70\%$ and would find out the mean squared error for those predictions. For shots considered for the validation purpose (let us denote them by $i \in V=\{v_1,\hdots,v_m\}$), we would predict the probability of them being converted to goals and let the predicted probabilities be denoted by $\hat p_i$ while $Y_i$'s are the actual outcomes. Then, for each $\phi$, we calculate the mean squared error of the predictions using the following formula:
\begin{equation}
\label{eqn:validation-mse}
\text{MSE}_\phi = \frac{1}{m}\sum_{i\in V} (Y_i - \hat p_i)^2.
\end{equation}

The value of $\phi$ with the least mean squared validation error will be our optimal choice for the analysis. The prediction procedure is discussed in detail in \Cref{subsec:prediction}.


\subsection{Posterior distribution and Gibbs sampling}
\label{subsec:gibbs}

Throughout the discussion below, $K$ will always indicate a constant term, that may vary from time to time. Using the full model described in \Cref{eqn:full-model} and the prior distributions described in the paragraph following that equation, we can write the joint posterior distribution as:
\begin{eqnarray}
\log \pi(r,\theta,\sigma^2,w\given Y) &=& K + \sum_{i=1}^N\{Y_i\log(P(r_i>0))+(1-Y_i)\log(P(r_i\leq0))\} - \frac{b}{\sigma^2} \nonumber \\
\label{eqn:joint-posterior}
&&   - \frac{\norm{r-X\theta-w}^2}{2\sigma^2 } - \frac{w\dash \Sigma_{w}^{-1}w}{2\sigma^2}  -(a+N+1)\log \sigma^2. %- \frac{1}{2v^2}\norm{\theta}^2
\end{eqnarray}

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 $\sigma^2$, we can get the full conditional distribution as:
\begin{equation*}
\log \pi(\sigma^2\given Y,r,\theta,w) = K - (a+N+1)\log \sigma^2 - \frac{1}{\sigma^2}\biggl[b + \frac12\norm{r-X\theta-w}^2 + \frac12w\dash \Sigma_{w}^{-1}w\biggr] 
\end{equation*}
and thus, 
\begin{equation}
\label{eqn:sigma-posterior}
(\sigma^2\given Y,r,\theta,w) \sim IG\biggl(a+N,b+\frac12\norm{r-X\theta-w}^2 + \frac12w\dash \Sigma_{w}^{-1}w\biggr).
\end{equation}

For $\theta$, we get that
\begin{equation*}
\log \pi(\theta \given Y,r,\sigma^2,w) = K - \frac12 \theta\dash \biggl(\frac{X\dash X}{\sigma^2}\biggr)\theta + \frac12 \cdot 2\theta\dash  \cdot \frac{X\dash (r-w)}{\sigma^2}.
\end{equation*}
Straightforward calculations then tell us that the conditional posterior of $\theta$ is 
\begin{equation}
\label{eqn:theta-posterior}
(\theta \given Y,r,\sigma^2,w) \sim N\bigl((X\dash X)^{-1}X\dash (r-w),\sigma^2(X\dash X)^{-1}\bigr).
\end{equation}

On the other hand, using similar techniques as before, we can get the posterior distribution for $w$ as:
\begin{equation}
\label{eqn:w-posterior}
(w \given Y,r,\theta,\sigma^2) \sim N\bigl((I+\Sigma_{w}^{-1})^{-1}(r-X\theta),\sigma^2(I+\Sigma_{w}^{-1})^{-1}\bigr).\end{equation}

Note that  the means of the conditional distributions of both $\theta$ 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 $\sigma$ and so, one needs to compute the inverses of three big matrices, namely $X\dash X$, $\Sigma_w$ and $(I+\Sigma_w^{-1})$, only once throughout the MCMC analysis. Hence, directly using the above conditional distribution of $w$ 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 distribution in \Cref{eqn:joint-posterior} clearly shows 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.
\begin{equation}
\label{eqn:r-posterior}
(r_i\given Y,\theta,w,\sigma^2) \sim \begin{cases} TN(\alpha_i,\sigma^2;0,\infty) & \text{if } \ Y_i=1 \\ TN(\alpha_i,\sigma^2;-\infty,0) & \text{if } \ Y_i=0 \end{cases}
\end{equation}
In the above, $\alpha=X\theta+w$ and $TN(\mu,\tau^2;a,b)$ denotes a normal distribution with mean parameter $\mu$ and variance parameter $\tau^2$, truncated in the interval $[a,b]$.


\subsection{Future prediction}
\label{subsec:prediction}

In order to make a new prediction for a shot from location $s\dash$, we will take resort to the posterior estimates obtained from the Gibbs sampler. Let us denote the outcome of the shot from $s\dash$ by $Y(s\dash)$ ($X(s\dash)$, $r(s\dash)$ and $w(s\dash)$ are defined accordingly). Note that the posterior predictive distribution $f(Y(s\dash)\given 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 would draw samples for $r,w,\theta,\sigma^2$ using the conditional posterior distributions described in \Cref{subsec:gibbs}. After that, we would draw samples for $w(s\dash)$ using the conditional distribution of $(w(s\dash)\given w,\sigma^2)$ as discussed below. Using the above sampled observations, an estimate for $r(s\dash)$ can be obtained by taking a sample from a normal distribution with mean $X(s\dash)\dash \theta+w(s\dash)$ and variance $\sigma^2$. Finally, we set $Y(s\dash)=\ind(r(s\dash)>0)$. Alternately, one can also use $Y(s\dash) \sim \ber(p(s\dash))$, where $p(s\dash)=P(r(s\dash)>0)$ can be computed using the sampled values of $\theta$, $\sigma^2$ and $w(s\dash)$.

To compute the conditional distribution of  $(w(s\dash)\given w,\sigma^2)$, we start with the fact that
$$\left(\begin{array}{c} w(s\dash) \\ w\end{array}\right) \sim N\biggl(0,\sigma^2\left[\begin{array}{cc} 1 & \Sigma_{s-s\dash}\dash  \\ \Sigma_{s-s\dash} & \Sigma_{w} \end{array}\right]\biggr).$$

Here, $\Sigma_{s-s\dash}$ is the column vector denoting the covariance of $w(s\dash)$ with the elements of $w$. Following \Cref{eqn:spatial-covariance}, the $j$th element of $\Sigma_{s-s\dash}$ is $\exp(-\phi\norm{s_j-s\dash})$. Using the principle of conditional distribution for multivariate normal distribution, we can then say that
\begin{eqnarray}
\label{eqn:wnew-conditional}
w(s\dash)\given w &\sim& N(\mu_{\text{new}},\sigma^2_{\text{new}}), \\ 
\text{where } \ \sigma_{\text{new}}^2 &=& \sigma^2\bigl(1-(\Sigma\dash _{s-s\dash}\Sigma_w^{-1}\Sigma_{s-s\dash})\bigr), \nonumber \\ 
\mu_{\text{new}} &=& \Sigma\dash _{s-s\dash}\Sigma_w^{-1}w. \nonumber
\end{eqnarray}


\subsection{Measures to quantify player abilities}
\label{subsec:measures}

Here, we propose two measures for rating a player's scoring ability on the pitch based on our model. To proceed, at first, we remove the shots taken by particular player of interest from our data. This will be our training data to build the model and depending on the output from this model, we will quantify each player in two ways, as described below.

To start with, we would like to introduce, what we call, the {\it Shooting Prowess (SP)} of a player. Suppose, in a total of $g_k$ matches, a player $k$ has taken $m_k$ shots (calculated over a time period, possibly over a whole season or for the entire career) in total. The end result of the $i$th shot, $i=1,\hdots,m_k$, is denoted by $Y_{i}$ ($0$ for miss and $1$ for goal) and the expected probability of the shot predicted from the training model is denoted by $\hat p_i$. Then, ${SV}_{ki}$ (Shooting Value) of that shot is denoted by $(Y_i-\hat p_i)$. Note that this quantity lies in the interval $[-1,1]$; 1 denoting the best possible conversion and $-1$ denoting the worst miss. Using these shooting values for every shot, we define the measure ${SP}_k$ as follows:
\begin{equation}
\label{eqn:sp}
{SP}_k= \frac{1}{g_k}\sum_{i=1}^{m_k}{{SV}_{ki}} = \frac{1}{g_k}\sum_{i=1}^{m_k} {(Y_i-\hat{p_i})}
\end{equation}

A discussion of why we put forward such a measure to quantify a player's ability to score goals is warranted at this point. For any model to analyze such binary data, there are still some variances unaccounted for. It is a logical assumption that for a game like soccer, the conversion of a shot to goal will always depend on the player's abilities. So, we suggest that the difference between the end result and the modeled probability of a shot being converted to goal captures the shooter's interference in that shot. On one hand, a high negative value will indicate a player missing a golden opportunity, thereby putting a negative impact on his rating while on the other, successfully converting a shot with less theoretical probability would capture the player's more than average ability to score goals. 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. For example, a player taking $2$ shots in a match, with a SV of $0.5$ each, will have lower ratings than another player who takes $10$ shots in the same match, with a SV of $0.5$ each, thereby proving that the second player is better and more consistent than the first player. A more sophisticated way to quantify the SP is to consider a weighted average of the sum of shooting values per game, the weights being proportional to the number of minutes played in every match. Since we did not have enough data on this, we chose to define it in the above fashion.

Next, we would like to quantify a player's {\it 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. Hence, summing up the predicted probability of shot conversion for a particular player can give us a perfect idea about his positioning sense. Necessarily, the ${PS}_k$ (Positioning Sense) of the $k$th player is defined by
\begin{equation}
\label{eqn:ps}
{PS}_k=\frac{1}{g_k}\sum_{i=1}^{m_k} {\hat p_i}.
\end{equation}

Once again, this measure captures the average effect per game and hence, has the aforementioned nice property similar to the other measure.




\section{Results}
\label{sec:results}


\subsection{Fitting the model}

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

For other shots, let $Y_i=\ind(r_i>0)$ denote whether $i$th shot was a goal. Then, we set $r_i=\mu_i+w_i+\epsilon_i$, where $w_i,\epsilon_i$ are as described in the previous section and $\mu_i$ is the mean function defined by
\begin{eqnarray}
\mu_i &=& \beta_0 + \beta_1\log(D_i) + \beta_2 \cos(A_i) + \beta_3 K_i + \beta_4 P_i + \beta_5\ind\{H_i=1\} + \beta_6\ind\{F_i=1\}  \nonumber \\
\label{eqn:model-mean}
&&  + \sum_{j=0}^1\beta_{7,j}\ind\{GD_i=j\} + \beta_8\ind\{S_i=1\} + \sum_{k=2}^3 \beta_{9,j}\ind\{BP_i=j\}.
\end{eqnarray}

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 $\beta_{7,-1}=\beta_{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 (\Cref{eqn:model-mean}), we fit the model described in \Cref{sec:methods}. First of all, using the procedure described towards the end of \Cref{subsec:proposed-model}, we evaluate the best choice of $\phi$ to be used in the spatial covariance function (\Cref{eqn:spatial-covariance}). It is obtained to be 0.1 (for both cases) and following the approximate relation $e^{-\phi d}\approx 0.05$, we can say that the spatial correlation is insignificant for locations approximately 30 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 500 to evaluate the estimates, standard errors and the confidence intervals of the parameters. Summary of the results for the `headers' are shown in \Cref{tab:estimates-headers}, while the same for `other shots' are displayed in \Cref{tab:estimates-others}.

\begin{table}[!hbt]
\caption{Results from fitting our model to headers (total number of events was 961)}
\centering
\label{tab:estimates-headers}
\vspace{0.1in}
\begin{tabular}{|lccc|}
\hline
Parameter & Estimate & Standard error & Confidence interval \\
\hline
$\beta_0$ (Intercept)  		    & $-0.544$ 	& 1.341 & $(-3.212,2.098)$ \\
$\beta_1$ (log(Distance)) 		& $-0.642$ 	& 0.486 & $(-1.588,0.314)$ \\
$\beta_2$ (cos(Angle)) 			& $0.454$ 	& 1.070 & $(-1.617,2.549)$ \\
$\beta_3$ (Keeper's reach)  	& $-0.168$ 	& 0.052 & $(-0.276,-0.068)$ \\
$\beta_4$ (Opponent)     		& $-0.688$ 	& 1.379 & $(-3.333,1.959)$ \\
$\beta_5$ (Home)		 	    & $0.024$ 	& 0.142 & $(-0.266,0.305)$ \\
$\beta_6$ (First half) 		    & $-0.031$ 	& 0.153 & $(-0.317,0.264)$ \\
$\beta_{7,0}$ (score level) 	& $0.056$ 	& 0.152 & $(-0.226,0.346)$ \\
$\beta_{7,1}$ (leading) 	    & $0.232$ 	& 0.194 & $(-0.155,0.610)$ \\
$\beta_8$ (Stoppage time)	    & $0.439$ 	& 0.258 & $(-0.082,0.918)$ \\
$\sigma^2$ (Variance) 	    	& $0.996$ 	& 0.069 & $(0.868,1.134)$ \\
\hline
\end{tabular}
\end{table}

\begin{table}[!hbt]
\caption{Results from fitting our model to other shots (total number of events was 2996)}
\centering
\label{tab:estimates-others}
\vspace{0.1in}
\begin{tabular}{|lccc|}
\hline
Parameter & Estimate & Standard error & Confidence interval \\
\hline
$\beta_0$ (Intercept)  		    & $1.401$ 	& 1.061 & $(-0.580,3.512)$ \\
$\beta_1$ (log(Distance)) 		& $-1.513$ 	& 0.296 & $(-2.146,-0.948)$ \\
$\beta_2$ (cos(Angle)) 			& $1.163$ 	& 0.717 & $(-0.170,2.483)$ \\
$\beta_3$ (Keeper's reach)  	& $-0.052$ 	& 0.031 & $(-0.110, 0.007)$ \\
$\beta_4$ (Opponent)     		& $-0.544$ 	& 0.783 & $(-2.065,0.974)$ \\
$\beta_5$ (Home)		 	    & $0.001$ 	& 0.077 & $(-0.154,0.153)$ \\
$\beta_6$ (First half) 		    & $0.026$ 	& 0.079 & $(-0.134,0.176)$ \\
$\beta_{7,0}$ (score level) 	& $0.001$ 	& 0.095 & $(-0.172,0.182)$ \\
$\beta_{7,1}$ (leading) 	    & $0.192$ 	& 0.091 & $(0.014,0.370)$ \\
$\beta_8$ (Stoppage time)	    & $-0.068$ 	& 0.170 & $(-0.410,0.258)$ \\
$\beta_{9,2}$ (Left foot) 		& $0.451$ 	& 0.423 & $(-0.310,1.366)$ \\
$\beta_{9,3}$ (Right foot) 		& $0.474$ 	& 0.422 & $(-0.234,1.330)$ \\
$\sigma^2$ (Variance) 	    	& $0.998$ 	& 0.040 & $(0.925,1.077)$ \\
\hline
\end{tabular}
\end{table}

From the results, we can see that the effect of logarithm of distance is significantly negative for non-headers, which 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. While 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, we have found that the angle of the shot is not a significant variable in our model. This happens probably because the effect of the angle is already being explained by the spatial process. 

For headers, however, neither the distance nor the angle is significant. In fact, no other covariate except keeper's reach is significant for headers. The chance of converting from a header increases when the keeper's reach is less. But, barring that, it is clear that the inclusion of the spatially correlated process can explain most of the variability in the data.

%On the other hand, it is worth mention that if we use the values of the angles instead of the cosine transformation, the usual logistic regression model results in an unrealistic conclusion which says that a bigger angle increases the goal-scoring chance. But, even in that case, conclusions from fitting our model remain unchanged. 

Further, for non-headed shots, so far as the other covariates are concerned, we can see that the chances increase significantly when the shooter's team leads. It is in line with the notion that a positive goal difference can boost the morale of the team to a great extent. However, we observed that the opponent effect, the home advantage, the half or the stoppage time does not play a significant role in explaining the shot conversion rates. The conclusion is same for left or right footed shots, though it is worth noting that the coefficients corresponding to these two covariates are very close to each other (0.451 and 0.474), thereby supporting the nature of symmetry in the game of soccer. 


\subsection{Diagnostics}

We next delved into the diagnostics of the model, and compared the performance of our model to two different models. First, we considered 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 \citet{klier08}. 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.

At first, we calculated 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 $\sum_{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 computed it to be 0.076, while for the headers it was approximately 0.073. Note that for both cases, this score is approximately 20\% less than what we observed in SLRM or KMM. 

Another standard method to evaluate probabilistic prediction is the log score, which is defined by $\sum_{t=1}^N\{o_t\log f_t+(1-o_t)\log(1-f_t)\}/N$. 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 (and KMM) yields log score of approximately $-0.302$ (and $-0.339$), and that is nearly $34\%$ (and 51\%) less than what our method gives (log score of about $-0.225$). Similar phenomena was observed for non-headed shots as well.

Next, we evaluated the proportion of times the models could classify the outcomes correctly. Here, the improvement with our model was only about 10\%, as compared to SLRM or KMM.

The above results are summarized in table \ref{tab:comparison}.

\begin{table}[!hbt]
\caption{Comparison of three candidate models}
\centering
\label{tab:comparison}
\vspace{0.1in}
\begin{tabular}{|cl|ccc|}
\hline
& & SLRM & KMM & Our model \\
\hline
        & Brier Score  	& 0.091 	 & 0.101        & 0.073 \\
Headers & log score     & $-301.738$ & $-338.876$   & $-224.735$ \\
        & error \% 	    & $11.03\%$  & $11.654\%$   & $10.406\%$ \\
\hline
            & Brier Score  	& 0.094 	 & 0.098        & 0.076 \\
Other shots & log score     & $-958.212$ & $-1043.567$  & $-728.228$ \\
            & error \% 		& $11.983\%$ & $11.983\%$   & $10.748\%$ \\
\hline
\end{tabular}
\end{table}


As a final piece, we wanted to see how accurately we can predict 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 \citet{king01}. Since we did not have that information, we had to rely on usual prediction methods (refer to \Cref{subsec:prediction}). We used cross-validation techniques to evaluate how good the predictions are, by taking a subset of the full data (approximately 90\%) to learn the model and applying it on the rest (approximately 10\%). Afterwards, in order to evaluate the prediction efficiency, we followed the general scoring rules outlined by \citet{merkle13}. 

In the aforementioned study, the authors discussed the beta family of proper scoring rules, proposed by \citet{buja05} for binary regression problems, in detail. Using different combinations of the parameters $(\alpha,\beta)$ for the beta distribution, one can choose the most appropriate scoring rule. For example, $\alpha=\beta=0$ corresponds to log score while $\alpha=\beta=1$ corresponds to Brier score. On the other hand, \citet{merkle13} presented another way of thinking about the parameters, where one can assume that the value of $\alpha/(\alpha+\beta)$ (or $1-\alpha/(\alpha+\beta)$) is the cost of a false positive (or false negative). In light of this, one should note that scoring rules with $\alpha<\beta$ 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 $\alpha>\beta$ emphasize high-probability forecasts. Using the R package `scoring' (\citet{Scoring}), we have computed the scores for different values of $\alpha/(\alpha+\beta)$. The logarithm of the scores for the three models are plotted against different cost functions in \Cref{fig:score-cost}. 

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=0.8\textwidth]{R-score-cost.eps}
\end{center}
\caption{Logarithm of the scores for the two models, corresponding to different cost values.}
\label{fig:score-cost}
\end{figure}

We can see that the performances of SLRM and KMM are very similar. But, our model shows consistently better (higher score function) performance when the cost for a false positive is not less than 0.6. 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 $\alpha=\beta$, when the cost is equal for goal and miss, calculations showed that for our model, the score function is approximately $15\%$ (and 19\%) higher than that of the SLRM (and KMM).

To conclude this section, we emphasize that goals are example of rare events (\citet{anderson13}), 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 in such a rare events regression model 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 existing method (SLRM), and also, our specification of the spatial correlation in the probit model works better than other popular choice (KMM). 


\subsection{SP and PS of the players}
\label{subsec:sp-ps}

We now proceed to compute the two measures described in \Cref{subsec:measures}. 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 \Cref{fig:histogram-sp-ps}. 

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=\textwidth]{R-sp-ps.eps}
\end{center}
\caption{Histogram of the shooting prowess (left) and positioning senses (right) for 381 players, calculated based on fitted probability and actual outcomes.}
\label{fig:histogram-sp-ps}
\end{figure}

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. On the other hand, we can see some clear outliers towards the right side of the histogram. We confirmed that these outliers mostly represent the players who took only few shots (usually one or two) all of which were converted successfully. Similarly for the figures on the right, we can see that most players showed average 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 and to exclude the aforementioned outliers, we considered players with more than 10 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 \Cref{tab:sp-ps}. In that table, we also present the data for the top players according to SPH and PSH. 

\begin{table}[!hbt]
\caption{Real-life statistics of best players as identified by our model}
\centering
\label{tab:sp-ps}
\vspace{0.1in}
\begin{tabular}{|lccccc|}
\hline
\multicolumn{6}{|c|}{Players with best SPS (ability to score with a shot)} \\
\hline
Player 	& SPS & PSS & Shots (on target) & Goals & Assists \\
\hline
Michael Barrios 	& 0.294	& 0.173	& 50 (26) & 9 & 2 \\
Giles Barnes 	    & 0.291	& 0.109	& 51 (19) & 6 & 3 \\
Ignacio Piatti 	    & 0.277	& 0.370	& 95 (42) & 17 & 6 \\
Giovani Dos Santos 	& 0.168 & 0.216	& 63 (31) & 14 & 12 \\
Sebastian Giovinco 	& 0.112 & 0.288	& 177 (59) & 17 & 15 \\
\hline
\multicolumn{6}{|c|}{Players with best PSS (positioning sense for shooting)} \\
\hline
Player & PSS & SPS & Shots (on target) & Goals & Assists \\
\hline
Ola Kamara 		& 0.509	& 0.074	    & 80 (38) & 16 & 2 \\
David Accam     & 0.418 & $-0.132$	& 71 (28) & 9 & 5 \\
Yura Movsisyan 	& 0.411	& $-0.117$	& 68 (28) & 9 & 3 \\
C.J. Sapong 	& 0.406	& $-0.139$ 	& 60 (22) & 7 & 2 \\
Jack McInerney  & 0.403	& $-0.189$	& 37 (18) & 5 & 0 \\
\hline
\multicolumn{6}{|c|}{Players with best SPH (ability to score with a header)} \\
\hline
Player 	& SPH & PSH & Shots (on target) & Goals & Assists \\
\hline
Drew Moor               & 0.291 & 0.084	& 21 (6)   & 3  & 0 \\
Bradley Wright-Phillips & 0.204 & 0.068	& 103 (56) & 24 & 5 \\
Chris Wondolowski       & 0.150 & 0.100	& 78 (31)  & 12 & 3 \\
Dominic Dwyer           & 0.134 & 0.133 & 119 (44) & 16 & 3 \\
Clint Dempsey           & 0.104 & 0.063 & 65 (22)  & 8 & 2 \\
\hline
\multicolumn{6}{|c|}{Players with best PSH (positioning sense for heading)} \\
\hline
Player & PSH & SPH & Shots (on target) & Goals & Assists \\
\hline
Kendall Waston  & 0.298 & $-0.076$	& 21 (10)  & 3  & 1 \\
Quincy Amarikwa & 0.291 & $-0.201$	& 50 (12)  & 3  & 3 \\
Juan Agudelo    & 0.230 & $-0.088$  & 52 (18)  & 7  & 5  \\
Ola Kamara 		& 0.229	& 0.021	    & 80 (38)  & 16 & 2 \\
David Horst     & 0.230 & $-0.088$  & 21 (11)   & 3  & 0  \\
\hline
\end{tabular}
\end{table}
 
We noted that Bradley Wright-Phillips, Sebastian Giovinco, Ignacio Piatti, Dominic Dwyer and Giovani Dos Santos 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. On the other hand, Ola Camara was the sixth highest goal-scorer of the season and we found that he has one of the best positioning senses, both in terms of shooting and heading. 

The results for the heading abilities were something intriguing. At first, we like to point out the special cases of Drew Moor, Kendall Waston 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. On the other hand, Juan Agudelo and Quincy Amarikwa showed really good positioning senses when it comes to heading. But, their conversion abilities are poor and that has affected their overall stats. In fact, they were also not the first choice strikers of their respective teams. Finally, we note that Clint Dempsey, one of the finest players United States of America has ever produced, has brilliant heading prowess, but his SP quotient was not one of the best in the league.

It is worth mention here that the training data we used in our model contained information only from a part of the season and the above measures are always calculated based on other players' 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.



\section{Summary}
\label{sec:conclusion}

In conclusion, we have presented an alternate way of analyzing the conversion rate of shots in soccer. Using a spatially correlated error process, we have shown that our model fits the data well. Our proposed method has exhibited improved Brier score, log score and predictive abilities, as compared to the logistic regression model that is in practice at the moment. On the other hand, our way of specifying the spatial correlation term is also more appropriate than other possible way for probit regression models. Another key contribution of our work is to properly quantify the shooters' efficiency and positioning sense. To the best of our knowledge, this is the first paper to introduce such measures. We have also established that the quantitative measures are good enough to identify the best strikers and thus, the teams would certainly be benefited from using this concept. In a more specialized way, conditioning on different opponents, one can use relevant data to find out the best shooters against them and that would help the coaches in making strategies and tactics against particular teams.

On a related note, we also mention that using a similar measure for the goalkeepers, one can quantify the `saving prowess' to single out the best custodians in the league. We chose not to do that in this study, for we think the goalkeepers' abilities should not be evaluated solely based on the number of goals conceded, but also on the number of blocks and saves made by the goalkeeper. And that is one of the future directions we are interested in. We emphasize that the methods we described here are applicable to binary data (goals and misses), but it can be extended to the multinomial case (blocked, goal, missed, hit by post or saved) as well and that is what needs to be done to identify the best goalkeepers in business. 

A minor issue with our model is that if we have the exact same co-ordinates for two shot locations, then the covariance matrix associated with the error process will be singular and the model will fail. But, theoretically, it is a measure zero event and it happens very rarely in practice. For example, we did not encounter such an issue in our data. To understand it a bit further, observe that penalty kicks are one such kind of events, where the shots are taken from a fixed location inside the box. But, penalty kicks should not be considered in a study like ours because the outcome of the spot-kicks do not really depend on the distance, angle or other physiological attributes. Rather, it only depends on the penalty taker's abilities and usually every team has a specialist for this job. Hence, quite rightly, we left out these cases from the data before running the analysis. Now, apart from penalty kicks, even if we have two or more shots fired from the same exact location, we can easily make adjustments as follows. Suppose, $z_1,..,z_k$ are the shots fired from exact same location. Then, we randomly pick one shot from $z_1,..,z_k$  and leave out the other shots from the data. Next, we fit the model based on the available data and predict the outcomes of the left-out shots from the model. We calculate the sum of squared prediction errors in this case and continue this procedure for all of the $k$ shots. And finally, we pick the shot with the lowest error and include it, discarding the rest, in our final model.

We also note that the parameter $\phi$ (refer to \Cref{eqn:spatial-covariance}), which denotes the extent of spatial dependence, is not estimated directly from our model. Instead, we use a cross-validation approach to find out the best choice of $\phi$ from a grid. Naturally, this does not allow us to test the hypothesis of no spatial dependence, cf. \citet{anselin2001rao}. While it might seem a disadvantage, we want to point out that this is a common practice in many spatial studies, see, for example, \citet{sahu06}. In fact, in this study, we built the model after getting sufficient evidence in favour of spatial autocorrelation and hence, afterwards, we are more interested in explaining the variation among conversion of shots explained by the covariates and the spatial error process and then, how accurately we can predict the outcomes. The results establish that our model works well to serve that purpose.

Finally, we believe that we could have improved the predictions if we had more intense data containing information like the shot speed, the attacking speed and the position of the defenders in the shooting range of the shooter. As our model is really flexible, in future, we aim to incorporate those attributes in the discussion to conduct a more comprehensive study.


\section{Acknowledgement}
\label{sec:acknowledgement}

We thank Dr. Chris Anderson, the author of {\it The Numbers Game : Why Everything You Know About Soccer Is Wrong} and a pioneer in soccer analytics, for providing his valuable feedback about our approach to address the research problem. Our gratefulness also goes to American Soccer Analysis (ASA) for providing the data in this study. 


\bibliographystyle{agsm}
\bibliography{reference}



\end{document}
main.bbl
\begin{thebibliography}{xx}

\harvarditem{{American Soccer Analysis}}{2016}{asa16}
{American Soccer Analysis}  \harvardyearleft 2016\harvardyearright , `Expected
  goals 2016 - xgoals explanation'.
\newline\harvardurl{http://www.americansocceranalysis.com/explanation/}

\harvarditem{Anderson \harvardand\ Sally}{2013}{anderson13}
Anderson, C. \harvardand\ Sally, D.  \harvardyearleft 2013\harvardyearright ,
  {\em The numbers game: Why everything you know about soccer is wrong},
  Penguin.

\harvarditem{Anselin}{2001}{anselin2001rao}
Anselin, L.  \harvardyearleft 2001\harvardyearright , `Rao's score test in
  spatial econometrics', {\em Journal of statistical planning and inference}
  {\bf 97}(1),~113--139.

\harvarditem[Baltagi et~al.]{Baltagi, LeSage \harvardand\
  Pace}{2016}{baltagi2016spatial}
Baltagi, B., LeSage, J. \harvardand\ Pace, R.  \harvardyearleft
  2016\harvardyearright , {\em Spatial Econometrics: Qualitative and Limited
  Dependent Variables}, Advances in econometrics, Emerald Publishing Limited.
\newline\harvardurl{https://books.google.com/books?id=O4K9jwEACAAJ}

\harvarditem{Bertin}{2015}{bertin15}
Bertin, M.  \harvardyearleft 2015\harvardyearright , `Why soccer's most popular
  advanced stat kind of sucks', {\em Deadspin} .
\newline\harvardurl{http://bit.ly/2qtKK3E}

\harvarditem[Brooks et~al.]{Brooks, Kerr \harvardand\ Guttag}{2016}{brooks16}
Brooks, J., Kerr, M. \harvardand\ Guttag, J.  \harvardyearleft
  2016\harvardyearright , `Using machine learning to draw inferences from pass
  location data in soccer', {\em Statistical Analysis and Data Mining: The ASA
  Data Science Journal} {\bf 9}(5),~338--349.

\harvarditem[Buja et~al.]{Buja, Stuetzle \harvardand\ Shen}{2005}{buja05}
Buja, A., Stuetzle, W. \harvardand\ Shen, Y.  \harvardyearleft
  2005\harvardyearright , `Loss functions for binary class probability
  estimation and classification: Structure and applications', {\em Working
  draft, November} .

\harvarditem[Clark et~al.]{Clark, Johnson \harvardand\ Stimpson}{2013}{clark13}
Clark, T., Johnson, A. \harvardand\ Stimpson, A.  \harvardyearleft
  2013\harvardyearright , Going for three: Predicting the likelihood of field
  goal success with logistic regression, {\em in} `The 7th Annual MIT Sloan
  Sports Analytics Conference'.

\harvarditem{Cliff \harvardand\ Ord}{1981}{cliff81}
Cliff, A.~D. \harvardand\ Ord, J.~K.  \harvardyearleft 1981\harvardyearright ,
  {\em Spatial processes: models \& applications}, Taylor \& Francis.

\harvarditem{Dixon}{2002}{dixon02}
Dixon, P.~M.  \harvardyearleft 2002\harvardyearright , `Ripley's k function',
  {\em Encyclopedia of environmetrics} .

\harvarditem[Elhorst et~al.]{Elhorst, Heijnen, Samarina \harvardand\
  Jacobs}{2017}{elhorst2017transitions}
Elhorst, J.~P., Heijnen, P., Samarina, A. \harvardand\ Jacobs, J.~P.
  \harvardyearleft 2017\harvardyearright , `Transitions at different moments in
  time: A spatial probit approach', {\em Journal of Applied Econometrics} {\bf
  32}(2),~422--439.

\harvarditem{Goddard}{2005}{goddard05}
Goddard, J.  \harvardyearleft 2005\harvardyearright , `Regression models for
  forecasting goals and match results in association football', {\em
  International Journal of forecasting} {\bf 21}(2),~331--340.

\harvarditem{Hughes \harvardand\ Franks}{2005}{hughes05}
Hughes, M. \harvardand\ Franks, I.  \harvardyearleft 2005\harvardyearright ,
  `Analysis of passing sequences, shots and goals in soccer', {\em Journal of
  Sports Sciences} {\bf 23}(5),~509--514.

\harvarditem[Jordet et~al.]{Jordet, Bloomfield \harvardand\
  Heijmerikx}{2013}{jordet13}
Jordet, G., Bloomfield, J. \harvardand\ Heijmerikx, J.  \harvardyearleft
  2013\harvardyearright , The hidden foundation of field vision in english
  premier league (epl) soccer players, {\em in} `Proceedings of the MIT Sloan
  Sports Analytics Conference'.

\harvarditem{King \harvardand\ Zeng}{2001}{king01}
King, G. \harvardand\ Zeng, L.  \harvardyearleft 2001\harvardyearright ,
  `Logistic regression in rare events data', {\em Political analysis} {\bf
  9}(2),~137--163.

\harvarditem{Klier \harvardand\ McMillen}{2008}{klier08}
Klier, T. \harvardand\ McMillen, D.~P.  \harvardyearleft 2008\harvardyearright
  , `Clustering of auto supplier plants in the united states: generalized
  method of moments spatial logit for large samples', {\em Journal of Business
  \& Economic Statistics} {\bf 26}(4),~460--471.

\harvarditem[Lasek et~al.]{Lasek, Szl{\'a}vik \harvardand\
  Bhulai}{2013}{lasek13}
Lasek, J., Szl{\'a}vik, Z. \harvardand\ Bhulai, S.  \harvardyearleft
  2013\harvardyearright , `The predictive power of ranking systems in
  association football', {\em International Journal of Applied Pattern
  Recognition} {\bf 1}(1),~27--46.

\harvarditem[Leitner et~al.]{Leitner, Zeileis \harvardand\
  Hornik}{2010}{leitner10}
Leitner, C., Zeileis, A. \harvardand\ Hornik, K.  \harvardyearleft
  2010\harvardyearright , `Forecasting sports tournaments by ratings of (prob)
  abilities: A comparison for the euro 2008', {\em International Journal of
  Forecasting} {\bf 26}(3),~471--481.

\harvarditem[Lucey et~al.]{Lucey, Bialkowski, Monfort, Carr \harvardand\
  Matthews}{2014}{lucey14}
Lucey, P., Bialkowski, A., Monfort, M., Carr, P. \harvardand\ Matthews, I.
  \harvardyearleft 2014\harvardyearright , Quality vs quantity: Improved shot
  prediction in soccer using strategic features from spatiotemporal data, {\em
  in} `Proc. 8th Annual MIT Sloan Sports Analytics Conference', pp.~1--9.

\harvarditem{Merkle}{2014}{Scoring}
Merkle, E.  \harvardyearleft 2014\harvardyearright , {\em scoring: Proper
  scoring rules}.
\newblock R package version 0.5-1.
\newline\harvardurl{https://CRAN.R-project.org/package=scoring}

\harvarditem{Merkle \harvardand\ Steyvers}{2013}{merkle13}
Merkle, E.~C. \harvardand\ Steyvers, M.  \harvardyearleft 2013\harvardyearright
  , `Choosing a strictly proper scoring rule', {\em Decision Analysis} {\bf
  10}(4),~292--304.

\harvarditem[Moura et~al.]{Moura, Santiago, Misuta, de~Barros \harvardand\
  Cunha}{2007}{moura07}
Moura, F.~A., Santiago, P. R.~P., Misuta, M.~S., de~Barros, R. M.~L.
  \harvardand\ Cunha, S.~A.  \harvardyearleft 2007\harvardyearright , Analysis
  of the shots to goal strategies of first division brazilian professional
  soccer teams, {\em in} `ISBS-Conference Proceedings Archive', Vol.~1.

\harvarditem[Sahu et~al.]{Sahu, Gelfand \harvardand\ Holland}{2006}{sahu06}
Sahu, S.~K., Gelfand, A.~E. \harvardand\ Holland, D.~M.  \harvardyearleft
  2006\harvardyearright , `Spatio-temporal modeling of fine particulate
  matter', {\em Journal of Agricultural, Biological, and Environmental
  Statistics} {\bf 11}(1),~61--86.

\harvarditem{Williams \harvardand\ Walters}{2011}{williams11}
Williams, T. \harvardand\ Walters, C.  \harvardyearleft 2011\harvardyearright ,
  The effects of altitude on soccer match outcomes, {\em in} `Proceedings of
  the MIT Sloan Sports Analytics Conference'.

\end{thebibliography}
