# Build a predictive model for number of hockey goals in a season for an individual player

I'm looking to build a predictive model for hockey players individual statistics. My goal is to predict how many points a player can be expected to have at the end of the season.

To do so, I thought I could base myself on the progression of the player's statistics over the year. I currently have 10 years worth of data and, for example, a player could have played in 5 of those 10 years. In this case, I'm able to see if this particular player is on an increasing or decreasing curve regarding his statistics.

So, to resume, for some thousands of players the data goes like this:

### Player #1

Season 2005-2006, Points : X, Game played : y, ...

Season 2006-2007, Points : X, Game played: Y, ...

Season 2007-2008, Points : X, Game played: Y, ...

...

### Player #2

...

My main problem is that I don't know what kind of data this is so I can't find more information on how to build my model. In my title I've written "multiple small multivariate time series" for a lack of better name. I've obviously looked into multivariate time series to fit them in an RNN architecture, but it doesn't seem to fit because of the multiple time series that aren't really related.

I've found this question which kind of relates but there were no answer and the related comments didn't help me.

All in all, I'm trying to find some guidance as to what kind of model I could use to work with my data. I'd also like to know if there's a name for this kind of data, so I can inform myself better.

"Predict how many points a player can be expected to have at the end of the season" could be framed as count-based regression problem (not time series).

A player will have a collection of features (e.g., age, position, team, minutes plays). A machine learning model can learn how to weight to the features to make a prediction. In this case, the prediction is the number of goals by the end of the season (e.g., 0, 1, 2, 3, …).

The data and features might have time component but the model does not have to be time-series.

• I agree with this, I don't know about hockey but I would assume that the performance of players in a match are independent on other matches? Or maybe if a player is in form then they will continue being in form? If that is the case maybe the time component needs to be taken into account, or at least some way of capturing this, through feature augmentation? Mar 15 at 13:31

This can be done fairly well with a Bayesian model. Let's first make an unrealistic and simplifying solution:

• Games are IID (False, poor performance may result in being moved to the farm team)
• There is no variation in the expected goals per game by position (again False, defensemen are expected to score fewer goals than forwards).

We will revisit the second assumption later on.

Consider a scenario in which you only have data from the last season. From your post, it looks like you have the following information:

• Player identifier $$i$$
• Games played $$n_i$$
• Goals scored $$y_i$$

You can estimate the number of goals per game for player $$i$$ as $$\lambda_i$$. Naturally, those players with fewer games played will have a noisier estimate than those who played every game. We can regularize those estimates by specifying a hierarchical model. Let's simulate some data to exposit the approach.

library(tidyverse)
n_player <- 100

# Simulate 100 player's goal rate
# The average goal per game is 0.4
lam <- rgamma(n_player, 1, scale=0.4)

# Simulate how many games they played
n_games <- round(runif(n_player, 1, 80))

# The goals per game can be considered as a poisson rate
# if games are iid, then simply muliply the rate by the number of games
n_goals_per_season <- rpois(n_player, n_games*lam)

d <- tibble(lam, n_goals_per_season, n_games, player = 1:n_player)


## Fitting a Model

Now, we specify a Bayesian model for these data. We will assume that the $$\lambda_i$$ are draws from some larger population. Let's use the brms library to specify this model, as we will need it later. The model we fit is

$$y_i \mid \lambda_i \sim \mbox{Poisson}(\lambda_i; n_i)$$ $$\log(\lambda_i) = \beta_{0,i} + \log(n_i)$$

$$\beta_{0, i} \sim \mbox{Normal}(., .)$$

Here $$\beta_{0, i}$$ are the log expected goals per game, and they are assumed to be draws from a normal distribution with some mean and variance. We get to specify that in our model. For now, I'll just use brms' defaults since I don't have good priors on what these should be. The model in brms is

library(brms)

fit <- brm(n_goals_per_season ~ 1 + (1|player) + offset(log(n_games)),
data=d,
family = poisson(),
backend = 'cmdstanr')


Let's plot the fitted vs estimated rates from the raw data per player. I will plot them on the log scale so we can see some interesting patterns and facet by the number of games played.

d %>%
ggplot(aes(log(n_goals_per_season/n_games), log(.epred/n_games))) +
stat_pointinterval() +
geom_abline() +
labs(
x=expression(log(y) - log(n)),
y = expression(log(hat(lambda)))
) +
facet_wrap(~cut(n_games, c(0, 10, 20, 40, 80)))



Let's take a look at the top left facet -- those players playing between 1 and 10 games. Note how these estimates are regularized. Because we expect the lowest and highest average scorers to score more/less than we saw because we saw them in so few games. There is also some regularization in the other facets for the lowest scoring players.

Ok, let's talk about how to make predictions.

## Making Predictions

So we have an estimated goals per game $$\hat{\lambda}_i$$. How do we make predictions. Well, first we need some idea of how many games each player will play at the end of the season. This in itself may require a model since injuries and promotions/demotions from farm teams can modulate this. But let's say we are going to focus on the first 5 games of play. All we need to do then is draw samples from a poisson distribution with rate $$5\hat{\lambda}_i$$. Because our model is Bayesian, we get a distribution of expected goals. Luckily tidybayes has some functions to make this easy. I'll demonstrate how to do this with one player.


tibble(player=1, n_games=5) %>%
ggplot(aes(.prediction)) +
geom_bar(aes(y = after_stat(count)/sum(after_stat(count)))) +
labs(y = 'Probability', x = 'Goals Scored After 5 Games')


This player is expected to probably score 0 or 1 total goals in the first 5 games of the next season.

But your question was doing this over time. How do we extend this model?

## Model Extensions

You have player data over seasons, not just for one season. We can extend this data so that the goals per game can evolve over time $$t$$.

$$y_i \mid \lambda_i \sim \mbox{Poisson}(\lambda_i; n_i)$$ $$\log(\lambda_i) = \beta_{0,i} + \beta_{1,i} t + \log(n_{i, t})$$

This is essentially the same model, but now there is a linear effect of time (so players can get better and score more...or less).

It is also easy to account for player position by adding appropriate covariates to the model. For an additive effect of position $$p$$, maybe we could specify

$$y_i \mid \lambda_i \sim \mbox{Poisson}(\lambda_i; n_i)$$ $$\log(\lambda_i) = \beta_{0,i} + \beta_{1,i} t +\beta_{2, i}p + \log(n_{i, t})$$

with appropriate priors on all parameters. Or maybe the rate at which players improve depends on their position (interaction). Or maybe there are team level effects, or shift level effects, etc etc.

there are many extensions to this model, but starting with a fairly simple model (predicting the next season from the last) is a good way to start. Then, we can extend that model slowly to the one we actually want.