An Improved Analysis of NBA Foul Calls with Python

%matplotlib inline

import datetime
from itertools import product
import logging
import pickle

from matplotlib import pyplot as plt
from matplotlib.offsetbox import AnchoredText
from matplotlib.ticker import FuncFormatter, StrMethodFormatter
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
from sklearn.preprocessing import LabelEncoder
from theano import tensor as tt

pct_formatter = StrMethodFormatter('{x:.1%}')

sns.set()
blue, green, *_ = sns.color_palette()

plt.rc('figure', figsize=(8, 6))

LABELSIZE = 14
plt.rc('axes', labelsize=LABELSIZE)
plt.rc('axes', titlesize=LABELSIZE)
plt.rc('figure', titlesize=LABELSIZE)
plt.rc('legend', fontsize=LABELSIZE)
plt.rc('xtick', labelsize=LABELSIZE)
plt.rc('ytick', labelsize=LABELSIZE)

SEED = 207183 # from random.org, for reproducibility

# keep theano from complaining about compile locks for small models
(logging.getLogger('theano.gof.compilelock')
        .setLevel(logging.CRITICAL))

%%bash
DATA_URI=https://raw.githubusercontent.com/polygraph-cool/last-two-minute-report/32f1c43dfa06c2e7652cc51ea65758007f2a1a01/output/all_games.csv
DATA_DEST=/tmp/all_games.csv

if [[ ! -e $DATA_DEST ]];
then
    wget -q -O $DATA_DEST $DATA_URI
fi

USECOLS = [
    'period',
    'seconds_left',
    'call_type',
    'committing_player',
    'disadvantaged_player',
    'review_decision',
    'play_id',
    'away',
    'home',
    'date',
    'score_away',
    'score_home',
    'disadvantaged_team',
    'committing_team'
]

orig_df = pd.read_csv(
    '/tmp/all_games.csv',
    usecols=USECOLS,
    index_col='play_id',
    parse_dates=['date']
)

orig_df.shape[0]

orig_df.head(n=2).T

(orig_df['call_type']
        .value_counts()
        .head(n=15))

(orig_df['call_type']
        .str.split(':', expand=True)
        .iloc[:, 0]
        .value_counts()
        .plot(
            kind='bar',
            color=blue, logy=True, 
            title="Call types"
        )
        .set_ylabel("Frequency"));

foul_df = orig_df[
    orig_df['call_type']
           .fillna("UNKNOWN")
           .str.startswith("Foul")
]

(foul_df['call_type']
        .str.split(': ', expand=True)
        .iloc[:, 1]
        .value_counts()
        .plot(
            kind='bar',
            color=blue, logy=True,
            title="Foul Types"
        )
        .set_ylabel("Frequency"));

FOULS = [
    f"Foul: {foul_type}"
        for foul_type in [
            "Personal",
            "Shooting",
            "Offensive",
            "Loose Ball",
            "Away from Play"
        ]
]

TEAM_MAP = {
    "NKY": "NYK",
    "COS": "BOS",
    "SAT": "SAS",
    "CHi": "CHI",
    "LA)": "LAC",
    "AT)": "ATL",
    "ARL": "ATL"
}

def correct_team_name(col):
    def _correct_team_name(df):
        return df[col].apply(lambda team_name: TEAM_MAP.get(team_name, team_name))
    
    return _correct_team_name

def date_to_season(date):
    if date >= datetime.datetime(2017, 10, 17):
        return '2017-2018'
    elif date >= datetime.datetime(2016, 10, 25):
        return '2016-2017'
    elif date >= datetime.datetime(2015, 10, 27):
        return '2015-2016'
    else:
        return '2014-2015'

clean_df = (foul_df.where(lambda df: df['period'] == "Q4")
                   .where(lambda df: (df['date'].between(datetime.datetime(2016, 10, 25),
                                                         datetime.datetime(2017, 4, 12))
                                     | df['date'].between(datetime.datetime(2015, 10, 27),
                                                          datetime.datetime(2016, 5, 30)))
                   )
                   .assign(
                       review_decision=lambda df: df['review_decision'].fillna("INC"),
                       committing_team=correct_team_name('committing_team'),
                       disadvantged_team=correct_team_name('disadvantaged_team'),
                       away=correct_team_name('away'),
                       home=correct_team_name('home'),
                       season=lambda df: df['date'].apply(date_to_season)
                   )
                   .where(lambda df: df['call_type'].isin(FOULS))
                   .dropna()
                   .drop('period', axis=1)
                   .assign(call_type=lambda df: (df['call_type']
                                                   .str.split(': ', expand=True)  
                                                   .iloc[:, 1])))

clean_df.shape[0] / orig_df.shape[0]

clean_df.head(n=2).T

call_type_enc = LabelEncoder().fit(
    clean_df['call_type']
)
n_call_type = call_type_enc.classes_.size

player_enc = LabelEncoder().fit(
    np.concatenate((
        clean_df['committing_player'],
        clean_df['disadvantaged_player']
    ))
)
n_player = player_enc.classes_.size

season_enc = LabelEncoder().fit(
    clean_df['season']
)
n_season = season_enc.classes_.size

df = (clean_df[['seconds_left']]
              .round(0)
              .assign(
                call_type=call_type_enc.transform(clean_df['call_type']),
                foul_called=1. * clean_df['review_decision'].isin(['CC', 'INC']),
                player_committing=player_enc.transform(clean_df['committing_player']),
                player_disadvantaged=player_enc.transform(clean_df['disadvantaged_player']),
                score_committing=clean_df['score_home'].where(
                    clean_df['committing_team'] == clean_df['home'],
                    clean_df['score_away']
                ),
                score_disadvantaged=clean_df['score_home'].where(
                    clean_df['disadvantaged_team'] == clean_df['home'],
                    clean_df['score_away']
                ),
                season=season_enc.transform(clean_df['season'])
              ))

df.head(n=2).T

def make_foul_rate_yaxis(ax, label="Observed foul call rate"):
    ax.yaxis.set_major_formatter(pct_formatter)
    ax.set_ylabel(label)
    
    return ax

make_foul_rate_yaxis(
    df.pivot_table('foul_called', 'season')
      .rename(index=season_enc.inverse_transform)
      .rename_axis("Season")
      .plot(kind='bar', rot=0, legend=False)
);

import pymc3 as pm

with pm.Model() as base_model:
    β_season = pm.Normal('β_season', 0., 5., shape=n_season)
    p = pm.Deterministic('p', pm.math.sigmoid(β_season))

season = df['season'].values

with base_model:
    y = pm.Bernoulli(
        'y', p[season],
        observed=df['foul_called']
    )

NJOBS = 3

SAMPLE_KWARGS = {
    'draws': 1000,
    'njobs': NJOBS,
    'random_seed': [
        SEED + i for i in range(NJOBS)
    ],
    'nuts_kwargs': {
        'target_accept': 0.9
    }
}

with base_model:
    base_trace = pm.sample(**SAMPLE_KWARGS)

bfmi = pm.bfmi(base_trace)

max_gr = max(
    np.max(gr_stats) for gr_stats in pm.gelman_rubin(base_trace).values()
)

CONVERGENCE_TITLE = lambda: f"BFMI = {bfmi:.2f}\nGelman-Rubin = {max_gr:.3f}"

(pm.energyplot(base_trace, legend=False, figsize=(6, 4))
   .set_title(CONVERGENCE_TITLE()));

base_trace['p']

resid_df = (df.assign(p_hat=base_trace['p'][:, df['season']].mean(axis=0))
              .assign(resid=lambda df: df['foul_called'] - df['p_hat']))

resid_df[['foul_called', 'p_hat', 'resid']].head()

(resid_df.pivot_table('resid', 'season')
         .rename(index=season_enc.inverse_transform))

def make_time_axes(ax,
                   xlabel="Seconds remaining in game",
                   ylabel="Observed foul call rate"):
    ax.invert_xaxis()
    ax.set_xlabel(xlabel)
    
    return make_foul_rate_yaxis(ax, label=ylabel)

make_time_axes(
    resid_df.pivot_table('resid', 'seconds_left')
            .reset_index()
            .plot('seconds_left', 'resid', kind='scatter'),
    ylabel="Residual"
);

df['trailing_committing'] = (df['score_committing']
                               .lt(df['score_disadvantaged'])
                               .mul(1.)
                               .astype(np.int64))

make_time_axes(
    df.pivot_table(
        'foul_called',
        'seconds_left',
        'trailing_committing'
      )
      .rolling(20)
      .mean()
      .rename(columns={
          0: "No", 1: "Yes"
      })
      .rename_axis(
          "Committing team is trailing",
          axis=1
      )
      .plot()
);

ax = (df.pivot_table('foul_called', 'call_type')
        .rename(index=call_type_enc.inverse_transform)
        .rename_axis("Call type", axis=0)
        .plot(kind='barh', legend=False))

ax.xaxis.set_major_formatter(pct_formatter);
ax.set_xlabel("Observed foul call rate");

with pm.Model() as poss_model:
    β_season = pm.Normal('β_season', 0., 5., shape=2)

def hierarchical_normal(name, shape, σ_shape=1):
    Δ = pm.Normal(f'Δ_{name}', 0., 1., shape=shape)
    σ = pm.HalfNormal(f'σ_{name}', 5., shape=σ_shape)
    
    return pm.Deterministic(name, Δ * σ)

with poss_model:
    β_call = hierarchical_normal('β_call', n_call_type)

df['score_diff'] = (df['score_disadvantaged']
                      .sub(df['score_committing']))

df['trailing_poss'] = (df['score_diff']
                         .div(3)
                         .apply(np.ceil))

trailing_poss_enc = LabelEncoder().fit(df['trailing_poss'])
trailing_poss = trailing_poss_enc.transform(df['trailing_poss'])
n_trailing_poss = trailing_poss_enc.classes_.size

make_time_axes(
    df.pivot_table(
        'foul_called',
        'seconds_left',
        'trailing_poss'
      )
      .loc[:, 1:3]
      .rolling(20).mean()
      .rename_axis(
          "Trailing possessions\n(committing team)",
          axis=1
      )
      .plot()
);

make_time_axes(
    df.pivot_table('foul_called', 'seconds_left', 'call_type')
      .rolling(20).mean()
      .rename(columns=call_type_enc.inverse_transform)
      .rename_axis(None, axis=1)
      .plot()
);

df['remaining_poss'] = (df['seconds_left']
                          .floordiv(25)
                          .add(1))

remaining_poss_enc = LabelEncoder().fit(df['remaining_poss'])
remaining_poss = remaining_poss_enc.transform(df['remaining_poss'])
n_remaining_poss = remaining_poss_enc.classes_.size

ax = sns.heatmap(
    df.pivot_table(
        'foul_called',
        'trailing_poss',
        'remaining_poss'
      )
      .rename_axis(
          "Trailing possessions\n(committing team)",
          axis=0
      )
      .rename_axis("Remaining possessions", axis=1),
    cmap='seismic',
    cbar_kws={'format': pct_formatter}
)

ax.invert_yaxis();
ax.set_title("Observed foul call rate");

call_name_df = df.assign(
    call_type=lambda df: call_type_enc.inverse_transform(
        df['call_type'].values
    )
)

diff_df = (pd.merge(
                call_name_df,
                call_name_df.groupby('call_type')
                            ['foul_called']
                            .mean()
                            .rename('avg_foul_called')
                            .reset_index()
             )
             .assign(diff=lambda df: df['foul_called'] - df['avg_foul_called']))

def plot_foul_diff_heatmap(*_, data=None, **kwargs):
    ax = plt.gca()

    sns.heatmap(
        data.pivot_table(
            'diff',
            'trailing_poss',
            'remaining_poss'
        ),
        cmap='seismic', robust=True,
        cbar_kws={'format': pct_formatter}
    )
    
    ax.invert_yaxis()
    ax.set_title("Observed foul call rate")

(sns.FacetGrid(diff_df, col='call_type', col_wrap=3, aspect=1.5)
    .map_dataframe(plot_foul_diff_heatmap)
    .set_axis_labels(
        "Remaining possessions",
        "Trailing possessions\n(committing team)"
    )
    .set_titles("{col_name}"));

with poss_model:
    β_poss = hierarchical_normal(
        'β_poss',
        (n_trailing_poss, n_remaining_poss, n_call_type),
        σ_shape=(1, 1, n_call_type)
    )

call_type = df['call_type'].values

with poss_model:
    η_game = β_season[season] \
                + β_call[call_type] \
                + β_poss[
                    trailing_poss,
                    remaining_poss,
                    call_type
                ]

with poss_model:
    p = pm.Deterministic('p', pm.math.sigmoid(η_game))
    y = pm.Bernoulli('y', p, observed=df['foul_called'])

with poss_model:
    poss_trace = pm.sample(**SAMPLE_KWARGS)

bfmi = pm.bfmi(poss_trace)

max_gr = max(
    np.max(gr_stats) for gr_stats in pm.gelman_rubin(poss_trace).values()
)

(pm.energyplot(poss_trace, legend=False, figsize=(6, 4))
   .set_title(CONVERGENCE_TITLE()));

resid_df = (df.assign(p_hat=poss_trace['p'].mean(axis=0))
              .assign(resid=lambda df: df.foul_called - df.p_hat))

ax = sns.heatmap(
    resid_df.pivot_table(
                'resid',
                'trailing_poss',
                'remaining_poss'
            )
            .rename_axis(
                "Trailing possessions\n(committing team)",
                axis=0
            )
            .rename_axis(
                "Remaining possessions",
                axis=1
            )
            .loc[-3:3],
    cmap='seismic', 
    cbar_kws={'format': pct_formatter}
)

ax.invert_yaxis();
ax.set_title("Observed foul call rate");

N_BIN = 20

bin_ix, bins = pd.qcut(
    resid_df.p_hat, N_BIN,
    labels=np.arange(N_BIN),
    retbins=True
)

ax = (resid_df.groupby(bins[bin_ix])
              .resid.mean()
              .rename_axis('p_hat', axis=0)
              .reset_index()
              .plot('p_hat', 'resid', kind='scatter'))

ax.xaxis.set_major_formatter(pct_formatter);
ax.set_xlabel(r"Binned $\hat{p}$");

make_foul_rate_yaxis(ax, label="Residual");

ax = (resid_df.groupby('seconds_left')
              .resid.mean()
              .reset_index()
              .plot('seconds_left', 'resid', kind='scatter'))
make_time_axes(ax, ylabel="Residual");

MODEL_NAME_MAP = {
    0: "Base",
    1: "Possession"
}

comp_df = (pm.compare(
                (base_trace, poss_trace),
                (base_model, poss_model)
             )
             .rename(index=MODEL_NAME_MAP)
             .loc[MODEL_NAME_MAP.values()])

comp_df

fig, ax = plt.subplots()

ax.errorbar(
    np.arange(len(MODEL_NAME_MAP)),
    comp_df.WAIC,
    yerr=comp_df.SE, fmt='o'
);

ax.set_xticks(np.arange(len(MODEL_NAME_MAP)));
ax.set_xticklabels(comp_df.index);
ax.set_xlabel("Model");

ax.set_ylabel("WAIC");

with pm.Model() as irt_model:
    β_season = pm.Normal('β_season', 0., 5., shape=n_season)
    β_call = hierarchical_normal('β_call', n_call_type)
    β_poss = hierarchical_normal(
        'β_poss',
        (n_trailing_poss, n_remaining_poss, n_call_type),
        σ_shape=(1, 1, n_call_type)
    )
    
    η_game = β_season[season] \
                + β_call[call_type] \
                + β_poss[
                    trailing_poss,
                    remaining_poss,
                    call_type
                ]

player_disadvantaged = df['player_disadvantaged'].values
n_player = player_enc.classes_.size

with irt_model:
    θ_player = hierarchical_normal(
        'θ_player', (n_player, n_season)
    )
    θ = θ_player[player_disadvantaged, season]

player_committing = df['player_committing'].values

with irt_model:
    b_player = hierarchical_normal(
        'b_player', (n_player, n_season)
    )
    b = b_player[player_committing, season]

with irt_model:
    η_player = θ - b

with irt_model:
    η = η_game + η_player

with irt_model:
    p = pm.Deterministic('p', pm.math.sigmoid(η))
    y = pm.Bernoulli(
        'y', p,
        observed=df['foul_called']
    )

with irt_model:
    irt_trace = pm.sample(**SAMPLE_KWARGS)

bfmi = pm.bfmi(irt_trace)

max_gr = max(
    np.max(gr_stats) for gr_stats in pm.gelman_rubin(irt_trace).values()
)

(pm.energyplot(irt_trace, legend=False, figsize=(6, 4))
   .set_title(CONVERGENCE_TITLE()));

resid_df = (df.assign(p_hat=irt_trace['p'].mean(axis=0))
              .assign(resid=lambda df: df['foul_called'] - df['p_hat']))

N_BIN = 50

bin_ix, bins = pd.qcut(
    resid_df.p_hat, N_BIN,
    labels=np.arange(N_BIN),
    retbins=True
)

ax = (resid_df.groupby(bins[bin_ix])
              .resid.mean()
              .rename_axis('p_hat', axis=0)
              .reset_index()
              .plot('p_hat', 'resid', kind='scatter'))

ax.xaxis.set_major_formatter(pct_formatter);
ax.set_xlabel(r"Binned $\hat{p}$");

make_foul_rate_yaxis(ax, label="Residual");

ax = (resid_df.groupby('seconds_left')
              .resid.mean()
              .reset_index()
              .plot('seconds_left', 'resid', kind='scatter'))
make_time_axes(ax, ylabel="Residual");

MODEL_NAME_MAP[2] = "IRT"

comp_df = (pm.compare(
                (base_trace, poss_trace, irt_trace),
                (base_model, poss_model, irt_model)
             )
             .rename(index=MODEL_NAME_MAP)
             .loc[MODEL_NAME_MAP.values()])

comp_df

fig, ax = plt.subplots()
ax.errorbar(
    np.arange(len(MODEL_NAME_MAP)), comp_df.WAIC,
    yerr=comp_df.SE, fmt='o'
);
ax.set_xticks(np.arange(len(MODEL_NAME_MAP)));
ax.set_xticklabels(comp_df.index);
ax.set_xlabel("Model");
ax.set_ylabel("WAIC");

def varname_to_param(varname):
    return varname[0]

def varname_to_player(varname):
    return int(varname[3:-2])

def varname_to_season(varname):
    return int(varname[-1])

irt_df = (pm.trace_to_dataframe(
                irt_trace, varnames=['θ_player', 'b_player']
            )
            .rename(columns=lambda col: col.replace('_player', ''))
            .T
            .apply(
               lambda s: pd.Series.describe(
                   s, percentiles=[0.055, 0.945]
               ),
               axis=1
            )
            [['mean', '5.5%', '94.5%']]
            .rename(columns={
                 '5.5%': 'low',
                 '94.5%': 'high'
            })
            .rename_axis('varname')
            .reset_index()
            .assign(
                param=lambda df: df['varname'].apply(varname_to_param),
                player=lambda df: df['varname'].apply(varname_to_player),
                season=lambda df: df['varname'].apply(varname_to_season)
            )
            .drop('varname', axis=1))

irt_df.head()

player_irt_df = irt_df.pivot_table(
    index='player',
    columns=['param', 'season'],
    values='mean'
)

player_irt_df.head()

def plot_latent_params(df):
    fig, ax = plt.subplots()
    
    n, _ = df.shape
    y = np.arange(n)

    ax.errorbar(
        df['mean'], y,
        xerr=(df[['high', 'low']]
                .sub(df['mean'], axis=0)
                .abs()
                .values.T),
        fmt='o'
    )

    ax.set_yticks(y)
    ax.set_yticklabels(
        player_enc.inverse_transform(df.player)
    )
    ax.set_ylabel("Player")
    
    return fig, ax

fig, axes = plt.subplots(
    ncols=2, nrows=2, sharex=True,
    figsize=(16, 8)
)
(θ0_ax, θ1_ax), (b0_ax, b1_ax) = axes

bins = np.linspace(
    0.9 * irt_df['mean'].min(),
    1.1 * irt_df['mean'].max(),
    75
)

θ0_ax.hist(
    player_irt_df['θ', 0],
    bins=bins, normed=True
);
θ1_ax.hist(
    player_irt_df['θ', 1],
    bins=bins, normed=True
);

θ0_ax.set_yticks([]);
θ0_ax.set_title(
    r"$\hat{\theta}$ (" + season_enc.inverse_transform(0) + ")"
);

θ1_ax.set_yticks([]);
θ1_ax.set_title(
    r"$\hat{\theta}$ (" + season_enc.inverse_transform(1) + ")"
);

b0_ax.hist(
    player_irt_df['b', 0],
    bins=bins, normed=True, color=green
);
b1_ax.hist(
    player_irt_df['b', 1],
    bins=bins, normed=True, color=green
);

b0_ax.set_xlabel(
    r"$\hat{b}$ (" + season_enc.inverse_transform(0) + ")"
);

b0_ax.invert_yaxis();
b0_ax.xaxis.tick_top();
b0_ax.set_yticks([]);

b1_ax.set_xlabel(
    r"$\hat{b}$ (" + season_enc.inverse_transform(1) + ")"
);

b1_ax.invert_yaxis();
b1_ax.xaxis.tick_top();
b1_ax.set_yticks([]);

fig.suptitle("Disadvantaged skill", size=18);
fig.text(0.45, 0.02, "Committing skill", size=18)
fig.tight_layout();

top_bot_irt_df = (irt_df.groupby('param')
                    .apply(
                        lambda df: pd.concat((
                                df.nlargest(10, 'mean'),
                                df.nsmallest(10, 'mean')
                            ),
                            axis=0, ignore_index=True
                        )
                    )
                    .reset_index(drop=True))

top_bot_irt_df.head()

fig, ax = plot_latent_params(
    top_bot_irt_df[top_bot_irt_df['param'] == 'θ']
                  .sort_values('mean')
)
ax.set_xlabel(r"$\hat{\theta}$");
ax.set_title("Top and bottom ten");

fig, ax = plot_latent_params(
    top_bot_irt_df[top_bot_irt_df['param'] == 'b']
                  .sort_values('mean')
)
ax.set_xlabel(r"$\hat{b}$");
ax.set_title("Top and bottom ten");

def p_val_to_asterisks(p_val):
    if p_val < 0.0001:
        return "****"
    elif p_val < 0.001:
        return "***"
    elif p_val < 0.01:
        return "**"
    elif p_val < 0.05:
        return "*"
    else:
        return ""

def plot_corr(x, y, **kwargs):
    corrcoeff, p_val = sp.stats.pearsonr(x, y)
    asterisks = p_val_to_asterisks(p_val)
    
    artist = AnchoredText(
        f'{corrcoeff:.2f}{asterisks}',
        loc=10, frameon=False,
        prop=dict(size=LABELSIZE)
    )
    plt.gca().add_artist(artist)
    plt.grid(b=False)

PARAM_MAP = {
    'θ': r"$\hat{\theta}$",
    'b': r"$\hat{b}$"
}

def replace_label(label):
    param, season = eval(label)

    return "{param}\n({season})".format(
        param=PARAM_MAP[param],
        season=season_enc.inverse_transform(season)
    )
    
def style_grid(grid):
    for ax in grid.axes.flat:
        ax.grid(False)
        ax.set_xticklabels([]);
        ax.set_yticklabels([]);
        
        if ax.get_xlabel():
            ax.set_xlabel(replace_label(ax.get_xlabel()))

        if ax.get_ylabel():
            ax.set_ylabel(replace_label(ax.get_ylabel()))
            
    return grid

player_all_season = set(df.groupby('player_disadvantaged')
                          .filter(lambda df: df['season'].nunique() == n_season)
                          ['player_committing']) \
                        & set(df.groupby('player_committing')
                                .filter(lambda df: df['season'].nunique() == n_season)
                                ['player_committing'])

style_grid(
    sns.PairGrid(
        player_irt_df.loc[player_all_season],
        size=1.75
       )
       .map_upper(plt.scatter, alpha=0.5)
       .map_diag(plt.hist)
       .map_lower(plot_corr)
);

MIN = 10

player_has_min = set(df.groupby('player_disadvantaged')
                       .filter(lambda df: (df['season']
                                             .value_counts()
                                             .gt(MIN)
                                             .all()))
                       ['player_committing']) \
                    & set(df.groupby('player_committing')
                            .filter(lambda df: (df['season']
                                                  .value_counts()
                                                  .gt(MIN)
                                                  .all()))
                            ['player_committing'])

grid = style_grid(
    sns.PairGrid(
        player_irt_df.loc[player_has_min],
        size=1.75
       )
       .map_upper(plt.scatter, alpha=0.5)
       .map_diag(plt.hist)
       .map_lower(plot_corr)
)