「データ解析のための統計モデリング入門」をPythonで写経 Vol.15 ～ ７章「一般化線形混合モデル（GLMM）」②一般化線形混合モデルの実装～パラメータ推定

７章「一般化線形混合モデル（GLMM）」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」７章「一般化線形混合モデル（GLMM）」の Python写経活動記録 です。 

書籍は第７章で 一般化線形混合モデル（GLMM） をやり切ります！
この記事は前回記事の流れを汲み、 一般化線形混合モデルのパラメータ推定を Python で実装 します。

前回記事の「一般化線形混合モデルの骨格」（GLMMのイントロ）も合わせてご覧ください！

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

準備

■ 記事の範囲
この記事はテキスト７章の以下の節を取り扱います。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd
from scipy.special import expit

# 統計計算
import scipy.stats as stats

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

統計モデリング・サマリー

この記事で扱う統計モデリングの概要です。

■ 統計モデル
一般化線形混合モデルです。

■ モデリング手続き

一般化線形混合モデル

データの確認

前回記事と同じ例題データを利用します。
今回記事のデータ確認はざっくり行います。
細かなデータ確認は前回記事をぜひご覧ください。

■ データの読み込み
data.csv ファイルを pandas データフレームの data に読み込みます。

# データの読み込み
data = pd.read_csv('./data/ch07/data.csv')
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 100 です。
植物の個体（個体数 100）に関する仮想の観測データです。
葉数 x の値ごとに 20 の個体を調査しています。

【変数の説明】
植物の個体ごとに採取した調査種子数 N、生存種子数 y（目的変数）、葉数 x です。

■ 生存種子数 y の可視化
y のヒストグラムを描画します。

# yのヒストグラムの描画

# 描画領域の設定
plt.figure(figsize=(5, 4))
# ヒストグラムの描画
sns.histplot(data=data, x='y', bins=np.arange(-0.5, 8.6), edgecolor='white',
             alpha=0.7)
# 修飾
plt.xlabel('生存種子数 $y_i$', fontsize=12)
plt.ylabel('頻度（個体数）', fontsize=12);

【実行結果】
Ｕ字型の分布になっています。

■ データの散布図と真の生存確率に基づく予測値の重ね描き
テキスト p.146 の図 7.2 に相当する散布図を描画します。
「真の生存確率」に基づく生存種子数の予測値的な曲線を重ねます。
真のパラメータ値は $${\beta_1=-4,\ \beta_2=1}$$ です。

# GLMではうまくあつかえない生存種子数の例題 p.146 図7.2(B)

## 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# x軸の値
x_val = np.linspace(data.x.min(), data.x.max(), 100)
# 真の生存確率の一例（点線）のy軸の値（ロジスティック関数）
y_val = expit(beta1 + beta2 * x_val) * data.y.max()

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
# 真の生存確率の一例の点線の描画
ax.plot(x_val, y_val, color='tab:red', ls='--', label='真の生存確率の例')
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存種子数 $y_i$', fontsize=12);
ax.set_title('(B) 全100個体の $x_i$ と $y_i$')
ax.legend(loc='upper left');

【実行結果】

散布図のデータ点の重なりに注目します。
葉数が大きくなるにつれて生存種子数の大きい個体が増える印象です。
とは言うものの、どの葉数でも生存種子数が０～８まで個体はばらついています。

赤点線は、真の生存確率に基づく生存種子数の曲線です。
個人の主観ですが、曲線を見ただけでは、個体のばらつきとの関係を読み取れません（モヤッとします…）。

◆ ◆ ◆

■ データの特徴まとめ
データの特徴を整理します。

データのばらつきは「確率分布」で表現します。
生存種子数の分布の形状がＵ字型をしており、二項分布だけではデータをうまく表現できない感じです。

前回記事で生存種子数に「個体差」があることを確認しました。
二項分布・ロジットリンク関数を適用しつつ、線形予測子に個体差 $${r_i}$$ を追加し、個体差 $${r_i}$$ のばらつきを平均０、標準偏差 $${s}$$ の正規分布と仮定する GLMM（一般化線形混合モデル）でモデリングします。

統計モデルをデータに当てはめ（モデルの理解）

GLMM でモデリングします。

■ 確率分布、リンク関数、線形予測子、ランダム効果
GLMMの４要素である「確率分布」、「リンク関数」、「線形予測子」、「ランダム効果」を定義します。

• 個体差 $${r_i}$$ で条件付けられた生存種子数 $${y_i}$$ は、観測種子数 $${N_i}$$、生存確率 $${q_i}$$ の二項分布に従います。

• リンク関数は「ロジット」、線形予測子は「$${\beta_1 + \beta_2 x_i + r_i}$$」です。

• $${r_i}$$ は平均０、標準偏差 $${s}$$ の正規分布に従います。

■ 無限混合分布
最尤推定では個体差 $${r_i}$$ の推定が難しいため、以下のように積分を行って $${r_i}$$ を消します。
結果は無限混合分布になります。

$${f(r_i \mid 0, s^2)}$$ は個体差が従う正規分布の確率密度関数です。

統計モデルをデータに当てはめ（当てはめと評価）

■ 統計モデルをデータに当てはめ、の準備
GLMMの統計モデルをデータに当てはめします。
統計モデルの対数尤度 $${\log L}$$ が最大になるパラメータ $${\beta_1, \beta_2, s}$$ を推定します。

🔷 今回の統計モデルの尤度関数 $${L}$$ と対数尤度関数 $${\log L}$$

二項分布の確率質量関数 $${p(y_i \mid \beta_1, \beta_2, r_i)}$$、正規分布の確率密度関数 $${f(r_i \mid s)}$$ を用いると、個体 $${i}$$ ごとの尤度 $${L_i}$$ は次のように表せます。

📊 個体ごとの尤度

積分記号 $${\int}$$ の内側で「二項分布の確率 $${\times}$$ 正規分布の確率密度」を計算しています。

全体の尤度関数 $${L}$$ は個体ごとの尤度 $${L_i}$$ の積です。

📊 尤度関数

対数尤度関数は全体の尤度関数 $${L}$$ の対数をとったものになります。

📊 対数尤度関数

◆ ◆ ◆

■ Python で統計モデルをデータに当てはめ
二項分布の GLMM をスマートに実装できる Python ライブラリが見つかりませんでした（泣

そこで、Web記事を参考にさせていただき、スクラッチ実装します。
記事ライターさま、ありがとうございます！

Python 実装では、対数尤度関数を３つのパートに分解して Python 関数化します。

加えて４番目 (4) にて、 Python 関数で負の対数尤度を最小化するパラメータ推定を行います。

数式と Python 関数の結びつきのわかりやすさを重視してコード化します。
処理速度は遅いです。

(2) の個体 $${i}$$ の対数尤度の積分計算には、scipy の 数値積分 quad を利用します。
(4) のパラメータ推定の最適化計算には、scipy の minmize を利用します。

４つの関数を定義します。

# Pythonを使ってGLMMのパラメータを推定 p.159~160
# コードの引用元: https://qiita.com/warper/items/9b53db1b35f3b6c0bbef

# == 追加インポート
from scipy import integrate
from scipy.optimize import minimize

# == 関数定義

# (1)個体差r_iのある値の二項分布pmf×正規分布pdfを算出する関数
def binom_prob_times_norm_prob(r_i, params, x_i, y_i, N_i):
    # パラメータβ1, β2, σをセット
    beta1, beta2, sigma = params
    # 生存確率q_iの算出：q_i = logistic(β1 + β2 * x_i + r_i)
    q_i = expit(beta1 + beta2 * x_i + r_i)
    # 個体差r_iのある値の二項確率×正規確率密度の算出　※p.156の数式の右辺の被積分関数
    prob_r_i = (stats.binom.pmf(k=y_i, n=N_i, p=q_i)
                * stats.norm.pdf(x=r_i, loc=0, scale=sigma))
    # 戻り値：個体差r_iのある値の二項確率×正規確率密度
    return prob_r_i

# (2)個体iの対数尤度を算出する関数
def log_likelihood_for_each_individual(params, x_i, y_i, N_i):
    # 個体iの尤度L_iの算出：個体差r_iを積分消去するための数値積分を実行
    # 数値積分に scipy.integrate.quad を使用 ※p.156のL_iの積分
    L_i, _ = integrate.quad(
        func=binom_prob_times_norm_prob,   # 被積分関数
        a=-np.inf,                         # 積分区間の下端
        b=np.inf,                          # 積分区間の上端
        args=(params, x_i, y_i, N_i),      # 被積分関数のパラメータ
    )
    # 戻り値：個体iの対数尤度
    return np.log(L_i)

# (3)（全体の）対数尤度を算出する関数の定義: 最小化の目的関数
def log_likelihood_binom_glmm(params, x, y, N):
    # 個体ごとの対数尤度L_iを加算して対数尤度を算出
    log_L = sum(
        [log_likelihood_for_each_individual(params, x[i], y[i], N[i])
         for i in range(len(y))]
    )
    # 戻り値：対数尤度（最小化するため、負の対数尤度を返す）
    return -log_L

# (4)最尤推定によるパラメータ推定の実行
def maximum_likelihood_estimation(x, y, N, x0):

    # 初期値設定
    step_count = 0            # 実行ステップカウンター
    start_time = time.time()  # スタート時間

    # callback関数の定義：最適化処理の途中経過の表示
    def callback(xk):
        nonlocal step_count
        step_count += 1
        print(f'{step_count:2d}: Current sol {xk}, '          # Step番号と現在の解
              f'elapsed {time.time() - start_time:6.2f} sec.' # 経過時間（秒）
        )
    
    # 最適化（最小化）の実行
    result = minimize(
        fun=log_likelihood_binom_glmm,                      # 最小化する目的関数
        x0=x0,                                              # パラメータの初期値
        args=(x, y, N),                                     # 目的関数の追加引数
        method='L-BFGS-B',                                  # 最適化方法
        bounds=[(None, None), (None, None), (1e-8, None)],  # パラメータの推定範囲
        callback=callback,                                  # 途中経過を表示
        options={'maxiter': 1000, 'ftol': 1e-9},            # L-BFGS-Bの引数設定
    )

    # 戻り値：最適化の結果（minimizeの戻り値）
    return result

【実行結果】なし

ではパラメータの最尤推定を実行します。

%%time
# 最尤推定によるパラメータ推定の実行

# パラメータ[β1, β2, sigma]の初期値の設定
initial_params = [0, 0, 1]

# データセットの設定
x, y, N = data[['x', 'y', 'N']].values.T

# 最尤推定の実行
result = maximum_likelihood_estimation(x, y, N, initial_params)

【実行結果】
処理時間は 3 分 30 秒ほど掛かりました。
最適化計算のステップごとに推定値の途中経過を表示しています。

パラメータ推定結果を表示します。

# 結果の表示

## パラメータ推定値・最大対数尤度の結果表示
# パラメータ推定値
beta1_hat, beta2_hat, sigma_hat = result.x
print('beta0_hat  :', beta1_hat)  # coef intercept  ※真値は-4
print('beta1_hat  :', beta2_hat)  # coef x          ※真値は 1
print('sigma_hat  :', sigma_hat)  # Scale parameter ※真値は 3
# 最大対数尤度
max_llf = -result.fun # 目的関数で算出した負の対数尤度の符号を反転
print('最大対数尤度:', max_llf)

## 残差逸脱度とAICの算出・表示
# フルモデルの最大対数尤度 ※生存確率pは観測データの生存種子割合
max_llf_full = sum(stats.binom.logpmf(k=data.y, n=data.N, p=data.y / data.N))
# 残差逸脱度
resid_deviance = 2 * (max_llf_full - max_llf)
print('残差逸脱度  :', resid_deviance)
# AIC
k = len(initial_params)       # k=3
aic = resid_deviance + 2 * k  # ※p.76のAICの数式
print('AIC        :', aic)

【実行結果】

テキストの R 言語・glmmML 関数によるパラメータ推定値とほぼ同じになりました。
【テキストのパラメータ推定値など】

なお、パラメータ推定値の標準誤差の算出はうまくいきませんでした（泣

生存種子数の予測

テキストにならって、生存種子数 y の予測を行って可視化します。

■ 予測値を可視化
p.161 図 7.10 に相当します。

# GLMM化したロジスティック回帰の推定にもとづく予測 p.161 図7.10

#### 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# 調査種子数
N = 8

#### 描画処理
## 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

### 左の(A)の描画
## 設定
x_val_a = np.linspace(1, 7, 100)              # x軸の値
y_val_a = expit(beta1 + beta2 * x_val_a) * N  # 真の生存確率の場合の生存種子数y

## 描画
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax1)
# 真の生存確率の赤い点線の描画
ax1.plot(x_val_a, y_val_a, color='tab:red', ls='--', label='真の生存確率')
# GLMMの予測値の緑実線の描画 ※r=0である
ax1.plot(x_val_a, expit(beta1_hat + beta2_hat * x_val_a + 0) * N,
         color='green', lw=2)
# 修飾
ax1.set_xlabel('葉数 $x_i$', fontsize=12)
ax1.set_ylabel('生存種子数 $y_i$', fontsize=12)
ax1.set(xlim=(1.6, 6.4), title='(A) 葉数と生存種子数の関係')

### 右の(B)の描画
## 設定
xi = 4                                           # xiの値=4
num_xi = len(data[data.x==xi])                   # x=4のデータ個数
xi_data = data[data.x==xi]['y'].value_counts().to_frame() # 生存種子数yごとの個体数
x_val_b = np.repeat(xi, N + 1)                   # xi=4をリピートしたデータ
y_val_b = list(range(N + 1))                     # 二項分布の横軸のyの値

## 計算
# 正規分布に従うrの200個のパーセンタイル点の取得
rs2 = stats.norm.ppf(q=np.linspace(0.001, 0.999, 200), loc=0, scale=sigma_hat)
# 生存確率logistic(β1+β2*x+r)の二項分布に従う生存種子数yの取得 shape=(200, 9)
y_mix = np.array(
    [stats.binom.pmf(k=y_val_b, n=N, p=expit(beta1_hat + beta2_hat * x_val_b + r))
     for r in rs2]
) * num_xi

## 描画
# x=4のときの生存種子数yの個体数の描画（青い点）
sns.scatterplot(data=xi_data, x='y', y='count', s=60, alpha=0.5, ax=ax2)
# x=4のときの種子数分布の描画（緑実線） ※y_minの平均値
ax2.plot(y_val_b, y_mix.mean(axis=0), '-o', ms=7, color='green', mec='white')
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax2.set_xlabel('生存種子数 $y_i$', fontsize=12)
ax2.set_ylabel('個体数', fontsize=12)
ax2.set_title(f'(B) 葉数 $x_i=${xi} での種子数分布')
ax2.set_ylim(0, 6.2);

【実行結果】

【左のチャート】
個体差 $${r_i}$$ の影響を調整済みのパラメータ $${\beta_1, \beta_2}$$ を用いて、平均的な個体差 $${r_i = 0}$$ による生存種子数の平均予測値を緑実線で描画しています。
パラメータ真値による赤点線と近い結果を得られていることが分かります。

【右のチャート】
葉数 $${x_i=4}$$ のときの生存種子数の分布を緑実線で描画しています。
コード中の変数 $${\mathtt{y\_min}}$$ で 個体差 $${r_i}$$ を加味した 200 点を計算し、チャート描画時に $${\mathtt{y\_min.mean()}}$$ で平均値を計算して描画しています。
観測値の分布（青い点）に近い結果を得られていることが分かります。

◆ ◆ ◆

GLMM モデルを構築する動機のひとつは、個体差・グループ差を排除した純粋な平均値を求めたい、という点にあるようです。
しかし、GLMM の特徴である「個体差・グループ差」もせっかくなので確認したいです。
そこでパパパっと可視化しよう、とコードを書き始めた途端、厳しい現実を突きつけられました…

積分で $${r_i}$$ を消去したから、個体差が推定されていない！

こんなときは、生成AI におねだりです。

◆ ◆ ◆

（注意）
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。

■ 個体差 $${r_i}$$ の可視化

推定したパラメータ $${\hat{\beta}_1, \hat{\beta}_2, \hat{s}}$$ と観測値を用いて、個体差 $${r_i}$$ を推定します。
ざっくり、個体差 $${r_i}$$ の事後分布の最頻値（MAP推定値）をニュートン法で算出します（と ChatGPT が申しております）。

こちらは個体差 $${r_i}$$ のMAP推定関数です。

# -------------------------------------------------------------------
# ランダム切片 r_i の MAP 推定（= EBLUP）を『やさしく』計算する関数
#    - 1変数のニュートン法（勾配/曲率を解析式で）
#    - r が「データの声」と「0に戻す圧力（事前）」の綱引きのバランス点
# -------------------------------------------------------------------
def map_r_i(beta1, beta2, sigma, x_i, y_i, N_i, tol=1e-8, maxiter=50,
            step_shrink=1.0):
    '''
    個体 i の ランダム切片 r_i をMAP（事後分布の最頻値）推定する
      ℓ(r) = log Binom(y_i | N_i, logistic(β1 + β2 x_i + r)) - r^2/(2σ^2) + 定数
    方法：
      ニュートン法（1変数の二次近似）で解を探す。
      1階微分と2階微分を計算して r を更新。

    引数：
        beta1, beta2, sigma : パラメータ
        x_i, y_i, N_i : 観測データ
        tol : 精度（デフォルト 1e-8）
        maxiter : ニュートン法の反復繰り返し数の最大値（デフォルト 50）
        step_shrink : 1.0 以下にすると一歩を小さくして安定化(極端なy=0/Nの時に有効)
    '''
    
    # 初期値：0（事前の平均）で十分安定
    r = 0.0
    # ニュートン法の実行
    for _ in range(maxiter):
        # 線形予測子ηと生存確率pの算出
        eta = beta1 + beta2 * x_i + r
        p   = expit(eta)
        # 対数事後分布の1階微分・2階微分の算出
        g = (y_i - N_i * p) - r / (sigma**2)           # 1階微分（勾配）
        H = -(N_i * p * (1.0 - p)) - 1.0 / (sigma**2)  # 2階微分（曲率, 常に負）
        # rの更新
        step = step_shrink * (g / H)  # 更新ステップをstep_shrinkで調整
        r_new = r - step              # rの更新
        if abs(r_new - r) < tol:      # rの更新差分が精度以下になったら終了
            return float(r_new)
        r = r_new
    # 戻り値：r ※収束しなくても最後の値を返す
    return float(r)

def map_r_all(beta1, beta2, sigma, x, y, N, step_shrink=1.0):
    '''
    全個体ぶんの MAP 推定値（r_hat）をまとめて返すヘルパー関数
    '''
    return np.array(
        [map_r_i(beta1, beta2, sigma, xi, yi, Ni, step_shrink=step_shrink)
         for xi, yi, Ni in zip(x, y, N)]
    )

【実行結果】なし

関数を使って $${r_i}$$ を推定します。
処理時間はとても短いのでご安心ください。

# ランダム効果r_i推定の実行
est_r_hat = map_r_all(beta1_hat, beta2_hat, sigma_hat, data.x, data.y, data.N)
# 結果の表示（先頭10個）
est_r_hat[:10]

【実行結果】
全個体 100 のうちの先頭 10 個を表示しました。

100 個体の生存種子数の予測値を可視化します。
先ほどの予測値チャートの (A) に重ね描きします。

# GLMM化したロジスティック回帰の推定にもとづく予測 p.161 図7.10(A) の改造

### 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# 調査種子数
N = 8

### 描画処理
## 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

## 設定
x_val_a = np.linspace(1, 7, 100)              # x軸の値
y_val_a = expit(beta1 + beta2 * x_val_a) * N  # 真の生存確率の場合の生存種子数y

## 描画
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15, seed=42,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
# 真の生存確率の赤い点線の描画
ax.plot(x_val_a, y_val_a, color='tab:red', lw=2, ls='--', label='真の生存確率の例')
# GLMMの予測値の緑実線の描画 ※r=0である
ax.plot(x_val_a, expit(beta1_hat + beta2_hat * x_val_a + 0) * N,
        color='green', lw=2, label='予測値(個体差除く)')
# GLMMの予測値のオレンジ実線の描画
for r in est_r_hat:
   ax.plot(x_val_a, expit(beta1_hat + beta2_hat * x_val_a + r) * N,
           color='tab:orange', lw=0.5, alpha=0.3)
ax.plot(0, 0, color='tab:orange', lw=1, alpha=0.5, label='予測値(個体差含む)')
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存種子数 $y_i$', fontsize=12)
ax.set(xlim=(1.6, 6.4), title='(A) 葉数と生存種子数の関係')
ax.legend(bbox_to_anchor=(1, 1), loc='upper left');

【実行結果】
薄オレンジ線が個体ごとの予測値です。

【感想】
チャートを薄オレンジ線が埋め尽くしています。
下側の線は生存種子数の小さな値に、上側の線は生存種子数の大きな値に対応しているように見えます。
個体差の存在を実感できました。
また、中央部に描かれた「個体差 $${r_i = 0}$$ の緑実線」が個体平均の位置づけであることも再認識できました。

◆ ◆ ◆

■ 個体差 $${r_i}$$ の推定コードの解説 by ChatGPT

ChatGPTの解説をご堪能ください！

コードにつながる形で“やさしく＋ていねい”に整理します。
流れは「全体像 → 数式 → 数式の意味 → 数式⇔コードの対応」です。

🔷 全体像（何をしている？）

• 目的：
  個体 i のランダム切片 $${r_i}$$ の MAP（事後分布の最頻値） を求める

• 着想：
  その個体の対数事後分布

　を $${r}$$ について最大化する

• 方法：
  1変数なので ニュートン法 を用います（二次近似でサクッと解く）

🍀🍀🍀

🔷 数式（対数事後・勾配・曲率）

1) モデルの中身

• 線形予測子：
  $${\eta_i=\beta_1+\beta_2 x_i+r_i}$$

• ロジットリンク：
  $${p_i=\sigma(\eta_i)=\dfrac{1}{1+e^{-\eta_i}}}$$

• 尤度（二項分布）：
  $${\log p(y_i\mid r_i,\beta)=y_i\log p_i+(N_i-y_i)\log(1-p_i)}$$

• 事前分布（正規分布）：
  $${\log p(r_i\mid\sigma)=-\dfrac{r_i^2}{2\sigma^2}+C}$$

よって対数事後分布（定数 $${C}$$ は無視）は

2) 一階微分（勾配）
導出のポイント：$${\frac{dp}{d\eta}=p(1-p),\ \frac{d\eta}{dr}=1}$$。

3) 二階微分（曲率；ヘッセ）
常に負なので上に「凸＝一意の最大点」がありニュートン法が安定。

4) ニュートン更新式（最大化）

勾配＝0 を解く近似：

🍀🍀🍀

🔷 数式の“意味合い”（直感）

• 勾配 $${g=(y-Np)-r/\sigma^2}$$ は
  「観測 − 期待（=Np）」 と 「事前の引き戻し（=r/σ²）」 のバランス。

  • もし $${y > Np}$$（観測が想定より多い）
    → $${g>0}$$ → $${r}$$ を 上げる 方向へ。

  • もし $${y < Np}$$ → $${g<0}$$
    → $${r}$$ を 下げる 方向へ。

  • ただし大きすぎる $${r}$$ には 事前罰則 $${r/\sigma^2}$$ がかかり、0 に縮む（シュリンク）。

• 曲率 $${H=-N p(1-p)-1/\sigma^2<0}$$ は「上に凸の強さ」。
  データ量（$${N}$$）が多いほど $${|H|}$$ が大きく、解はキリッと安定。
  $${y=0}$$ や $${y=N}$$ の極端ケースでは $${p}$$ が 0/1 に寄りやすく、$${|H|}$$ が小さくなりがち
  →ステップを縮めると安全。

🍀🍀🍀

🔷 数式 ⇔ コードの対応表

関数 $${\mathtt{map\_r\_i(...)}}$$ と一対一で対応づけます。

eta = beta1 + beta2 * x_i + r
p   = expit(eta)

• 数式の $${\eta=\beta_1+\beta_2 x_i + r}$$、$${p=\sigma(\eta)}$$ に対応。

g = (y_i - N_i * p) - r / (sigma**2)
H = -(N_i * p * (1.0 - p)) - 1.0 / (sigma**2)

• そのまま 勾配 $${g=(y-Np)-r/\sigma^2}$$ と 曲率 $${H=-N p(1-p)-1/\sigma^2}$$。

step = step_shrink * (g / H)
r_new = r - step

• ニュートン更新 $${r_{\text{new}}=r - g/H}$$ に一致。

• $${\mathtt{step\_shrink \leq 1.0}}$$ は ステップ幅 を縮めるダンピング（極端データでの暴走防止）。

if abs(r_new - r) < tol: return float(r_new)

• 収束判定（$${|\Delta r|<\text{tol}}$$ なら終了）。

r = 0.0

• 初期値。事前平均 0 に置くのが素直で十分安定。

🍀🍀🍀

🔷 ちょい上級メモ（現場で効くコツ）

• 極端比率（$${y=0}$$ や $${y=N}$$）が多いときは $${\mathtt{step\_shrink=0.5}}$$ などにすると安定します。

• $${\sigma}$$ は必ず 正。推定の上流で $${\sigma>0}$$ を確保（下限 $${\mathtt{1e-8}}$$ など）。

• 1個体に複数観測 $${j=1, \ldots, m_i}$$ がある場合は

  1. $${g=\sum_j(y_{ij}-N_{ij}p_{ij})-\frac{r}{\sigma^2},\quad H=-\sum_j N_{ij}p_{ij}(1-p_{ij})-\frac{1}{\sigma^2}}$$

  2. と 合算して同じコード（ループ内で和を取る）でOK。

数式だらけでしたね (^_^;;;;)

アディショナルタイム

（注意）
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。

処理時間が短い GLMM コードを ChatGPT に所望したところ…
チャッピー（一般呼称）が爆速コードを書いてくれました！！！

Gauss–Hermite 法で積分を精度良く近似し、かつ、ベクトル化した演算を用いて、高速化を実現しているそうです。

🍀🍀🍀

🔷 GLMM 関数の定義
今回のモデルに「特化した」GLMM 関数を定義します。
汎用的に GLMM をモデリングできるわけではないので、ご了承ください。
コードの概説は後ほど。

# ==============================================
# GLMM（二項×正規ランダム切片, ロジットリンク）
#  ベクトル化 Gauss–Hermite で高速 MLE
#  ・完全版（データ生成→推定→SE→r_hat(MAP)）
#  ・やさしい日本語コメントつき
# ==============================================

# import numpy as np
# from scipy import optimize
# from scipy.special import expit

from scipy.special import logsumexp, gammaln
from statsmodels.tools.numdiff import approx_hess

# --------------------------------------------------------------
# 0) Gauss–Hermite のノード/重み（標準形 e^{-x^2}）を一度だけ作る
#    - K を変えない限り nodes, weights は使い回せます
# --------------------------------------------------------------
def gh_nodes_weights_standard(K: int):
    """
    K : Gauss–Hermite の積分点（大きいほど精度↑, 計算コストも↑）
    戻り値：
      nodes, weights ・・・ e^{-x^2} 用の標準Gauss–Hermiteノード/重み
    """
    nodes, weights = np.polynomial.hermite.hermgauss(K)
    return nodes, weights


# --------------------------------------------------------------------
# 1) ベクトル化した『負の対数尤度』関数（GH 近似）
#    - x,y,N は長さ n の配列（各 i が1グループ）
#    - nodes, weights は標準 GH（0で作って後で σ に合わせてスケール）
# --------------------------------------------------------------------
def nll_glmm_binom_logit_GH_vec(params, x, y, N, *, K, nodes, weights):
    """
    params = [beta1, beta2, log_sigma]
      beta1, beta2 : 固定効果（切片・傾き）
      log_sigma    : ランダム切片の σ を対数で推定（σ>0の制約のため）

    x, y, N        : 各グループ i の説明変数・成功回数・試行回数
    K              : GH の積分点数
    nodes, weights : 標準GHのノード・重み（事前に一度だけ作成）

    戻り値：負の対数尤度（最小化対象）
    """
    beta1, beta2, log_sigma = params
    sigma = np.exp(log_sigma)                       # σ>0 制約を満たす

    # --- GHノードと重みを N(0, σ^2) 用にスケール ---
    #     r_j = √2 * σ * nodes,   w_j = weights / √π
    rj = np.sqrt(2.0) * sigma * nodes              # 形=(K,)
    wj = weights / np.sqrt(np.pi)                  # 形=(K,) かつ合計 ≈ 1
    log_wj = np.log(wj + 1e-300)                   # log(0) 防止の微小値

    # --- 線形予測子を n×K 行列で一気に作る（完全ベクトル化）---
    #     eta[i, j] = (beta1 + beta2 * x[i]) + rj[j]
    eta0 = beta1 + beta2 * x                       # 形=(n,)
    eta  = eta0[:, None] + rj[None, :]             # 形=(n, K)

    # --- 成功確率 p = logistic(eta) を一括計算 ---
    p = expit(eta)                                  # 形=(n, K)

    # --- 二項の log pmf を自前実装（ベクトル化・高速・安定）---
    #     log C(N,y) + y log p + (N-y) log(1-p)
    Ncol = N[:, None].astype(float)                 # 形=(n, 1)
    ycol = y[:, None].astype(float)                 # 形=(n, 1)
    logC = gammaln(Ncol + 1.0) - gammaln(ycol + 1.0) - gammaln(Ncol - ycol + 1.0)
    eps = 1e-300
    logpmf = (                                      # (n,K)
        logC + ycol * np.log(p + eps) + (Ncol - ycol) * np.log(1.0 - p + eps)
      )

    # --- 各個体 i で log Σ_j w_j * pmf(y_i | r_j) を log-sum-exp で安定計算 ---
    log_terms = log_wj[None, :] + logpmf            # 形=(n, K)
    log_Li = logsumexp(log_terms, axis=1)           # 形=(n,)

    # --- 全体の負の対数尤度 ---
    return -np.sum(log_Li)


# ----------------------------------------------------------------
# 2) MLE のメイン関数（完全ベクトル化版）
#    - L-BFGS-B で推定
#    - 近似標準誤差（ヘッセ逆）も返す（σはデルタ法でSEに変換）
# ----------------------------------------------------------------
def fit_glmm_binom_logit_MLE_vec(
    x, y, N, *, K=80, init=(0.0, 0.0, np.log(1.0)), return_hess=False
):
    """
    x, y, N : 1行=1グループの配列（長さ n）
    K       : GH の積分点数（σが大きいほど増やすのがコツ：20〜120目安）
    init    : 初期値 [beta1, beta2, log_sigma]
    return_hess : True でヘッセと共分散行列も返す

    戻り値：
      results : 推定値・SE・最適化情報などの辞書
      res     : scipy.optimize の結果オブジェクト
      （return_hess=True のときは 'H', 'cov' も含める）
    """
    # --- 標準 GH ノード・重みは一度だけ作って使い回し ---
    nodes, weights = gh_nodes_weights_standard(K)

    # --- 最小化を実行（L-BFGS-B）---
    obj = lambda th: nll_glmm_binom_logit_GH_vec(
        th, x, y, N, K=K, nodes=nodes, weights=weights)
    res = optimize.minimize(
        fun=obj, x0=np.array(init, dtype=float), method='L-BFGS-B')

    # --- 推定値（σは log から元スケールへ）
    beta1_mle, beta2_mle, log_sigma_mle = res.x
    sigma_mle = float(np.exp(log_sigma_mle))

    # --- 近似標準誤差（観測情報 ≈ ヘッセ）の逆行列で共分散 ---
    H = cov = None
    try:
        H = approx_hess(res.x, lambda th: nll_glmm_binom_logit_GH_vec(
            th, x, y, N, K=K, nodes=nodes, weights=weights))
        cov = np.linalg.inv(H)
        se_beta1, se_beta2, se_logsig = np.sqrt(np.diag(cov))
        # デルタ法：Var(σ) ≈ (dσ/d logσ)^2 Var(logσ) = σ^2 Var(logσ)
        se_sigma = se_logsig * sigma_mle
    except Exception:
        se_beta1 = se_beta2 = se_sigma = np.nan

    results = {
        'success':  bool(res.success),
        'message':  res.message,
        'n_iter':   int(res.nit),
        'beta1':    float(beta1_mle),
        'beta2':    float(beta2_mle),
        'sigma':    float(sigma_mle),
        'se_beta1': float(se_beta1),
        'se_beta2': float(se_beta2),
        'se_sigma': float(se_sigma),
        'K':        int(K),
    }
    if return_hess:
        results['H'] = H
        results['cov'] = cov
    return results, res

【実行結果】なし

🍀🍀🍀

🔷 GLMM のパラメータ推定の実行
ではでは GLMM のパラメータを最尤推定します！
パラメータ推定値の標準誤差＆個体差も計算してくれます！！

%%time
# Gauss–Hermite法によるGLMMのパラメータ推定

# データセットの設定
x, y, N = data[['x', 'y', 'N']].values.T

# GLMMのパラメータ推定の実行
out, res = fit_glmm_binom_logit_MLE_vec(
    x, y, N, K=100, init=(0.0, 0.0, np.log(1.0)))
print('success:', out['success'], '|', out['message'])
print(f"beta1 MLE={out['beta1']:.4f} (se={out['se_beta1']:.4f})")
print(f"beta2 MLE={out['beta2']:.4f} (se={out['se_beta2']:.4f})")
print(f"sigma MLE={out['sigma']:.4f} (se={out['se_sigma']:.4f}),  K={out['K']}")

# 各個体の r_hat（MAP = EBLUP）の算出
r_hat = map_r_all(out['beta1'], out['beta2'], out['sigma'], x, y, N)
print(f'r_hat summary | mean: {np.mean(r_hat):.4f}, std: {np.std(r_hat):.4f}')
print('-'*50)

【実行結果】
処理時間は 100 ミリ秒未満！

パラメータ推定値と標準誤差はテキストの結果とほぼほぼ一致しています。

100 個体の生存種子数の予測値を可視化します。
薄オレンジ線がたくさん描画されるあのチャートです。

# GLMM化したロジスティック回帰の推定にもとづく予測 p.161 図7.10(A) の改造

### 設定と準備
# 真の係数パラメータ
beta1, beta2 = -4, 1
# 調査種子数
N = 8

### 描画処理
## 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

## 設定
x_val_a = np.linspace(1, 7, 100)              # x軸の値
y_val_a = expit(beta1 + beta2 * x_val_a) * N  # 真の生存確率の場合の生存種子数y

## 描画
# 観測値の散布図の描画（jitter処理込み）
sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
            scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
# 真の生存確率の赤い点線の描画
ax.plot(x_val_a, y_val_a, color='tab:red', lw=2, ls='--', label='真の生存確率の例')
# GLMMの予測値の緑実線の描画 ※r=0である
ax.plot(x_val_a, expit(out['beta1'] + out['beta2'] * x_val_a + 0) * N,
        color='green', lw=2, label='予測値(個体差除く)')
# GLMMの予測値のオレンジ実線の描画
for r in est_r_hat:
   ax.plot(x_val_a, expit(out['beta1'] + out['beta2'] * x_val_a + r) * N,
           color='tab:orange', lw=0.5, alpha=0.3)
ax.plot(0, 0, color='tab:orange', lw=1, alpha=0.5, label='予測値(個体差含む)')
# 修飾
ax.set_xlabel('葉数 $x_i$', fontsize=12)
ax.set_ylabel('生存;種子数 $y_i$', fontsize=12)
ax.set(xlim=(1.6, 6.4), title='(A) 葉数と生存種子数の関係')
ax.legend(bbox_to_anchor=(1, 1), loc='upper left');

【実行結果】
いい感じです！

🍀🍀🍀

🔷 コードの概説
ChatGPT にコードの概要を訊きました。

この「ベクトル化 Gauss–Hermite（GH）で高速 MLE」のコードが“なにをしているか”を、概要 → 数式 → 数式の意味 → 数式とコードの対応の順でやさしく解説します。

１．GH法によるGLMMの概要（まず全体像）

• モデル：二項×正規ランダム切片・ロジットリンク
  個体（グループ）ii の成功回数

• 目的：周辺尤度（ランダム効果 $${r_i}$$ を積分消去）で $${\theta=(\beta_1,\beta_2,\sigma)}$$ を 最尤推定。

• 課題：周辺尤度は積分を含む $${\int \text{Binom} \times \text{Normal}\ dr}$$。
  解析解がないので数値積分が必要。

• 解決：Gauss–Hermite（GH）で積分を精度よく近似、さらにベクトル化で一気に計算 → 高速＆安定。

📐📐

２．数式（周辺尤度とGH近似）

2.1 周辺尤度
個体 $${i}$$ の尤度は

全体の対数尤度は $${\log L(\theta)=\sum_i \log L_i(\theta)}$$。
最小化対象は $${- \log L(\theta)}$$。

2.2 GH近似のコア
標準の GH は

（$${z_j, w_j}$$ は $${\mathtt{hermgauss(K)}}$$ のノードと重み）。

ここで $${r=\sqrt{2} \sigma z}$$ とおくと

よって

📐📐

３．数式の“意味合い”（なにをやっている？）

• ランダム切片 $${r}$$ の事前平均 0、分散 $${\sigma^2}$$ に従う期待値（周辺化）を“点の和”で近似している。
  つまり「$${r}$$ が取りそうな代表点 $${r_j = \sqrt{2} \sigma z_j}$$」での二項尤度を、「その点の重み $${w_j/\sqrt{\pi}}$$」で足し合わせ → 周辺尤度の近似。

• $${K}$$ を増やすほど精度は上がる（ただし計算量も増える）。
  $${\sigma}$$ が大きい／データが極端（$${y=0}$$ or $${N}$$ が多い）ほど、少し大きめの $${K}$$ が効く。

• 数値安定の要：log-sum-exp で $${\log\sum w_j \mathrm{pmf}}$$ を安全に計算、二項係数は $${ \log \Gamma }$$（$${\mathtt{gammaln}}$$）で。

📐📐

４．数式 ⇔ コードの対応（どの式がどの行？）
以下、コード断片を式に対応づけて読み解きます。

4.1 GHノードと重み（標準形 $${e^{-z^2}}$$）

nodes, weights = np.polynomial.hermite.hermgauss(K)

• 数式の $${\{z_j,w_j\}_{j=1}^K}$$ に対応（標準 GH）。

4.2 Normal(0,σ²) へのスケーリング

sigma = np.exp(log_sigma)                 # σ>0 を保証（再パラ）
rj = np.sqrt(2.0) * sigma * nodes         # r_j = √2 σ z_j
wj = weights / np.sqrt(np.pi)             # w_j' = w_j / √π
log_wj = np.log(wj + 1e-300)              # log(0) 回避の微小量

• 置換 $${r=\sqrt{2} \sigma z}$$ と $${f(r)dr=(1/\sqrt{\pi})e^{-z^2}dz}$$ をそのままコード化。

4.3 線形予測子とロジット確率（ベクトル化）

eta0 = beta1 + beta2 * x                  # 形 (n,)
eta  = eta0[:, None] + rj[None, :]        # 形 (n, K)：η_{i,j}=β1+β2 x_i + r_j
p    = expit(eta)                         # p_{i,j} = logistic(η_{i,j})

• 全グループ $${i=1, \ldots, n}$$ と全ノード $${j=1, \ldots, K}$$ を一発で作る（ブロードキャスト）。

4.4 二項の対数 pmf（安定に計算）

Ncol = N[:, None].astype(float) ; ycol = y[:, None].astype(float)
logC = gammaln(Ncol + 1) - gammaln(ycol + 1) - gammaln(Ncol - ycol + 1)
eps  = 1e-300
logpmf = logC + ycol*np.log(p + eps) + (Ncol - ycol)*np.log(1.0 - p + eps)

• $${\log{N\choose y} + y\log p + (N-y)\log(1-p)}$$

• $${\mathtt{gammaln}}$$ で $${\log\Gamma}$$ を使い、桁落ち回避。

4.5 $${\log \sum_j w_j \,\mathrm{pmf}}$$ を log-sum-exp で

log_terms = log_wj[None, :] + logpmf      # 各 (i,j) の log(w_j*pmf)
log_Li    = logsumexp(log_terms, axis=1)  # log Σ_j w_j*pmf (各 i)
return -np.sum(log_Li)                    # 負の対数尤度

• $${\log L_i \approx \log\sum_j w'_j\,\mathrm{pmf}_{i,j}}$$

• $${\mathtt{logsumexp}}$$ は $${\log\sum \exp(\cdot)}$$ を安定に計算。

4.6 MLE（最小化）と SE

obj = lambda th: nll_glmm_binom_logit_GH_vec(
    th, x, y, N, K=K, nodes=nodes, weights=weights)
res = optimize.minimize(fun=obj, x0=init, method='L-BFGS-B')
beta1_mle, beta2_mle, log_sigma_mle = res.x
sigma_mle = np.exp(log_sigma_mle)

• $${\theta=\arg\min_\theta\{-\log L(\theta)\}}$$。

• $${\sigma}$$ は $${\log \sigma}$$ で推定→元に戻す。

H = approx_hess(res.x, lambda th: nll_glmm_binom_logit_GH_vec(...))
cov = np.linalg.inv(H)
se_beta1, se_beta2, se_logsig = np.sqrt(np.diag(cov))
se_sigma = se_logsig * sigma_mle          # デルタ法：dσ/d(logσ)=σ

• 近似ヘッセ $${H\approx -\partial^2\log L/\partial\theta^2}$$ の逆で共分散。

• $${\text{se}(\sigma)\approx \sigma \cdot \text{se}(\log \sigma)}$$（デルタ法）。

📐📐

５．ここが効く！実務＆数値安定のコツ

• K の選び方：$${\sigma}$$ が大きいほど $${K}$$ を増やす（20〜120目安）。

• log-sum-exp と gammaln は必須：極端な $${y=0 , y=N}$$ でも安定。

• ベクトル化で $${O(nK)}$$ を Python ループ無しに回す → 桁違いに速い。

• $${\sigma}$$ は $${\log \sigma}$$ で推定（正の制約）。L-BFGS-B は滑らかな問題に強い。

数式だらけでした (^_^;;;;)

まとめ

今回は GLMM を実践しました。

🔷 確率分布、リンク関数、線形予測子、ランダム効果

🔷 Python 実装
statsmodels には GLMM 関数がありません。
スクラッチで今回モデルのパラメータ推定を実行しました。

📐📐

テキストの見どころ

7.5 節、7.6 節には GLMM のヒントが欠かれています。
ぜひテキストをお読みください！

🔷 7.5 節「現実のデータ解析にはGLMMが必要」
「疑似反復」のあるデータこそ GLMM が使えることが分かります。

• 同じ個体から複数のデータをとる
  ⇒個体差のランダム効果を推定可能

• 同じ植木鉢から複数のデータをとる
  ⇒植木鉢差のランダム効果を推定可能

• 同じ店舗から複数のデータをとる
  ⇒店舗差のランダム効果を推定可能
  etc.

🍀

🔷 7.6 節「いろいろな分布のGLMM」
目的変数のばらつきが二項分布以外の GLMM を紹介しています。

• ポアソン分布

• 負の二項分布

• ガンマ分布

• 正規分布（恒等リンク関数の場合は線形混合モデルLMM）

🍀

また、7.5 節、7.7 節では複雑化した GLMM のパラメータ推定が簡単でないことに触れられています。
個体差に加えて植木鉢差を考慮しながら最尤推定するのは「数値計算の問題として難しいものになる」ようです。
８章ではより具体的に「ランダム効果の発生源が $${K}$$ 個の場合、$${K}$$ 回の多重積分が必要とされ、$${K}$$ が大きいほど、最尤推定に要する時間が長くなったり、最尤推定値の探索自体が困難になります」と説明されています。

そこで次章からベイズ統計モデルを取り組もう！という流れが作られます。
「GLMM と最尤推定」から「階層ベイズモデルとMCMC」へ移行します。

今回のブログは以上です。

次回は一般化線形混合モデルを Python の既存ライブラリで「なんとかして」攻略します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。