Data Collection

The dataset examined is an English Premier League (EPL) 2021-2022 season database. The dataset came from Kaggle.com and is a CSV file containing 106 variables. It contains general soccer statistics (goals, fouls, shots, cards, etc.) and has them divided between the home team and the away team. There is a total of 380 observations, with one row dedicated to every game played throughout the season. The data itself comes from readily available game statistics and they are put into the database. Some of the variables also contained analytics for sports betting, but they were omitted for the analysis.

Research Question

The goal is to find a model that can find the most likely predictors of the total number of goals in a game when looking at combined statistics between the home and away team. A larger number of goals in a game tends to boost excitement amongst fans and can create more thrilling matches. The other variables of interest that will be compared to the total number of goals are fouls, shots, shots on target, cards, and corners. Intuitively, some predictors make sense to have a correlation to the number of goals scored, such as shots and shots on target. Corners is likely to have some relationship to goals scored, but may not be significant, and cards and fouls are difficult to predict.

The first part of the analysis will be creating a visualization of the data. For this, the variables must be combined into cumulative statistics as opposed to being the individual teams. This can be done very easily as such:

mydata$goals <- mydata$FTHG + mydata$FTAG
mydata$fouls <- mydata$HF + mydata$AF
mydata$shots <- mydata$HS + mydata$AS
mydata$corners <- mydata$HC + mydata$AC
mydata$cards <- mydata$HY + mydata$AY + mydata$HR + mydata$AR
mydata$shots_OT <- mydata$HST + mydata$AST

Now that the data is in the correct form, we can create a pairwise plot to see relationships between the different variables as they relate to the total number of goals.

pairs(~goals + fouls + shots + corners + cards + shots_OT, data = mydata, 
      cex = 0.4, oma = c(2, 2, 2, 2), mar = c(4, 4, 2, 2))

Much of the data is very discrete, and this makes it difficult to see if there are clear relationships. Moving from left to right, we can observe the relationships of variables in relation to goals.

One note about the data being discrete is that there are likely duplicates of certain data points, which can make the true form of the data difficult to observe by eye. It may be, for example, that there are many duplicates of high numbers of corners with respect to a high number of goals scored, but there is no way to tell.

Weighted Least Squares Regression (WLS)

The method examined on the data mentioned above is Weighted Least Squares Regression, also known as WLS. It is very similar to Ordinary Least Squares Regression (OLS), but has addition of weights. Such weights are given to points with respect to their variance, and the weight determines the influence of the point. Lower variance correlates to a higher weight and these points will have a greater influence on the fitting of the model. The addition of the weights should create a “better” model that more reasonably represents the data.

Ordinary least squares regression (OLS) assumes that there is constant variance in the errors (homoscedastic), whereas weighted least squares regression can be used when there is not constant variance in the errors (heteroscedastic). The model has one response variable, as well as \(p\) predictor variables, and follows the form \[Y = \beta_0 + \beta_1X_{i1} + \dots + \beta_{p-1}X_{i,p-1} + \epsilon_i\] where the errors are independent \(\epsilon \sim N(0,\sigma_i^2)\). The notion of least squares regression is that all of the observations are given the same influence over the model, despite their relationship. This means that outliers are given the same amount of influence as points clustered around the line of fit, and may not always be the best approach. It is also assumed for WLS that there is linearity with the data, the observations are independent of one another, the errors are normally distributed, the errors are centered at zero, and there is heteroscedasticity.

The weights for WLS can be calculated through the reciprocal of the variance \[w_i = 1/\sigma_i^2\] where a larger variance results in a smaller weight. This gives a matrix \(W\) of the weights represented as \[ W = \begin{bmatrix} w_1 & 0 & 0 & \cdots & 0 \\ 0 & w_2 & 0 & \cdots & 0 \\ 0 & 0 & w_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & w_n \end{bmatrix} \]

that has a weight for each observation.

The goal is to minimize the sum of squares \[ \sum_{i=1}^{n} w_i \left(y_i - \beta_0 - \beta_1 x_{i1} - \cdots - \beta_j x_{ij}\right)^2 \]

To create an OLS model to compare to WLS, the lm function to create the linear model using the variables we have specified can be used.

epl_lm = lm(goals ~ fouls + shots + corners + shots_OT + cards, data = mydata)
summary(epl_lm)
## 
## Call:
## lm(formula = goals ~ fouls + shots + corners + shots_OT + cards, 
##     data = mydata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2988 -0.8599 -0.0893  0.8605  4.2226 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.01729    0.48075   4.196  3.4e-05 ***
## fouls       -0.02339    0.01413  -1.655  0.09872 .  
## shots       -0.03047    0.01561  -1.952  0.05168 .  
## corners     -0.06248    0.02109  -2.963  0.00324 ** 
## shots_OT     0.31556    0.02604  12.119  < 2e-16 ***
## cards       -0.02248    0.03783  -0.594  0.55281    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.332 on 374 degrees of freedom
## Multiple R-squared:  0.3385, Adjusted R-squared:  0.3296 
## F-statistic: 38.27 on 5 and 374 DF,  p-value: < 2.2e-16

The resulting model is \[\widehat{Goals_{OLS}} = 2.02 - 0.023 \cdot Fouls - 0.03 \cdot Shots - 0.06 \cdot Corners + 0.32 \cdot ShotsOT - 0.02 \cdot Cards\] It is interesting to see that all variables but shots on target decreased the number of goals scored. This is most surprising for shots and corners, as even though they have very small coefficients, this goes against the predictions.

Now that the OLS model has been created, the model should be checked for heteroscedasticity. The Breusch-Pagan test can be used to determine heteroscedasticity. The test investigates the variances of the residuals to determine if there is or is not constant variance. The null hypothesis (\(H_0\)) is that homoscedasticity is present (the residuals have equal variance). The alternative hypothesis (\(H_A\)) is that heteroscedasticity is present (the residuals do not have equal variance). If the null hypothesis is rejected, then conditions for using WLS have been met to use the weights. The Breusch-Pagan test is run through the lmtest library with the bptest method. The arguments for the test are a fitted OLS model, and this has been defined above as epl_lm.

bptest(epl_lm)
## 
##  studentized Breusch-Pagan test
## 
## data:  epl_lm
## BP = 11.84, df = 5, p-value = 0.03705

We can see the resulting p-value is 0.037, so we can reject the null hypothesis at a significance level of \(\alpha = 0.05\) and assume the data is heteroscedastic, which allows us to use weighted least squares.

We can manually create the weights using the residuals from the OLS model, and assign those to a new model called wls_model.

resids <- abs(residuals(epl_lm))
resids_sq <- resids^2
variance_model <- lm(resids_sq ~ fouls + shots + corners + 
                  shots_OT + cards, data = mydata)
predicted_var <- predict(variance_model)
mydata$predicted_var <- predicted_var
wls_model <- lm(goals ~ fouls + shots + corners + 
             shots_OT + cards, data = mydata, 
             weights = 1 / predicted_var)

summary(wls_model)
## 
## Call:
## lm(formula = goals ~ fouls + shots + corners + shots_OT + cards, 
##     data = mydata, weights = 1/predicted_var)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.31714 -0.67060 -0.08344  0.68470  3.02414 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.96729    0.46515   4.229 2.95e-05 ***
## fouls       -0.02256    0.01293  -1.745  0.08173 .  
## shots       -0.03152    0.01481  -2.129  0.03391 *  
## corners     -0.06455    0.01999  -3.230  0.00135 ** 
## shots_OT     0.32257    0.02484  12.988  < 2e-16 ***
## cards       -0.01679    0.03638  -0.461  0.64475    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.004 on 374 degrees of freedom
## Multiple R-squared:  0.3698, Adjusted R-squared:  0.3614 
## F-statistic: 43.89 on 5 and 374 DF,  p-value: < 2.2e-16

Working from top to bottom, the code is going through the steps of creating the weights from the residuals of the OLS model, and then assigning them to a new model for WLS. The predict function estimates the variance based on the OLS model. This gives a model equation of \[\widehat{Goals_{WLS}} = 1.97 - 0.023 \cdot Fouls - 0.03 \cdot Shots - 0.06 \cdot Corners + 0.32 \cdot ShotsOT - 0.02 \cdot Cards\] which is slightly different from that of the OLS model.

Results

To check the model conditions, we can make residual vs. fitted plots and QQ plots. Heteroscedasticity is already confirmed from the Breusch-Pagan test, but seeing patterns visually will assist with confirmation.

par(cex = 0.6, cex.axis = 0.6, cex.lab = 0.7, cex.main = 0.8)
plot(epl_lm, which = 1, main = "OLS")

plot(epl_lm, which = 2, main = "OLS")

Looking at the residuals vs. fitted plot, we can see that the errors are roughly centered around zero, but there is a pattern of lines that suggests a lack of homoscedasticity. The observations may be independent and appear to be linear, and in the QQ plot it can be seen that there are not any major concerns. There is some flaring at both ends of the graph, but nothing significant. Despite the fact that the errors meet some of the conditions, a lack of homoscedasticity renders the model conditions unfulfilled, and a new model should be used. Let’s take a look at the WLS plots.

par(cex = 0.6, cex.axis = 0.6, cex.lab = 0.7, cex.main = 0.8)
plot(wls_model, which = 1, main = "WLS")

plot(wls_model, which = 2, main = "WLS")

We can see that these plots exhibit minimal changes when compared to the OLS plots. The residuals vs. fitted plot is nearly identical to that of the OLS residuals vs. fitted plot. However, the patterns in this case are not a concern, as heteroscedasticity is a condition of WLS. The other conditions of linearity, zero-centered errors, and independence of errors appear to be validated. The QQ plot does exhibit some very slightly different behavior. The flaring on the bottom left is about the same, but the deviation on the upper right end is slightly less than the OLS QQ plot. We can see that the end of the points comes back towards the line, which did not happen prior to the introduction of the weights. It is hard to tell if there is any noticeable improvement by WLS on these plots, but the summaries of the models can be examined to find other differences.

## 
## Call:
## lm(formula = goals ~ fouls + shots + corners + shots_OT + cards, 
##     data = mydata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2988 -0.8599 -0.0893  0.8605  4.2226 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.01729    0.48075   4.196  3.4e-05 ***
## fouls       -0.02339    0.01413  -1.655  0.09872 .  
## shots       -0.03047    0.01561  -1.952  0.05168 .  
## corners     -0.06248    0.02109  -2.963  0.00324 ** 
## shots_OT     0.31556    0.02604  12.119  < 2e-16 ***
## cards       -0.02248    0.03783  -0.594  0.55281    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.332 on 374 degrees of freedom
## Multiple R-squared:  0.3385, Adjusted R-squared:  0.3296 
## F-statistic: 38.27 on 5 and 374 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = goals ~ fouls + shots + corners + shots_OT + cards, 
##     data = mydata, weights = 1/predicted_var)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.31714 -0.67060 -0.08344  0.68470  3.02414 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.96729    0.46515   4.229 2.95e-05 ***
## fouls       -0.02256    0.01293  -1.745  0.08173 .  
## shots       -0.03152    0.01481  -2.129  0.03391 *  
## corners     -0.06455    0.01999  -3.230  0.00135 ** 
## shots_OT     0.32257    0.02484  12.988  < 2e-16 ***
## cards       -0.01679    0.03638  -0.461  0.64475    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.004 on 374 degrees of freedom
## Multiple R-squared:  0.3698, Adjusted R-squared:  0.3614 
## F-statistic: 43.89 on 5 and 374 DF,  p-value: < 2.2e-16

There are some elements of the summary for the WLS model that suggest it created a better fit for the data. First, the residual standard error drops from 1.332 to 1.004. The adjusted \(R^2\) increased from 0.3296 to 0.3614, and the F-statistic went from 38.27 to 43.89. Additionally, the p-value of cards, which was not significant in either model, increased 0.55 to 0.64, suggesting that the model gave more value to the predictors with higher significance and cards became less important. All of these results are expected, since the data being heteroscedastic means that a WLS model will give more promising results to that of an OLS model.

Discussion and Critique

Shots, corners, and shots on target were significant predictors, which was expected. If a team is closer to the goal, they are statistically more likely to score. However, there is a very slight negative correlation between shots and corners which is surprising. Fouls and cards were not found to be significant predictors. The absence of fouls as a significant predictor is interesting, as fouls towards an opponent’s goal can lead to a free kick or a penalty kick inside the 18-yard box, which almost always results in a goal. However, there are also matches where there are a large number of fouls with many of them towards the center of the field. Since more fouls leads to more cards, it also makes sense that there would not be a correlation with number of cards and number of goals. Defenders may tend to be more careful when they are closer to their own goal, and that could reduce the number of fouls near the goal that could in turn result in goal-scoring opportunities. The WLS model appeared to also create a better model compared to the OLS model as expected.

There are some drawbacks and difficulties of using the WLS model. One note is that the weights can be calculated by ways other than the reciprocal of the variance, and this brings on the point of how to choose weights. Should variance be prioritized? Or should something else? Regardless, there is interpretation to be had about the choice of weights. Another drawback of the WLS model is that the differences can be very hard to spot. With the residuals vs. fitted plot and the QQ plot, it is very difficult to visually see the differences between the methods. The data also needs to be heteroscedastic, which is not always the case with datasets. The weights do not have much use if the data has equal variance amongst the residuals.

These findings are from the entire 2021-2022 season of the English Premier League, and may be similar to findings to that of other seasons. Understanding the statistics of what can lead to a larger number of goals in a game is something that can be examined by rule makers to make games more exciting, as well as for individual teams when looking at their own performances.

Sources