STAT542 ISP Final Project: Bayesian Analysis of NHL Goalie Performance
By: Matt Connor & Jeff Ellington
5/5/2017
1. Introduction: Questions arising from 2015-2016 Stanley Cup Playoffs on NHL Goaltending
2016 marked the first Stanley Cup appearance in the San Jose Sharks’ 25 year history. For me, a die-hard Sharks fan born in San Jose in 1987, this marked another Sharks ‘first’: my first Sharks game at the Cow Palance in 1991, first ever Sharks game in San Jose in 1993, and now first ever Sharks Stanley Cup home game: Game 3 of the 2016 Stanley Cup Finals against the Pittsburgh Penguins.
There’s a common quip that goaltending wins championships. In 2016 the Sharks seemingly had an edge, with Martin Jones registering Top 5 in Shutouts, Wins, and Games Played during the regular season. In comparison, Pittsburgh’s rookie goaltender Matt Murray had started just 13 regular season games. He replaced veteran Marc-Andre Fleury who suffered a concussion with 5 games remaining in the regular season. Yet, Pittsburg eventually won the Stanley Cup. Matt Murray played a critical role, holding the Sharks to 1 goal in 3 of Pittsburgh’s 4 Stanley Cup wins.
Looking back at Pittsburgh’s 2015-2016 season is a good microcosm of NHL goaltending and the decisions surrounding it. For quick background, Marc-Andre Fleury had already won the Stanley Cup as Pittsburgh’s starting goalie in 2008-2009. In addition, he’d been selected as a 2x All-Star (2011 & 2015). During the 2015-2016 season, Fleury suffered 2 concussions and as a result played in only 58 of 82 games. The remaining games were split by tenured backup Jeff Zatkoff (11 starts/14 appearances) and rookie Matt Murray (13 starts/13 appearances). Jeff Zatkoff started the first 2 games of the playoffs, going 1-1. Murray replaced Zatkoff, starting 21 of the final 22 games. Fleury saw just 1 additional game in the Eastern Conference Semifinals, where he let in 4 goals in a OT loss to Tampa Bay before being replaced by Murray in the next game.
The natural question that arises: can goalies get “hot”? For his 13 regular seasons games, Murray had the 2nd highest SV% and lowest GAA of all goalies with over 12 games. In retrospect, it seems Pittsburgh Coach Mike Sullivan was wise to switch from Zatkoff to Murray and ride Murray’s hot streak through the playoffs even after Fleury was healthy.
For the remainder of the paper, we will use Bayesian methods to explore the following questions: do goalies have “hot” streaks? How does age effect performance? Are there truly “elite” goalies? How much does a winning team impact a goalie’s performance? We will take a model based approach using publicly available data. The model will be composed of 4 levels. We will model the average number of goals allowed per minute for the season as binomial, rate of goals per minute as a function of team, wins, age, and save percentage, while separting for “elite” vs. “non-elite” goalies using expectation-maximization, and finally use a Hidden Markov Model to identify “hot” or “cold” based on prior years.
Our goal is see if the data supports our belief that Pittsburgh Coach Mike Sullivan’s intuitive decision to ride Matt Murray’s hot streak during the 2015-16 playoffs was correct, despite having his starting, All-Star caliber goalie healthy on the bench. The ability to analyze situations and make these goaltending decisions would be instrumental to both NHL coaches and general managers.
2. Data Sources and Model Selection
2.1 Data and Summary
Our data comes from scraping the publicly available Hockey-Reference website (http://www.hockey-reference.com/). It contains information on every NHL player from 1922 through today. We will be using data from the past 20 seasons (1996 to 2016) for our model estimation to predict on 2017 season goalie data. Our analysis will focus exclusively on NHL goalies. From this 20 year period, we will examine 1,287 goalie-seaons (summary stats of a goalie’s total season). We’ve limited the minimum number of minutes in a goalie-season at 400 to prevent goalies with very small sample sizes from polluting the data. For 2017, there will be 64 goalie-seasons to evaluate our model against. Within each season j, we will use the following predictors for goalie i:
1. Goals Against Total: Y(ij)
2. Total Minutes: M(ij)
3. Age: A(ij)
4. Team: T(ij)
5. Win Percentage: W(ij)
6. Save Percentage: S(ij)The two 2016 Stanley Cup goalies, Matt Murray (PIT) and Martin Jones (SJS) had Y(ij) = 25 & 143 goals against in M(ij) = 749 & 3,786 Total Minutes, respectively. Player ages ranged from 18 to 43, with most plyers between 21 and 38. There were 35 different possible NHL teams during the 20 year window, with an additional level “TOT” for players who were traded during a season. Win percentage ranges from 0% to 83% and appears to be a normal distribution centered at 41%. Save percentage ranges from 83% to 94% with a median of 90.8%. There is a large clustering around 90-92%.
2.2 Exploratory Data Analysis
If we were looking at the data from an empirical Bayesian perspective, we’d note that the average goals per minute (Yij/Mij) for the 1,287 goalie-seasons is 0.0452 goals/minute. The distribution appears to have an underlying normal structure with some type of clustering to the left of center. Based on previous analyses, we’d expect there might be two underlying groups, elite goals and non-elite goalies that could be distoring the distribution. We will examine this later in the modeling phase.
hist(gs.train$GPM, main = "Fig 2.2.1: Histogram of Goals Per Minute", xlab = "Goals Per Minute", col = "turquoise4")
The 1,287 goalie-seasons are comprised of 238 individual goalies, who have ranges of seasons played from 1 to 18. In 2017, 5 new goalies registering over 400 minutes entered the NHL. Comparing the histograms of seaons (Fig 2.2.2) and age (Fig 2.2.3), we can see that they follow a similar shape on the right side. On the left, younger players are effectively truncated into the lower number of total seasons (0-5), which makes sense given not every young player will enter the league at exactly the same age.
par(mfrow=c(2,1))
hist(table(gs.train$Player), main = "Fig 2.2.2: Histogram of Goalie-Seasons by Player", xlab = "# of Seasons", col = "turquoise4")
hist(gs.train$Age, main = "Fig 2.2.3: Histogram of Goalie Age", xlab = "Age", col = "turquoise4")
2.3 Expectation-Maximization
In Section 2.1 we mentioned that save percentage seemed to have a large clustering around 90-92%, distorting what would otherwise look like a normal distribution. Our hunch was that this implied there were multiple classes of goalies, some “elite” with higher save percentages and some “non-elite” with lower save percentages. We used a normal-normal Expectation-Maximization function to pull apart the two groups. The results can be seen below in Fig 2.3.1. What surprised us most, was that instead of a small elite group, there really appears to be a large set of solid goalies (we are still calling them “elite”) and a smaller set of pretty poor goalies. It seems that goalies who are unable to have save percentages over 90% do not last long in the NHL. We will try to pull out more signal in our complete Bayesian model, but it may have been interesting to try to tease this into 3 sub-groupings.
finalparam<-itermat[201,]
alpha <- finalparam[1]
mu0 <- finalparam[2]
mu1 <- finalparam[3]
sigsq0 <- finalparam[4]
sigsq1 <- finalparam[5]
par(mfrow=c(1,1))
hist(sv.per.sa,prob=T, main = "Fig 2.2.1: Histogram of Save %", xlab = "Save Percentage", col = "turquoise4")
x <- ppoints(10000)
y1 <- (1-alpha)*dnorm(x,mu0,sqrt(sigsq0))
y2 <- alpha*dnorm(x,mu1,sqrt(sigsq1))
lines(x,y1,col="orangered")lines(x,y2,col="black")
legend('topleft',c('Non-Elite','Elite'),col=c("orangered","black"),lty=1, cex=0.75)
Doing a similar normal-normal EM analysis for Shutout % (for goalies with at least 1 shutout in a season) hints at a smaller group of elite goalies. As can be seen in Fig 2.3.2 below, few goalies are able to achieve shutout percentages above 10%. The majority of goalies are clustered between 0-7.5%, with a tail of elite goalies extending out toward 20%. Hall-of-Famer Dominik Hasek recorded the highest single season shutout total in our dataset, 13 shutouts at a rate of 18.05%. Furthermore, future Hall-of-Fame goalie Martin Brodeur logged 7 seasons with 9 or more shutouts (out of his All-time NHL leading 125 total shutouts). Over those seasons he averaged a 13.67% shutout rate. Not surprisingly given the anecdotal accounts above, the histogram of elite vs. non-elite goalies by Shutout Percentage seem to sub-set a smaller number of elite goalies.
hist(y,prob=T,main="Fig 2.3.3: Histogram of Shutout %", xlab = "Shutout %", col = "turquoise4")
x <- seq(0,0.3,len=1001)
notelite.dens.gamma <- (theta.norm.gamma[1])*dgamma(x,theta.norm.gamma[2],theta.norm.gamma[3])
elite.dens2 <- (1 - theta.norm.gamma[1])*dnorm(x,theta.norm.gamma[4],sqrt(theta.norm.gamma[5]))
lines(x,notelite.dens.gamma,col="black",lwd=2)lines(x,elite.dens2,col="orangered",lwd=2)
legend('topright',c('Non-Elite','Elite'),col=c("orangered","black"),lty=1, cex=0.75)
mb.soNow we turn to our response, goals per minute (GPM). Again running a normal-normal EM function, we can see a split into two groups in Fig 2.3.4. Like Save %, GPM seems to split into a large “elite” group (more like above average) and a smaller group of poor performers. While this helps describe the distribution shape, there is likely more we can do building a complete Bayesian hierarchical model that will forecast performance taking into account hot and cold streaks via additional model parameters.
finalparam<-itermat[201,]
alpha <- finalparam[1]
mu0 <- finalparam[2]
mu1 <- finalparam[3]
sigsq0 <- finalparam[4]
sigsq1 <- finalparam[5]
par(mfrow=c(1,1))
hist(gpm,prob=T, main = "Fig 2.3.4: Histogram of Goals per Min", xlab = "Goals Per Minute", col = "turquoise4")
x <- ppoints(10000)
y1 <- (1-alpha)*dnorm(x,mu0,sqrt(sigsq0))
y2 <- alpha*dnorm(x,mu1,sqrt(sigsq1))
lines(x,y1,col="black")lines(x,y2,col="orangered")
legend('topright',c('Non-Elite','Elite'),col=c("orangered","black"),lty=1, cex=0.75)
2.4 Revisiting Matt Murray and Martin Jones
Before turning to our hierarchical model, it’s worth re-visiting last year’s two Stanley Cup Final’s goalies: Matt Murray (PIT) and Martin Jones (SJS). Both players had a 90% chance of being elite. Furthermore, Martin Jones has had an 87-90% of being elite from 2014-2016.
finalparam <- itermat[201,]
prob.elite <- 1-Estep(gpm, finalparam[1], finalparam[2], finalparam[3], finalparam[4], finalparam[5])
plot(gpm,prob.elite,main="Prob of Elite Goalie")
## Matt Murray
mm.seasons <- gs.train[gs.train$Player == "Matthew Murray", c('Player', 'year')]
mm.elite <- prob.elite[gs.train$Player == "Matthew Murray"]
mm.gpm <- gpm[gs.train$Player == "Matthew Murray"]
final.mm <- data.frame(mm.seasons, mm.elite, mm.gpm)
final.mm## Martin Jones
mj.seasons <- gs.train[gs.train$Player == "Martin Jones", c('Player', 'year')]
mj.elite <- prob.elite[gs.train$Player == "Martin Jones"]
mj.gpm <- gpm[gs.train$Player == "Martin Jones"]
final.mj <- data.frame(mj.seasons, mj.elite, mj.gpm)
final.mj