Massey Example

Calculating the Massey Rating¶

(This is an experiment in posting a Python notebook. The notebook doesn't seem to play well with Blogger and (at least for me) the notebook overlaps the sidebar unless I make my browser full-screen. There are also some problems with line breaks. Let me know in the comments how it works for you. And if anyone wants to figure out what I need to do in the HTML to make the notebook behave I'd be grateful!)

In this posting, I'll demonstrate how to calculate the Massey rating as described in (Massey 1997), which you can download from Papers. (His current rating system may be somewhat different.) I was a little surprised to see that I haven't ever written about the Massey ratings, so this will be a chance to do that. This is also an opportunity for me to experiment with creating an iPython Notebook!

You can read the details of the Massey rating system in his paper, but the basic notion is that the expected outcome of a game (in terms of Margin of Victory) is equal to the difference in ratings between the two teams:

$e_{ij} = R_i - R_j$ So how do we figure out the ratings?
After a few games, we have some actual outcomes and we can use these to calculate the team ratings. If we think of every actual game as equation with a value for $e$ (the actual margin by which $T_i$ won the game), then we have a set of linear equations with a variable ($R_i$) for each team. Once we have more games than teams we can solve for $R_i$. In the rest of this notebook I'll walk through how to do this in Python.

The first thing we'll need is some sample data for testing. In this case I'm going to have 10 games for 6 teams. I'll use a simple representation of a game as a list of three numbers:

$[T_i, T_j, Win_i]$ Teams will be represented by id numbers from 0 to 5, and $Win_i$ will be the margin of victory for $T_i$. So if Team 0 played Team 3 and won by 3 points, the game would be represented like this:
$[ 0, 3, 3 ]$ (Obviously real game data is more complex.) So with that game representation, here's a list of games:

In [19]:

games = [[0, 3, 3], [1, 4, 21], [1, 5, 5], [2, 0, 13], [2, 3, 13],
         [3, 0, 25], [4, 5, 8], [4, 3, 15], [5, 1, 21], [5, 4, 8]]

The first thing we want to do is build up a matrix representing the right half of our system of linear equations. (The $R_i - R_j$ part.) The matrix will have one row for every game and one column for every team. To represent the equation for a game, we put a 1 in the column for the winning team (always the first team listed), and -1 in the column for the losing team. (The +1 represents the $R_i$ in the equation above, and the -1 represents the $-R_j$.) The way I've represented games, this is very easy. We just use the team id numbers as the indices for columns and put in 1 or -1.

First, let's initialize an array of zeros of the right size.

In [20]:

import numpy as np
M = np.zeros([10,6])
print M

[[ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]]

That looks about right. Now we'll create a simple function that will take a list of games and a matrix and add the games to the matrix:

In [21]:

def buildGamesMatrix(games,M):
    row = 0
    for g in games:
        M[row, g[0]] = 1
        M[row, g[1]] = -1
        row += 1
    return M

So now let's try this with our sample games.

In [22]:

buildGamesMatrix(games, M)
print M

[[ 1.  0.  0. -1.  0.  0.]
 [ 0.  1.  0.  0. -1.  0.]
 [ 0.  1.  0.  0.  0. -1.]
 [-1.  0.  1.  0.  0.  0.]
 [ 0.  0.  1. -1.  0.  0.]
 [-1.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  1. -1.]
 [ 0.  0.  0. -1.  1.  0.]
 [ 0. -1.  0.  0.  0.  1.]
 [ 0.  0.  0.  0. -1.  1.]]

That looks fine, but it's kind of awkward to have to pass M into the function, so let's modify it to create a new array, populate it, and then return it. The function needs to be told the number of teams, but it can figure out how many games there are:

In [23]:

def buildGamesMatrix(games, num_teams):
    M = np.zeros([len(games), num_teams])
    row = 0
    for g in games:
        M[row, g[0]] = 1
        M[row, g[1]] = -1
        row += 1
    return M

M = buildGamesMatrix(games,6)
print M

[[ 1.  0.  0. -1.  0.  0.]
 [ 0.  1.  0.  0. -1.  0.]
 [ 0.  1.  0.  0.  0. -1.]
 [-1.  0.  1.  0.  0.  0.]
 [ 0.  0.  1. -1.  0.  0.]
 [-1.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  1. -1.]
 [ 0.  0.  0. -1.  1.  0.]
 [ 0. -1.  0.  0.  0.  1.]
 [ 0.  0.  0.  0. -1.  1.]]

Now we need to build another matrix, this one representing the game outcomes. This is just a single column with the $WIN_i$ values for each game. (This is the $e$ side of the linear equations.)

In [24]:

def buildOutcomes(games):
    E = np.zeros([len(games)])
    row = 0
    for g in games:
        E[row] = g[2]
        row += 1
    return E

E = buildOutcomes(games)
print E

[  3.  21.   5.  13.  13.  25.   8.  15.  21.   8.]

(This looks like a row, but believe me it's really a column :-).

Together $M$ and $E$ represent a set of linear equations. Now how do we solve these equations to get the team ratings?

Generally speaking, a set of linear equations like this won't have an exact solution. But we can find the closest approximation through a process called "linear least squares". What this does is find the linear equation that minimizes the total error for all the data points. There are a number of ways to do this in various Python libraries, but one straight forward way is to use the least squares solver in Numpy:

In [25]:

bh = np.linalg.lstsq(M,E)[0]
print bh

[-18.55833333  11.74166667  -1.15833333  -9.75833333   5.24166667
  12.49166667]

These numbers are the team ratings. So $R_0 = -18.55$, $R_1 = 11.74$, etc. Bigger numbers indicate better teams, so you can see that $T_0$ is the worst team and $T_5$ is the best team.

How good are these ratings? We can test this by looking at how well the ratings predicted the actual outcomes. For example, in our first game, $T_0$ played $T_3$ and won by 3 points. According to our model, the expected outcome was $R_0 - R_3$ which is $-18.55 - -9.75$ or about -8.8 points. So in this case, the prediction was off by 11 points (!). We can find the Mean Absolute Error by doing this calculation for every game and then dividing by the number of games:

In [26]:

sum(map(abs, M.dot(bh)-E))/10

Out[26]:

9.3300000000000018

To unpack that cryptic line of Python: The innermost M.dot(bh) puts the calculated ratings bh into the system of linear equations M to calculate the predicted WINs. Then the actual WINs E are subtracted to determine the error. The map takes the absolute values of the errors, and the sum and /10 calculate the mean. The result is a MAE of about 9.3 points.

How can we use the ratings to make predictions? If we want to predict the outcome when $T_1$ plays $T_4$ we go back to our equation:

$e_{14} = R_1 - R_4$ Now we substitute in the values of $R_1=11.74$ and $R_4=5.24$ and predict that $T_1$ will win by 6.5 points.

We can do also do this directly. First we have to create a new matrix representing the game (or games) we want to predict. Suppose that we want to predict a game between $T_0$ and $T_5$. We create a corresponding matrix, and (arbitrarily) mark $T_0$ as the winner and $T_5$ as the loser:

In [27]:

P = np.array([[1, 0, 0, 0, 0, -1]])

And then we dot multiply that with our ratings matrix bh to get a prediction:

In [28]:

print P.dot(bh)

[-31.05]

The prediction is that $T_0$ is going to "win" the game by -31 points. Not a big surprise when the weakest team plays the strongest team!

If we made our prediction matrix the other way:

In [29]:

P = np.array([[-1, 0, 0, 0, 0, 1]])
print P.dot(bh)

[ 31.05]

We get the same result, just expressed as a positive win for $T_5$.

This covers the basic notion of the Massey ratings, but there are a few subtleties that I'll go over in the next installment.