データ解析のための統計モデリング入門 7章 Pythonによるコーディング

本シリーズのモチベーション

データ解析のための統計モデリング入門の内容について理解を深めることを目的に、本書で例示されていた統計モデリング・データ解析に関する処理をPythonでコーディングします（本書ではRでコーディングされていました）。
注：記事は簡潔に書こうと思いますので、本書をご購入・ご一読の上でないと、理解が難しいです。

２章：データ解析のための統計モデリング入門　2章　Pythonによるコーディング
３章：データ解析のための統計モデリング入門　3章　Pythonによるコーディング
６章：データ解析のための統計モデリング入門　6章　Pythonによるコーディング
７章（本ページ）：データ解析のための統計モデリング入門　7章　Pythonによるコーディング

7章を読むモチベーション

7章では、一般化線形混合モデル（Generalized Linear Mixed Model：GLMM)を理解する。GLMMでは、GLMを拡張して、目的変数となるデータのばらつきと、説明変数では表現できないばらつき（ランダム効果）、2種類のばらつきを別々の確率分布で表現する。ランダム効果は、データをサンプリングする対象について、測定できなかった個体差・場所差、もしくは、測定しなかった個体差・場所差を反映した確率的に振舞う量であり、線形予測子の項に含めることで取り扱う。

混合の由来

統計モデルに線形予測子が含まれている場合、その構成要素は固定効果とランダム効果に分類されてきた。ランダム効果は、混合効果とも呼ばれ、そこからランダム効果を線形予測子にもつモデルを混合モデルと呼ぶようになった。

ランダム効果をモデルにどのように取り込むか、具体例を見るとわかりやすい。

あるモデルでは、目的変数$y$の確率分布が2項分布に従い、説明変数が$x$ただ１つと仮定する。2項分布の母数は$0$ ~ $1$の値を取るので、線形予測子と母数を繋ぐリンク関数にはロジット関数を用いる。このモデルをGLMMへと拡張するために、ランダム効果として平均$0$分散$s^2$の正規分布に従う量を線形予測子の項として追加する。ランダム効果は、データごとに独立だと考える。

数式でまとめると、以下のとおりである。

P(y_i, r_i | \beta_0, \beta_1, s) = _NC y_i p_i ^{y_i} (1-p_i) ^{N-y_i}, \quad i = 1, 2, \dots

logit(p_i) = \beta_0 + \beta_1 x_i + r_i, \quad i = 1, 2, \dots

p(r_i | s) = \frac{1}{\sqrt{2\pi s^2}}  \exp (- \frac{r_i^2}{2s^2}) \quad i = 1, 2, \dots

$i$は、データの番号である。データごとに独立した個体差・場所差を考えるために導入した。
$y_i$は、$i$番目のデータの目的変数であり、Nはその上限値である。
$x_i$は、$i$番目のデータの目的変数であり、$\beta_0$（切片）、$\beta_1$（傾き）、$s$（標準偏差）は、GLMMのパラメータである。
$r_i$は、$i$番目のデータのランダム効果であり、$p_i$は$i$番目のデータにおける二項分布の母数の値である。

モデルを構築できれば、あとはパラメータを最尤推定するだけである。最尤推定の方法も、上記の具体例を用いると、理解しやすい。

最大化したい尤度は、$L(\beta_0 , \beta_1 , s) = \Pi _i P(y_i | \beta_0, \beta_1, s)$である。

P(y_i | \beta_0, \beta_1, s) = \int ^{inf} _{-inf} P(y_i , r_i | \beta_0, \beta_1, s) p(r_i | s) dr_i 

$P(y_i | \beta_0, \beta_1, s)$は、同時確率をすべての$r_i$の場合について足すことで求まる。

GLMMをより深く理解するために、次の問題を考える。

ある架空の植物について、同種の植物を複数用意して、各個体（一つの植物）から8つの種子を採取し、そのうちいくつの種子が生存しているかを調べる。この生存種子数がGLMMに従うと仮定して、モデルのパラメータを求める。つまり、線形予測子のパラメータ（説明変数と目的変数の確率分布の母数の関係性）および、ランダム効果を表す分布の母数を求める。モデルの説明変数には、個体の葉の数（$x$）を用いる。モデルにおいて、目的変数の確率分布、線形予測子、リンク関数、そしてランダム効果を表現する確率分布には、上述の具体例として用いたセッティングを使用する。

データの分析

以下からデータを取得し、データフレームに格納する。

# データの概要を表示
df.describe(include='all')

最尤推定

GLMMで最尤推定を行う。

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import binom, norm
from scipy.optimize import minimize
import pandas as pd
from scipy.special import expit
from scipy.integrate import quad


def log_likelihood_binom_glmm(params, x, y, N):
    """
    ロジスティック回帰モデルの対数尤度関数 (GLMM)

    params: [beta0, beta1, sigma] - 線形モデルのパラメータと正規分布の標準偏差
    x: 説明変数
    y: 成功の数
    N: 試行回数
    """
    beta0, beta1, sigma = params

    total_log_likelihood = 0
    for i in range(len(y)):
        def individual_likelihood(random_effect):
            """
            個体iの尤度を計算する関数 (個体差 random_effect を考慮)

            random_effect: 個体iのランダム効果
            """
            logit_p = beta0 + beta1 * x[i] + random_effect
            p = expit(logit_p)  # ロジスティック関数で確率に変換
            return binom.pmf(y[i], N, p) * norm.pdf(random_effect, loc=0, scale=sigma)

        # 個体差を積分消去するために数値積分 (scipy.integrate.quad を使用)
        def integrand(random_effect):
            return individual_likelihood(random_effect)

        random_effect_integral, _ = quad(integrand, -np.inf, np.inf)

        total_log_likelihood += np.log(random_effect_integral)  # 積分後に全体に対数を取る

    return -total_log_likelihood  # 最小化するため、負の対数尤度を返す

# dfからx, y, Nのデータを取り出す
x = df['x'].values
y = df['y'].values
N = df['N'].values[0]

# 最尤推定
initial_params = [0, 0, 1]  # 初期値
result = minimize(log_likelihood_binom_glmm, initial_params, args=(x, y, N), bounds=[(None, None), (None, None), (0, None)])

# 推定されたパラメータ
beta0_hat, beta1_hat, sigma_hat = result.x
print("beta0_hat:", beta0_hat)
print("beta1_hat:", beta1_hat)
print("sigma_hat:", sigma_hat)

max_log_likelihood = -result.fun
print("最大対数尤度:", max_log_likelihood)

beta0_hat: -4.135822320790446
beta1_hat: 0.991820126508299
sigma_hat: 2.497645792843768
最大対数尤度: -198.0739452049279

AICを計算する。
久保先生のWebページを見ると、RのgmmlMLでは、以下のようにAICを定義しているそうだ。ここで、residuak devianceは、フル（飽和）モデルと評価したいモデルの最大対数尤度の差である。

AIC = (residual deviance) + 2×(パラメーター数)

from scipy.stats import binom

def saturated_log_likelihood(y, N):
    """
    飽和モデルの対数尤度を計算
    y: 成功の数
    N: 試行回数
    """
    total_saturated_log_likelihood = 0
    for i in range(len(y)):
        # p_saturated は観測データの割合そのもの
        p_saturated = y[i] / N if N > 0 else 0
        total_saturated_log_likelihood += binom.logpmf(y[i], N, p_saturated)
    
    return total_saturated_log_likelihood

# 飽和モデルの対数尤度を計算
saturated_ll = saturated_log_likelihood(y, N)

# Residual Deviance を計算
residual_deviance = 2 * (saturated_ll - max_log_likelihood)

print("Residual Deviance:", residual_deviance)

# パラメータの数
k = len(result.x)

# AICの計算
AIC = residual_deviance + 2 * k

print("AIC:", AIC)

Residual Deviance: 264.37856703098066
AIC: 270.37856703098066