# 2. Massey’s Method

Massey’s Method refers to Kenneth Massey’s method for ranking sport teams. Denote the following.

• $$n$$ is the number of teams

• $$m$$ is the number games played (only 2 teams can play in a game)

• $$X$$ is a $$n$$ x $$m$$ matrix where columns corresponding to teams and rows correspond to games

• Each row of $$X$$ is sparse (mostly zeros)

• For each row, 1 is placed in the column corresponding to the winning team

• For each row, -1 is placed in the column corresponding to the losing team

• $$r$$ is a $$n$$ x 1 vector corresponding to the rating of each team

• $$r$$ is what will be estimated

• $$y$$ is a $$n$$ x 1 vector corresponding to the margins of victory

Then, we are trying to solve

$$Xr = y$$

Futhermore, denote coefficient matrix $$M$$ as $$M = X^TX$$. Thus, $$M$$ will be $$n$$ x $$n$$ and

• the diagonal elements $$M_{ii}$$ of $$M$$ will be the number of games played by the i-th team,

• the off-diagonal elements $$M_{ij}$$ will be the negation of the number of times the i-th team played the j-th team.

Instead of solving $$Xr = y$$, then we can solve the following,

$$Mr = p$$,

where $$p = X^Ty$$.

## 2.1. Data

Let’s illustrate Massey’s Method using data from the NCAAF ACC conference for the year 2005. Below, t1 and t2 are the 2 teams playing and s1 and s2 are the scores corresponding to t1 and t2. There is no information on who was at home or away.

:

import pandas as pd

f_df

:

t1 t2 s1 s2
0 Duke Miami 7 52
1 Duke UNC 21 24
2 Duke UVA 7 38
3 Duke VT 0 45
4 Miami UNC 34 16
5 Miami UVA 25 17
6 Miami VT 27 7
7 UNC UVA 7 5
8 UNC VT 3 30
9 UVA VT 14 52

We can swap the teams and scores and get a view like the following.

:

r_df = pd.DataFrame([{'t1': r.t2, 't2': r.t1, 's1': r.s2, 's2': r.s1} for _, r in f_df.iterrows()])

game_df = pd.concat([f_df, r_df]).reset_index(drop=True)
game_df['differential'] = game_df.s1 - game_df.s2
game_df.sort_values(['t1', 't2'])

:

t1 t2 s1 s2 differential
0 Duke Miami 7 52 -45
1 Duke UNC 21 24 -3
2 Duke UVA 7 38 -31
3 Duke VT 0 45 -45
10 Miami Duke 52 7 45
4 Miami UNC 34 16 18
5 Miami UVA 25 17 8
6 Miami VT 27 7 20
11 UNC Duke 24 21 3
14 UNC Miami 16 34 -18
7 UNC UVA 7 5 2
8 UNC VT 3 30 -27
12 UVA Duke 38 7 31
15 UVA Miami 17 25 -8
17 UVA UNC 5 7 -2
9 UVA VT 14 52 -38
13 VT Duke 45 0 45
16 VT Miami 7 27 -20
18 VT UNC 30 3 27
19 VT UVA 52 14 38
:

import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('ggplot')

fig, ax = plt.subplots(figsize=(10, 4))

_ = sns.boxplot(x='t1', y='differential', data=game_df, ax=ax)
_ = ax.set_title('Box plots of differentials by team')
_ = ax.set_xlabel('Team')
_ = ax.set_ylabel('Differential')

plt.tight_layout() ## 2.2. Ranking from ratings

A $$n$$ x $$n$$ matrix of point differential for each team can be created as follows. In the matrix below, UNC

• won Duke by 3 points,

• lost to Miami by 18 points,

• win UVA by 2 points, and

• lost to VT by 27 points.

When a team wins, the number in this matrix is positive, and when a team loses, the number in this matrix is negative.

:

import numpy as np

def get_differential(t1, t2):
if t1 == t2:
return np.nan
s = game_df[(game_df.t1 == t1) & (game_df.t2 == t2)].iloc
return s.s1 - s.s2

teams = sorted(list(set(game_df.t1) | set(game_df.t2)))
differentials = [[get_differential(t1, t2) for t2 in teams] for t1 in teams]

diff_df = pd.DataFrame(differentials, index=teams, columns=teams)
diff_df

:

Duke Miami UNC UVA VT
Duke NaN -45.0 -3.0 -31.0 -45.0
Miami 45.0 NaN 18.0 8.0 20.0
UNC 3.0 -18.0 NaN 2.0 -27.0
UVA 31.0 -8.0 -2.0 NaN -38.0
VT 45.0 -20.0 27.0 38.0 NaN

From this point differential matrix, we can then compute the total point differential (sum across the columns) and the number of wins (positive scores) and losses (negative scores).

:

stat_df = pd.DataFrame({
'differential': diff_df.sum(axis=1),
'wins': diff_df.apply(lambda r: len(r[r > 0]), axis=1),
'losses': diff_df.apply(lambda r: len(r[r < 0]), axis=1)
})

stat_df

:

differential wins losses
Duke -124.0 0 4
Miami 91.0 4 0
UNC -40.0 2 2
UVA -17.0 1 3
VT 90.0 3 1

If we let the sum of the point differentials be the ratings, then the ranking produced is as follows where Miami is first and Duke is last.

:

stat_df.differential.sort_values(ascending=False)

:

Miami     91.0
VT        90.0
UVA      -17.0
UNC      -40.0
Duke    -124.0
Name: differential, dtype: float64


We can also let the number of wins be the ratings, and Miami is still first and Duke is last.

:

stat_df.wins.sort_values(ascending=False)

:

Miami    4
VT       3
UNC      2
UVA      1
Duke     0
Name: wins, dtype: int64


Let’s compare the ranking by point differential and wins. Notice that the rankings are nearly identical, except that UVA and UNC switch places?

:

pd.DataFrame({
'by_differential': stat_df.differential.sort_values(ascending=False).index,
'by_wins': stat_df.wins.sort_values(ascending=False).index
})

:

by_differential by_wins
0 Miami Miami
1 VT VT
2 UVA UNC
3 UNC UVA
4 Duke Duke

## 2.3. Xr = y

Let’s build or $$X$$ matrix.

:

def get_vector(r):
v = [0 for i in range(len(teams))]

if r.s1 > r.s2:
v[t2i[r.t1]] = 1
v[t2i[r.t2]] = -1
else:
v[t2i[r.t1]] = -1
v[t2i[r.t2]] = 1
return v

t2i = {t: i for i, t in enumerate(teams)}

X = pd.DataFrame([get_vector(r) for _, r in f_df.iterrows()], columns=teams)
X

:

Duke Miami UNC UVA VT
0 -1 1 0 0 0
1 -1 0 1 0 0
2 -1 0 0 1 0
3 -1 0 0 0 1
4 0 1 -1 0 0
5 0 1 0 -1 0
6 0 1 0 0 -1
7 0 0 1 -1 0
8 0 0 -1 0 1
9 0 0 0 -1 1

Notice that $$X^TX = M$$.

:

X.T.dot(X)

:

Duke Miami UNC UVA VT
Duke 4 -1 -1 -1 -1
Miami -1 4 -1 -1 -1
UNC -1 -1 4 -1 -1
UVA -1 -1 -1 4 -1
VT -1 -1 -1 -1 4

Our $$y$$ will be the point differential.

:

y = f_df.s1 - f_df.s2
y

:

0   -45
1    -3
2   -31
3   -45
4    18
5     8
6    20
7     2
8   -27
9   -38
dtype: int64


We will not attempt to learn the ratings using linear regression.

:

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X, y)
model.intercept_, model.coef_

:

(4.100000000000035, array([ 28.08,  -3.08,   1.6 ,   1.04, -27.64]))


The ratings are the coefficients (see below where we negate them). However, because $$m >> n$$, the linear system is highly overdetermined and inconsistent.

:

sorted(list(zip(teams, -model.coef_)), key=lambda tup: tup, reverse=True)

:

[('VT', 27.640000000000022),
('Miami', 3.080000000000026),
('UVA', -1.0400000000000234),
('UNC', -1.5999999999999974),
('Duke', -28.08000000000003)]


Note how Miami is undefeated, but still in second place? It is likely VT came in first place since their point differentials average higher. In most sports ranking system, point differentials are take out of consideration because teams with a high lead may stop scoring points as a matter of sportsmanship.

## 2.4. Mr = p

Here, we will estimating the ratings and thus the ranking using the coefficient matrix $$M$$.

:

def get_games_played(t1, t2):
if t1 == t2:
return f_df[(f_df.t1 == t1) | (f_df.t2 == t2)].shape
else:
q1 = (f_df.t1 == t1) & (f_df.t2 == t2)
q2 = (f_df.t1 == t2) & (f_df.t2 == t1)
q = q1 | q2
return -f_df[q].shape

mat = [[get_games_played(t1, t2) for t2 in teams] for t1 in teams]
mat

:

[[4, -1, -1, -1, -1],
[-1, 4, -1, -1, -1],
[-1, -1, 4, -1, -1],
[-1, -1, -1, 4, -1],
[-1, -1, -1, -1, 4]]


To ensure that $$M$$ has full rank, we will convert the last row to all 1’s.

:

mat[-1] = [1 for _ in range(len(mat[-1]))]
mat

:

[[4, -1, -1, -1, -1],
[-1, 4, -1, -1, -1],
[-1, -1, 4, -1, -1],
[-1, -1, -1, 4, -1],
[1, 1, 1, 1, 1]]


Our final $$M$$ looks like the following.

:

M = pd.DataFrame(mat, index=teams, columns=teams)
M

:

Duke Miami UNC UVA VT
Duke 4 -1 -1 -1 -1
Miami -1 4 -1 -1 -1
UNC -1 -1 4 -1 -1
UVA -1 -1 -1 4 -1
VT 1 1 1 1 1

Our $$p$$ looks like the following.

:

p = stat_df.differential
p

:

Duke    -124.0
Miami     91.0
UNC      -40.0
UVA      -17.0
VT        90.0
Name: differential, dtype: float64


Applying linear regression and sorting the ratings (coefficients), we get the final ranking below.

:

model = LinearRegression()
model.fit(M, p)
model.intercept_, model.coef_

:

(8.881784197001252e-15, array([-6.8, 36.2, 10. , 14.6, 36. ]))

:

sorted(list(zip(teams, model.coef_)), key=lambda tup: tup, reverse=True)

:

[('Miami', 36.2),
('VT', 36.00000000000002),
('UVA', 14.599999999999998),
('UNC', 9.999999999999998),
('Duke', -6.799999999999999)]


## 2.5. NBA 2021

Let’s apply Massey’s Method to the NBA for the 2021 season and games up to Thanksgiving.

:

def get_nba():
.rename(columns={
'a_team': 't1',
'h_team': 't2',
'a_score': 's1',
'h_score': 's2'})
x = x[x.preseason == False]\
.drop(columns=['preseason'])\
.reset_index(drop=True)
return x

def get_nfl():
.rename(columns={
'team1': 't1',
'team2': 't2',
'score1': 's1',
'score2': 's2'})\
.drop(columns=['week'])
x['t1'] = x['t1'].apply(lambda s: s.strip())
x['t2'] = x['t2'].apply(lambda s: s.strip())

return x

def get_X(df):
def get_vector(r):
v = [0 for i in range(len(teams))]

if r.s1 > r.s2:
v[t2i[r.t1]] = 1
v[t2i[r.t2]] = -1
else:
v[t2i[r.t1]] = -1
v[t2i[r.t2]] = 1
return v

teams = sorted(list(set(df.t1) | set(df.t2)))
t2i = {t: i for i, t in enumerate(teams)}

X = pd.DataFrame([get_vector(r) for _, r in df.iterrows()], columns=teams)
X = X.T.dot(X)
X.iloc[-1,:] = 1

return X

def get_y(df):
def get_diff(t):
a = df[df.t1 == t]
b = df[df.t2 == t]

c = a.s1 - a.s2
d = b.s2 - b.s1

return c.sum() + d.sum()

teams = sorted(list(set(df.t1) | set(df.t2)))
diffs = [get_diff(t) for t in teams]

return pd.Series(diffs, index=teams)

def get_Xy(df):
return get_X(df), get_y(df)

:

X, y = get_Xy(get_nba())

model = LinearRegression()
model.fit(X, y)

pd.DataFrame(zip(X.index, model.coef_), columns=['Team', 'Rating'])\
.sort_values('Rating', ascending=False)

:

Team Rating
28 Warriors 12.899850
10 Jazz 8.549796
24 Suns 6.744632
8 Heat 6.488263
16 Nets 4.254646
2 Bulls 4.235154
5 Clippers 3.792290
7 Hawks 3.338752
18 Pacers 3.138730
27 Trail Blazers 2.456361
0 76ers 2.100543
1 Bucks 1.999356
3 Cavaliers 1.429301
26 Timberwolves 1.276979
9 Hornets 1.227789
29 Wizards 1.180180
21 Raptors 1.055674
17 Nuggets 0.755934
4 Celtics 0.579776
11 Kings 0.089514
12 Knicks -0.306992
15 Mavericks -0.688728
13 Lakers -3.357641
23 Spurs -3.433295
6 Grizzlies -4.525625
19 Pelicans -4.916126
25 Thunder -5.150785
20 Pistons -7.552551
22 Rockets -8.788710
14 Magic -9.818479

## 2.6. NFL

Let’s apply Massey’s Method to the NFL for the 2021 season and games up to Thanksgiving.

:

X, y = get_Xy(get_nfl())

model = LinearRegression()
model.fit(X, y)

pd.DataFrame(zip(X.index, model.coef_), columns=['Team', 'Rating'])\
.sort_values('Rating', ascending=False)

:

Team Rating
3 Bills 8.425259
7 Cardinals 8.235687
21 Patriots 6.349719
6 Buccaneers 6.321022
11 Cowboys 5.819441
10 Colts 3.115332
9 Chiefs 2.706593
29 Titans 2.584503
23 Rams 2.472228
13 Eagles 1.696797
0 49ers 1.565063
19 Packers 1.443518
25 Saints 0.808172
30 Vikings 0.737758
8 Chargers -0.551756
2 Bengals -0.607977
26 Seahawks -0.880815
24 Ravens -1.608072
20 Panthers -1.742287
4 Broncos -2.431469
5 Browns -2.723871
31 Washington -3.649900
27 Steelers -3.707054
22 Raiders -5.414601
15 Giants -6.524430
12 Dolphins -8.305769
1 Bears -9.322540
16 Jaguars -10.338981
28 Texans -11.365255
18 Lions -12.281915
14 Falcons -12.409125
17 Jets -14.139647