はじめに
数理統計学より実践的な内容の理解に向けて，久保 拓弥 著「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。
本記事は，第3章「一般化線形モデル(GLM)」の読書メモである。なおPythonコードは，ChatGPT 3.5を用いて作成し，動作確認を行なっている。
- 本書の紹介ページ
- 本書のサポートページ
第3章 一般化線形モデル(GLM) ―ポアソン回帰―
3.1 例題 : 個体ごとに平均種子数が異なる場合
本節では，第2章と第3章の問題設定の違いを説明している。
興味がある変数は，いずれも種子数である。
第2章では，この種子数を決めるのはポアソン分布のパラメータのみであり，全ての種子に共通であった。
一方第3章では，体の大きさと施肥の有無
が追加されており，これらによって個体差が表現できるようになっている。
3.2 観測されたデータの概要を調べる
本節では，データの準備および基本統計量の計算のためのコードを説明している。
本書では，コードはR言語を用いているが，今回は同様の挙動をするPythonコードをGoogle Colab上で作成した。
コード
from google.colab import drive drive.mount('/content/drive')
出力
Mounted at /content/drive
本書のサポートページから'data3a.csv'を入手し，読み込む。
コード
import pandas as pd d = pd.read_csv('/content/drive/My Drive/Colab Notebooks/StatModeling/data3a.csv') print(d.head())
出力
y x f 0 6 8.31 C 1 6 9.44 C 2 6 9.50 C 3 12 9.07 C 4 10 10.16 C
3.3 統計モデリングの前にデータを図示する
本節では，データの可視化の例として，体の大きさと種子数
の散布図を紹介している。この散布図に，さらに施肥の有無の情報も追加している。
コード
import matplotlib.pyplot as plt # 散布図を作成 plt.scatter(d[d['f'] == 'C']['x'], d[d['f'] == 'C']['y'], color='white', edgecolor='black', label='C') # 'C'に対して白い丸 plt.scatter(d[d['f'] == 'T']['x'], d[d['f'] == 'T']['y'], color='black', label='T') # 'T'に対して黒い丸 # ラベルと凡例を追加 plt.xlabel('x') plt.ylabel('y') plt.legend() # グラフを表示 plt.show()
出力
次に，箱ひげ図を描く。
コード
grouped_data = [d[d['f'] == 'C']['y'], d[d['f'] == 'T']['y']] # 箱ひげ図を描画 plt.figure(figsize=(8, 6)) plt.boxplot(grouped_data, labels=['C', 'T']) # ラベルを追加 plt.xlabel('f') plt.ylabel('y') # グラフを表示 plt.show()
出力
3.4 ポアソン回帰の統計モデル
本節では，本書のテーマの一つである一般化線形モデルが導入されている。
一般化線形モデルのうち，カウントデータを表現するための統計モデルとして，ポアソン回帰モデルが説明されている。
ポアソン回帰モデルでは，まず種子数の確率をポアソン分布
にしたがうと仮定する。
ポアソン回帰モデルでは，さらに線形予測子(linear predictor)とリンク関数(link function)を導入する。
ポアソン分布のパラメータを，説明変数
をもちいて
と表現する。この右辺を線形予測子と呼ぶ。この線形予測子を用いることで，ポアソン回帰のパラメータに，個体ごとの情報である
を含めることができる。
また，この左辺の関数を，リンク関数と呼ぶ。この例では対数関数が使われているので，対数リンク関数と呼ぶ。
対数リンク関数を用いることで，ポアソン分布のパラメータについてが満たされる。
ポアソン回帰では，観測データを用いて，ポアソン分布を用いた統計モデルのパラメータの推定を行なう。Pythonでポアソン回帰を行なうためには，statmodelsパッケージを用いる。
コード
import numpy as np import statsmodels.api as sm # 説明変数と目的変数を定義 X = d['x'] y = d['y'] # 説明変数に定数項を追加 X = sm.add_constant(X) # ポアソン回帰モデルをフィット poisson_model = sm.GLM(y, X, family=sm.families.Poisson()).fit() # モデルの要約を表示 print(poisson_model.summary())
出力
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -235.39
Date: Tue, 15 Oct 2024 Deviance: 84.993
Time: 13:34:09 Pearson chi2: 83.8
No. Iterations: 4 Pseudo R-squ. (CS): 0.04414
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 1.2917 0.364 3.552 0.000 0.579 2.005
x 0.0757 0.036 2.125 0.034 0.006 0.145
==============================================================================推定された回帰係数には，重回帰分析の出力のように，それぞれWald統計量の情報が示されている。最尤推定量が漸近的に正規分布にしたがうことから，信頼区間が算出されている，と解釈できる。
なお，最大対数尤度の値は出力の中に含まれており，
Log-Likelihood: -235.39
と表示されている。
3.5 説明変数が因子型の統計モデル
3.4節では，施肥効果を用いていなかったが，本節ではこれを説明変数として組み込んだモデルを説明している。
この時のモデルは，
\begin{equation}
d_i=
\begin{cases}
0 & \text{if $f_i = C $ (施肥なし)} \\
1 & \text{if $f_i = T $ (施肥あり)}
\end{cases}
\end{equation}
である。
施肥効果は，値がCまたはTである2水準のカテゴリ変数であるため，ダミー変数に置き換える必要がある。
コード
# カテゴリ変数fをダミー変数に変換 # drop_first=Trueにより、'C'を基準カテゴリとして'T'のみをダミー変数にします X = pd.get_dummies(d['f'], drop_first=True) # 説明変数をfloat型に変換 X = X.astype(float) # 定数項（切片）を追加 X = sm.add_constant(X) # 目的変数 y = d['y'] # ポアソン回帰モデルをフィット poisson_model = sm.GLM(y, X, family=sm.families.Poisson()).fit() # モデルの要約を表示 print(poisson_model.summary())
出力
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -237.63
Date: Tue, 15 Oct 2024 Deviance: 89.475
Time: 14:05:20 Pearson chi2: 87.1
No. Iterations: 4 Pseudo R-squ. (CS): 0.0003192
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 2.0516 0.051 40.463 0.000 1.952 2.151
T 0.0128 0.071 0.179 0.858 -0.127 0.153
==============================================================================施肥効果を説明変数にした場合，
を説明変数にしたときと比べて，最大対数尤度が小さくなっている。
が説明変数の場合 : Log-Likelihood: -235.39
が説明変数の場合 : Log-Likelihood: -237.63
これは，施肥効果を説明変数にしたモデルの方が当てはまりが悪いためと考えられる。
3.6 説明変数が数量型+因子型の統計モデル
本節では，個体の体サイズと施肥効果
の複数の説明変数を同時に組込んだ統計モデルを作っている。
コード
# カテゴリ変数fをダミー変数に変換 # drop_first=Trueにより、'C'を基準カテゴリとして'T'のみをダミー変数にします X = pd.get_dummies(d[['x', 'f']], drop_first=True) # 説明変数をfloat型に変換 X = X.astype(float) # 定数項（切片）を追加 X = sm.add_constant(X) # 目的変数 y = d['y'] # ポアソン回帰モデルをフィット poisson_model = sm.GLM(y, X, family=sm.families.Poisson()).fit() # モデルの要約を表示 print(poisson_model.summary())
出力
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -235.29
Date: Tue, 15 Oct 2024 Deviance: 84.808
Time: 21:45:52 Pearson chi2: 83.8
No. Iterations: 4 Pseudo R-squ. (CS): 0.04590
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 1.2631 0.370 3.417 0.001 0.539 1.988
x 0.0801 0.037 2.162 0.031 0.007 0.153
f_T -0.0320 0.074 -0.430 0.667 -0.178 0.114
==============================================================================
