Massey ratings

This web page outlines a technique based in linear algebra to rate teams playing games in a common league based on the score differences that we see in the games. The basic idea is to perform a linear regression where the input consists of games and the output is the score differences.

This technique dates back to Ken Massey’s 1997 senior honors thesis at Bluefield College. My understanding of the technique is based mostly on Chapter 2 of Langville and Meyer’s text on sport’s ranking.

Code

This webpage describes the mathematics of Massey’s algorithm and also shows how to implement it in Python. Each code block has a clipboard icon in the upper right that you can click on to copy the code. You should be able to run the code in a scientific Python environment.

Imports

Let’s go ahead and import the libraries that we’ll need straight away:

# Data manipulation
import pandas as pd

# Numerics
import numpy as np

# SciPy's linear system solver
from scipy.linalg import solve

Description of the algorithm

Let’s describe the algorithm and take a look at a relatively simple example

The general problem

In general, we’ve got \(m\) games played between \(n\) teams. Let’s suppose that in the \(k^{\text{th}}\) game, team \(i\) beat team \(j\) by \(y_k\) points. We’d like the teams ratings \(r_i\) and \(r_j\) to reflect this fact, i.e. we’d like \[r_i - r_j = y_k.\] To put it another way, we’d like our ratings of teams \(i\) and \(j\) to indicate the expected score difference between the teams, if they play again.

If \(\vec{r}\) is the vector of ratings and \(\vec{y}\) is the vector of point differences, then we can express this as the matrix equation \[X\vec{r} = \vec{y},\] where \(X_{ki} = 1\), \(X_{kj} = -1\), and all other entries are zero.

Now, there are generally many more games than teams so that the system we’ve written down is likely overdetermined. Nonetheless, we can attempt to find an approximate solution in the least squares sense and an efficient way to do so is to solve the so-called normal equations obtained by multiplying both sides by the transpose of \(X\): \[X^T X\vec{r} = X^T \vec{y},\]

The matrix \(M = X^T X\) formed in this process is called the Massey matrix and has some interesting properties. Note that each entry in column \(i\) of the matrix \(X\) is either \(1\) (if team \(i\) won that game), \(-1\) (if team \(i\) lost that game), or zero (if team \(i\) didn’t play in that game). As a result, \(M_{ii}\) (which is obtained by computing the dot product of the \(i^{\text{th}}\) column of \(X\) with itself) denotes the number of games played by team \(i\). Similarly, \(M_{ij}\) is the opposite of the number of games played by teams \(i\) and \(j\). Thus, the sum of each row is zero. Put another way, the columns sum to the zero vector and are therefore linearly dependent. This implies that \(M\) is singular so that the normal equations don’t have a unique solution.

A way to get around this problem is to set all the elements in some row of \(M\) equal to 1. This is conventionally done in the last row. We similarly replace the last element of \[\vec{p} = X^T y\] to zero. This implies the constraint that the sum of the ratings must be zero. Furthermore, the resulting system \[M \vec{r} = \vec{p}\] can be shown to have a unique solution.

A simple working example

As a simple example, let’s suppose that we have three teams labeled \(A\), \(B\), and \(C\) and that each team plays the other teams twice each. We might visualize the collection of games via a labeled, directed graph like so:

Note that a directed edge from team \(i\) to team \(j\) indicates that team \(i\) defeated team \(j\) in that game; the edge labels indicates the score differences. Given the directed graph shown above, how would you rank the teams?

To rank the teams the Massey way, we start with the general equation \[X \vec{r} = \vec{y}\] which takes the specific form

\[ \begin{pmatrix}1&-1&0\\1&-1&0\\0&1&-1\\0&-1&1\\1&0&-1\\-1&0&1\end{pmatrix} \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}3\\2\\1\\2\\3\\1\end{pmatrix}. \]

If we multiply both sides by \(X^T\) and modify the last rows according Massey’s strategy, the equation becomes

\[ \begin{pmatrix}4&-2&-2\\-2&4&-2\\1&1&1\end{pmatrix} \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}7\\-6\\0\end{pmatrix}. \]

The solution of this system is

\[ \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}7/6\\-1\\-1/6\end{pmatrix}. \]

In particular, note that \(r_1 < r_2 < r_0\), which is probably in accordance with the order you expected.

The Massey equation

It’s worth noting that the Massey equation \[M\vec{r} = \vec{p}\] can itself be written down and interpreted without reference to \(X\) or \(y\). As we’ve already noted, \(M_{ii}\) denotes the number of games played by team \(i\) and \(M_{ij}\) is the opposite of the number of games played by teams \(i\) and \(j\). Similarly, the entries of \(\vec{p}\) is the accumulated set of score differences for each team.

For example, if we expand the first term in the product \(M\vec{r}\) and set that equal to the first term of \(p\), we get \[4r_0 - 2r_1 - 2r_2 = 7.\] This indicates that team \(0\) played four games, two against each other team, and that it outscored its opponents by a combined score of seven. If we solve for \(r_0\), we get \[r_0 = \frac{1}{2}r_1 + \frac{1}{2}r_2 + \frac{7}{2}.\] In general, the rating of team \(i\) is a combination of the ratings of the other teams together with the number of points that team \(i\) outscored those other teams by.

Implementation

Now, we’ll implement the technique in code for both the simple example and a more realistic example.

The simple example

Here’s an implementation of the technique for the simple example described above in Python:

X = np.array([
    [1, -1, 0],
    [1, -1, 0],
    [0, 1, -1],
    [0, -1, 1],
    [1, 0, -1],
    [-1, 0, 1]
  ])
M = X.transpose().dot(X)
M[-1] = np.ones(len(M))

y = np.array([[3], [2], [1], [2], [3], [1]])
y = np.array([3, 2, 1, 2, 3, 1])
y = X.transpose().dot(y)
y[-1]=0

solve(M,y)
array([ 1.16666667, -1.        , -0.16666667])

This is exactly the [7/6,-1,-1/6] ratings that we mentioned above.

The Big South

Let’s take a look at a somewhat bigger example, namely the current Big South. I happen to have a CSV file on my web space that has the scores of all Big South games played this season up through last night (Thur, Feb 13). We can download it right off my web space and take a look:

games = pd.read_csv('https://marksmath.org/data/BigSouth2025.csv')  
m = len(games)
games.head(8)
WTeamName WScore LTeamName LScore
0 Winthrop 95 SC_Upstate 76
1 Presbyterian 68 Longwood 60
2 High_Point 76 Radford 58
3 Charleston_So 72 Gardner_Webb 63
4 Radford 87 Winthrop 67
5 UNC_Asheville 103 High_Point 99
6 Gardner_Webb 63 Presbyterian 61
7 Longwood 83 Charleston_So 78

We’re going to need to set up a matrix indexed by the teams, so let’s assign an index to each team.

wteams = games.WTeamName
lteams = games.LTeamName
teams = pd.concat([wteams, lteams]).drop_duplicates().sort_values()
team_to_idx_dict = dict([(v,i) for i,v in enumerate(teams)])
n = len(teams)

games['WTeamIdx'] = games.WTeamName.apply(
  lambda t: team_to_idx_dict[t]
)
games['LTeamIdx'] = games.LTeamName.apply(
  lambda t: team_to_idx_dict[t]
)
games = games.reset_index().drop('index', axis=1)
games.tail(6)
WTeamName WScore LTeamName LScore WTeamIdx LTeamIdx
46 High_Point 104 UNC_Asheville 100 2 7
47 Charleston_So 79 Gardner_Webb 72 0 1
48 UNC_Asheville 92 SC_Upstate 85 7 6
49 Winthrop 78 Radford 74 8 5
50 Charleston_So 71 Presbyterian 70 0 4
51 High_Point 83 Longwood 72 2 3

Looks like there are 9 teams that have played 52 games so far. It’s much more complicated to visualize than our little example with three teams, but here’s a try:

OK, to implement Massey’s algorithm, we’ll initiate a matrix X and a vector y, and iterate through the games setting the values in X and y as we go. We use those to setup and solve the Massey equation.

X = np.zeros([m,n])
y = np.zeros([m])
X.transpose().dot(y)

for i,r in games.iterrows():
    X[i, r.WTeamIdx] = 1
    X[i, r.LTeamIdx] = -1
    y[i] = r.WScore - r.LScore

M = X.transpose().dot(X)
M[-1,:] = 1
p =  X.transpose().dot(y)
p[-1] = 0

r = solve(M,p)
r
array([-2.49328205, -2.94713586, 10.92966588,  0.53566291, -3.28521153,
        1.9915711 , -9.85892221,  3.84363398,  1.28401777])

It’s a bit hard to interpret just the vector r, so we’ll match up the values in r with the teams and sort by rating.

teams = pd.DataFrame(teams)
teams.columns = ['TeamName']
teams['Rating'] = teams.TeamName.apply(
  lambda t: r[team_to_idx_dict[t]]
)
teams.sort_values('Rating', ascending=False)
TeamName Rating
2 High_Point 10.929666
5 UNC_Asheville 3.843634
4 Radford 1.991571
0 Winthrop 1.284018
7 Longwood 0.535663
3 Charleston_So -2.493282
6 Gardner_Webb -2.947136
1 Presbyterian -3.285212
8 SC_Upstate -9.858922

Hey, not bad!