[統計] GLMの最尤推定を自分でやってみる

はじめに

ロジスティック回帰

対数尤度

対数尤度の計算

statsmodelsでGLM

実際の最尤推定の手順

メモ

関連記事

はじめに

GLMで対数尤度を自分で計算して自分で最尤推定をやってみる、という回。

ロジスティック回帰

今回はロジスティック回帰を行ってみる。

ガス濃度とカブトムシの死亡数のデータを使う（Dobson, 2008, p.144）。

x	n	y
1.6907	59	6
1.7242	60	13
1.7552	62	18
1.7842	56	28
1.8113	63	52
1.8369	59	53
1.8610	62	61
1.8839	60	60

統計モデルは次のような式で表される。

$y_i \sim Binomial(n_i, \pi_i)$

$\pi_i \sim N(\beta_1 + \beta_2 x_i, \sigma^2)$

対数尤度

2パラメータのロジスティック回帰の場合の対数尤度関数を示す (Dobson, 2008, p.145)。

$l = \sum_{i=1}^N \left[y_i (\beta_1 + \beta_2 x_i) - n_i \log [1 + exp(\beta_1 + \beta_2 x_i) ] + \log\begin{pmatrix} n_i \\ y_i \end{pmatrix} \right]$

対数尤度の計算

パラメータの値を変えて対数尤度を計算してみる。

from itertools import product

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
from scipy.special import comb


def loglikelihood_function(Y, N, X, beta1, beta2):
    res = 0
    for y, n, x in zip(Y, N, X):
        res += y * (beta1 + beta2 * x) - n * np.log(1 + np.exp(beta1 + beta2 * x)) + np.log(comb(n, y))
    return res


def calc_loglikelihood():
    df = pd.read_csv('beetle.csv')
    Y = df['y'].values
    X = np.ones(shape=[len(df), 2]) * np.nan
    X = df['x'].values
    N = df['n'].values
    RESOLUTION = 100
    beta1_arr = np.linspace(-80, -40, RESOLUTION)
    beta2_arr = np.linspace(15, 55, RESOLUTION)
    res = []
    for beta1, beta2 in tqdm(list(product(beta1_arr, beta2_arr))):
        loglikelihood = loglikelihood_function(Y=Y, N=N, X=X, beta1=beta1, beta2=beta2)
        res.append([beta1, beta2, loglikelihood])
    res = pd.DataFrame(res, columns=['beta1', 'beta2', 'loglikelihood'])
    res.to_csv('loglikelihood.csv')
    print(res[res.index == res['loglikelihood'].idxmax()])
    return


def plot_heatmap():
    df = pd.read_csv('loglikelihood.csv')
    df = df.pivot(index='beta2', columns='beta1', values='loglikelihood')
    df.to_csv('heatmap.csv')
    fig, ax = plt.subplots()
    sns.heatmap(data=df, ax=ax, cmap='YlGnBu')
    ax.set_xticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_xticklabels()])
    ax.set_yticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_yticklabels()])
    fig.tight_layout()
    fig.savefig('heatmap_loglikelohood.jpg', dpi=300)
    return


if __name__ == '__main__':
    calc_loglikelihood()
    plot_heatmap()

         beta1      beta2  loglikelihood
4947 -60.20202  33.989899      -18.72801

f:id:cyanatlas:20201018141256j:plain

図では分かりにくいが、対数尤度が最大になるのは $\beta_1=-60, \beta_2=34$ あたり。

statsmodelsでGLM

statsmodelsライブラリでロジスティック回帰を行ってみる。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf


def glm_statsmodels():
    df = pd.read_csv('beetle.csv')
    df['z'] = df['n'] - df['y']
    glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
    fit = glm.fit()
    print(fit.summary())
    return


if __name__ == '__main__':
    glm_statsmodels()

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ['y', 'z']   No. Observations:                    8
Model:                            GLM   Df Residuals:                        6
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -18.715
Date:                Sun, 18 Oct 2020   Deviance:                       11.232
Time:                        13:56:22   Pearson chi2:                     10.0
No. Iterations:                     6
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -60.7175      5.181    -11.720      0.000     -70.871     -50.563
x             34.2703      2.912     11.768      0.000      28.563      39.978
==============================================================================

となって、対数尤度は $\beta_1=-60, \beta_2=34$ あたりで最大となることが分かる。この結果は、先に計算した結果と一致する。