Play with Statistical Modeling: Quantify J-League Team Strength with Stan and Python

Introduction

The other day, I read the book Bayesian Statistical Modeling with Stan and R, which made me want to try statistical modeling. I recommend it because it is a very good book.

Stan can be used not only from R, but also from the Python library PyStan. This time, I will make a model that estimates the strength (like value) of each team in the 2016 soccer J1 league from the match results.

This time's main version

• Python 3.5.1
• pystan==2.14.0.0

Collect data

You can check the J-League match results from http://www.jleague.jp/stats/SFTD01.html. Suppose you have the following data collected from each team in the 1st stage of 2016 in some way. https://gist.github.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b

Thinking

• The number of teams is 18
• In the 1st stage, each team will play only once, so the number of games is ** 153 ** (= 18 * 17/2)
• In soccer, it's almost home and away, and there seems to be a difference between the local and the enemy's land in the match.

Think about the mechanism by which match results are produced

Let's assume how the match result (score) is generated as follows.

• Suppose each team has ** attack power ** ** defense power ** ** home power **.
• The higher the attack power, the higher the score, and the higher the defense power, the lower the goal.
• In home games, both offensive and defensive power are given advantages (= home power)

will do. Well, it doesn't seem that strange story so far.

Next, let's enter the world of statistics and think about it. First, it is appropriately assumed that the "score" follows the Poisson distribution. In other words

Score of a game of a team ~ Poisson (Attack power of that team-Defense power of the opponent team)

It's like that. (By the way, ~ means "follow the distribution".) However, the parameters of the Poisson distribution cannot take negative values, so if you take ʻexp` for convenience and add home power,

Scoring a match for a home team~ Poisson(exp( 
		(Home power+Attack power of the home team) -Away team defense
	))

Suppose that Also symmetrically

Conceded a match in a home team~ Poisson(exp( 
Away team attack power- (Home power+The defense of the home team)
	))

will do.

Write the model in Stan

It turned out to be something like this.

model_code = """
data {
    int N;  // N Games
    int K;  // K Teams
    int Th[N]; // Home Team ID
    int Ta[N]; // Away Team ID
    int Sh[N]; // Home Team score point
    int Sa[N]; // Away Team score point
}

parameters {
    real atk[K];
    real def[K];
    real home_power[K];
    real<lower=0> sigma;
    real<lower=0> hp_sigma;
}

model {
    for (k in 1:K) {
        atk[k] ~ normal(0, sigma);
        def[k] ~ normal(0, sigma);
        home_power[k] ~ normal(0, hp_sigma);
    }

    for (n in 1:N) {
        Sh[n] ~ poisson(exp(
            (home_power[Th[n]] + atk[Th[n]]) - 
            (def[Ta[n]])
        ));

        Sa[n] ~ poisson(exp(
            (atk[Ta[n]]) - 
            (def[Th[n]] + home_power[Th[n]])
        ));
    }
}

generated quantities {
    real games[K, K, 2];
    
    for (th in 1:K) {
        for (ta in 1:K) {
            games[th, ta, 1] = poisson_rng(exp((home_power[th] + atk[th]) - (def[ta])));
            games[th, ta, 2] = poisson_rng(exp((atk[ta]) - (def[th] + home_power[th])));
        }
    }
}
"""

easy explanation

Regarding Stan, the book I introduced earlier is very helpful, but I will explain it briefly here as well.

data block

data {
    int N;  // N Games
    int K;  // K Teams
    int Th[N]; // Home Team ID
    int Ta[N]; // Away Team ID
    int Sh[N]; // Home Team score point
    int Sa[N]; // Away Team score point
}

This data block declares externally given data. The data collected earlier will be injected later with this name.

parameters block

parameters {
    real atk[K];
    real def[K];
    real home_power[K];
    real<lower=0> sigma;
    real<lower=0> hp_sigma;
}

This parameters block describes the parameters you want to estimate. This time, the main focus is attack power (atk), defense power (def), and home power (home_power). <lower = 0> represents the range that the parameter can take.

model block

model {
    for (k in 1:K) {
        atk[k] ~ normal(0, sigma);
        def[k] ~ normal(0, sigma);
        home_power[k] ~ normal(0, hp_sigma);
    }

    for (n in 1:N) {
        Sh[n] ~ poisson(exp(
            (home_power[Th[n]] + atk[Th[n]]) - 
            (def[Ta[n]])
        ));

        Sa[n] ~ poisson(exp(
            (atk[Ta[n]]) - 
            (def[Th[n]] + home_power[Th[n]])
        ));
    }
}

This model block uses a probability distribution to represent the relationship between the input data and the parameters you want to estimate. The "mechanism" mentioned earlier is also expressed here. In addition, it is assumed that the offensive power and defense power of each team follow the "normal distribution of the variance sigma with an average of 0 ". This sigma is also (incidentally) an estimated parameter. Otherwise, the calculation (MCMC) would not converge (because all values would converge in various range values). Well, it seems like I wrote a hope that "I want you to fit in the value around here for the time being". What's interesting is that it was easier for the calculation to converge if it was estimated rather than writing sigma as a fixed 1 (the converged sigma was much smaller than 1). Hard-coding 1 makes me feel that the range is limited, but it is interesting that it is not so.

generated quantities block

generated quantities {
    real games[K, K, 2];
    
    for (th in 1:K) {
        for (ta in 1:K) {
            games[th, ta, 1] = poisson_rng(exp((home_power[th] + atk[th]) - (def[ta])));
            games[th, ta, 2] = poisson_rng(exp((atk[ta]) - (def[th] + home_power[th])));
        }
    }
}

This generated quantities block is like a bonus. It is a function that makes a different calculation using the converged value. Here, the parameters of each team are used to sample the score of the match when each team played against each other. You can calculate it separately on the Python side, but the calculation logic must match the Model, so it seems easier to maintain it if you manage it here.

Other

There is another transformed parameters block, which allows you to define different variables using the data and parameters variables. Expressions such as (home_power [th] + atk [th])-(def [ta]) this time also appear twice in model and generated quantities, so it is better to put them together. I wonder, but ...

Calculate with model and data

Define some useful functions

from hashlib import sha256
import pickle
import time

import requests
import pandas as pd
from pystan import StanModel


def stan_cache(file=None, model_code=None, model_name=None, cache_dir="."):
    """Use just as you would `stan`"""
    # http://pystan.readthedocs.io/en/latest/avoiding_recompilation.html#automatically-reusing-models

    if file:
        if file.startswith("http"):
            model_code = requests.get(file).text
        else:
            with open(file) as f:
                model_code = f.read()
    code_hash = sha256(model_code.encode('utf8')).hexdigest()
    if model_name is None:
        cache_fn = '{}/cached-model-{}.pkl'.format(cache_dir, code_hash)
    else:
        cache_fn = '{}/cached-{}-{}.pkl'.format(cache_dir, model_name, code_hash)
    try:
        sm = pickle.load(open(cache_fn, 'rb'))
    except:
        start_time = time.time()
        sm = StanModel(model_code=model_code)
        print("compile time: %s sec" % (time.time() - start_time))
        with open(cache_fn, 'wb') as f:
            pickle.dump(sm, f)
    return sm


def sampling_summary_to_df(fit4model):
    """

    :param StanFit4model fit4model:
    :return:
    """
    s = fit4model.summary()
    return pd.DataFrame(s["summary"], columns=s['summary_colnames'], index=s['summary_rownames'])

The stan_cache function caches the stan model and retrieves the model code via HTTP as well. Convenient.

The sampling_summary_to_df function turns the result of fitting the model into a Pandas DataFrame. It has the advantage of being easier to see when working with Jupyter.

Execute

import numpy as np
import pandas as pd

df = pd.read_csv("https://gist.githubusercontent.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b/raw/90e1c26e172b92d5f52d53d0591a611e0258bd2b/2016-j1league-1st.csv")

team_list = list(sorted(set(df['away_team']) | set(df['home_team'])))
teams = dict(zip(range(1, len(team_list)+1), team_list))
for k, v in list(teams.items()):
    teams[v] = k
    
model = stan_cache(model_code=model_code) 

stan_input = {
    "N": len(df),
    "K": len(team_list),
    "Th": df['home_team'].apply(lambda x: teams[x]),
    "Ta": df['away_team'].apply(lambda x: teams[x]),
    "Sh": df['home_score'],
    "Sa": df['away_score'],
}

fitted = model.sampling(data=stan_input, seed=999)
sampling_summary_to_df(fitted).sort_values("Rhat", ascending=False)

result

For example, the summary of the parameter estimation results is as follows. In the book, "If this Rhat is less than 1.1 for all estimated parameters, it can be considered that the calculation has converged ", so I will consider it as converged for the time being. However, it is a little subtle that the "number of times the value can be sampled properly" called n_eff is less than 100 (especially when estimating the interval), so there are some subtleties, but this time I will leave it as good.

View the results

Estimated team strength

Collect parameters by team.

tk = {}

params = fitted.extract()
atks = params['atk']
defs = params['def']
home_power = params['home_power']
games = params['games']


for k, team in enumerate(team_list):
    tk[team] = {
            "atk": np.mean(atks[:, k]),
            "def": np.mean(defs[:, k]),
            "home_power": np.mean(home_power[:, k]),
    }

tdf = pd.DataFrame(tk).T

For the time being, I will arrange them in descending order of "home power".

tdf.sort_values("home_power", ascending=False)

have become. This time, the average strength should be 0, so whether it is higher or lower than 0 is more than the average of the J1 team.

By the way, if you do it with the data of the 2nd stage, you will get a completely different result (Kashima will be much worse). To the last, it's just like this when calculated from the data of the 1st stage, so I'm not sure.

If you think about it, it's natural, but it's close to the 2016 1st stage standings. https://ja.wikipedia.org/wiki/2016%E5%B9%B4%E3%81%AEJ1%E3%83%AA%E3%83%BC%E3%82%B0#.E9.A0.86.E4.BD.8D.E8.A1.A8 In particular, there seems to be a general correlation between "total points and total goals" and "attack power and defense power" (although it would be a problem if there were none).

at the end

I tried to simulate what would happen if I bought the 2nd stage Toto with this prediction model, but I stopped because it seems to be a terrible result (laugh). It seems that we can make good predictions if the competition density is high in a short period of time. Predicting 17 verses with data up to 16 verses, is it a little better ??

It's amazing that this kind of modeling has become easy, and it's a lot of fun.