2022-12-23

一般化線形モデル

Statistics

Statistical Model

Python

一般化線形モデル（GLM）とは

一般化線形モデル (Generalized Linear Model; GLM）とは、目的変数が正規分布以外の確率分布に従う場合であっても適用できる統計モデルです。

単回帰分析や重回帰分析は線形モデルと呼ばれるもので、目的変数が正規分布に従うと仮定した統計モデルです。目的変数が正規分布に従うと仮定できない場合、例えば目的変数がカテゴリカル、カウントデータのような離散変数の場合は線形モデルは適用できません。GLMでは正規分布だけでなく、二項分布、ガンマ分布、ポアソン分布といった指数型分布族と呼ばれる分布も扱うことができます。

GLMでは、目的変数と説明変数との関係式は簡単な線形式である必要はなく、以下のように表すことができます。

g(y) = \beta_0x_0 + \beta_1x_1 + \cdots + \beta_nx_n

$\beta$ はパラメータ、 $x$ は説明変数、 $n$ は説明変数の数です。上式の右辺は線形結合となっており、線形予測子 linear predictorと呼ばれます。対して右辺が線形予測子であるときの左辺の関数 $g(y)$ はリンク関数 link functionと呼ばれます。

確率分布とリンク関数の組み合わせ

よく使われる確率分布とリンク関数の組み合わせは以下になります。

	確率分布	リンク関数
離散	二項分布	logit
離散	ポアソン分布	log
離散	負の二項分布	log
連続	ガンマ分布	log
連続	正規分布	identity

GLM の構築

実際にGLMを構築してみます。ここでは、ある架空の植物100個体の種子数データについて考えます。データは以下のカラムで構成されています。

y: 個体の種子数
x: 個体のサイズ
f: 個体に施肥処理をしたかどうか（C: 施肥処理をしていない、T: 施肥処理をしている）

import pandas as pd

# number of seeds
y = [6,6,6,12,10,4,9,9,9,11,6,10,6,10,11,8,3,8,5,5,4,11,5,10,6,6,7,9,3,10,2,9,10,
     5,11,10,4,8,9,12,8,9,8,6,6,10,10,9,12,6,14,6,7,9,6,7,9,13,9,13,7,8,10,7,12,
     6,15,3,4,6,10,8,8,7,5,6,8,9,9,6,7,10,6,11,11,11,5,6,4,5,6,5,8,5,9,8,6,8,7,9]

# body size
x = [8.31,9.44,9.5,9.07,10.16,8.32,10.61,10.06,9.93,10.43,10.36,10.15,10.92,8.85,
     9.42,11.11,8.02,11.93,8.55,7.19,9.83,10.79,8.89,10.09,11.63,10.21,9.45,10.44,
     9.44,10.48,9.43,10.32,10.33,8.5,9.41,8.96,9.67,10.26,10.36,11.8,10.94,10.25,
     8.74,10.46,9.37,9.74,8.95,8.74,11.32,9.25,10.14,9.05,9.89,8.76,12.04,9.91,
     9.84,11.87,10.16,9.34,10.17,10.99,9.19,10.67,10.96,10.55,9.69,10.91,9.6,
     12.37,10.54,11.3,12.4,10.18,9.53,10.24,11.76,9.52,10.4,9.96,10.3,11.54,9.42,
     11.28,9.73,10.78,10.21,10.51,10.73,8.85,11.2,9.86,11.54,10.03,11.88,9.15,8.52,
     10.24,10.86,9.97]

# feed (C: control, T: treatment)
f = ['C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C',
     'C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C',
     'C','C','C','C','C','C','C','C','C','C','C','C','T','T','T','T','T','T','T',
     'T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T',
     'T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T','T',
     'T','T','T','T','T']

df = pd.DataFrame(list(zip(y, x, f)), columns =['y', 'x', 'f'])

display(df.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

データの確認

記述統計量を確認します。

df.describe()

	y	x
count	100.000000	100.000000
mean	7.830000	10.089100
std	2.624881	1.008049
min	2.000000	7.190000
25%	6.000000	9.427500
50%	8.000000	10.155000
75%	10.000000	10.685000
max	15.000000	12.400000

fについてC、Tの数をそれぞれ表示します。

print('num of C:', len(df[df['f'] == 'C']))
print('num of T:', len(df[df['f'] == 'T']))

num of C: 50
num of T: 50

データの記述統計量から以下のことが確認できます。

種子数y:
- 非負の整数となっている（2から15）
- 平均と分散がまあまあ近い
個体サイズx:
- 正の実数（連続値）となっている
- 分散は小さい
施肥処理の有無f:
- カテゴリ変数となっている
- 各カテゴリの要素数は均等になっている

次にデータの可視化を行います。まずは散布図を表示します。

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

plt.style.use('ggplot')
plt.figure(figsize=(10,5))
sns.set(rc={'figure.figsize':(10,5)})

sns.scatterplot(x='x', y='y', hue='f', data=df)

Scatter of f

箱ヒゲ図を表示します。

sns.boxplot(x='f', y='y', data=df)

Box of f

グラフから以下のことが確認できます。

サイズ $x$ が増加するにつれて種子数yも増加している（散布図）
肥料の効果は無いように見える（箱ヒゲ図）

ポアソン回帰で統計モデルの構築

ある架空の植物の個体 $i$ の種子数 $y_i$ はポアソン分布に従うと仮定します。ポアソン分布の確率分布は以下の式になります。

P(X=y| \lambda) = \frac{e^{-\lambda} \lambda^y_i}{y_i!}

ポアソン分布には分布の形状を決めるパラメータ $\lambda$ がありますが、 $\lambda$ （平均種子数）は個体ごとに異なり、個体ごとに異なる形状をしたポアソン分布から種子数が決まると考えます。今回は以下の3つのパターンの統計モデルを構築します。

種子数が個体のサイズに依存する統計モデル
種子数が施肥処理の有無に依存する統計モデル
種子数がサイズと施肥処理の有無に依存する統計モデル

種子数が個体のサイズに依存する統計モデル

個体ごとに異なる $\lambda_i$ を、まずは個体サイズ $x_i$ の関数として以下のように定義してみます。

\lambda_i = e^{\beta_1 + \beta_2x_i}

ここで、 $\beta_1$ や $\beta_2$ をパラメータ、もしくは係数といい、 $\beta_1$ は切片（intercept）といいます。

両辺対数を取ると以下のようになります。

\log (\lambda_i) = \beta_1 + \beta_2 x_i

このとき、右辺 $\beta_1 + \beta_2 x_i$ は線形結合になっているため、線形予測子 linear predictorと呼ばれます。

線形予測子に対する左辺の関数はリンク関数 link functionと呼ばれます。今回のリンク関数は $\log (\lambda_i)$ であり、対数となっているため、対数リンク関数 log link functionと呼ばれます。

ポアソン回帰のGLMではよく対数リンク関数が使われます。その理由としては、ポアソン分布の平均 $\lambda$ は非負値である必要があるためです。 $\lambda_i = e^{linear\_predictor} \geqq 0$ となり非負値であるため都合が良いということになります。

これから以下の対数尤度関数を最大化する $\beta_1$ 、 $\beta_2$ を推定します。

\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i}e^{-\lambda_i}}{y_i!}

パラメータの推定は以下のPythonコードで簡単に行うことができます。

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

fit_x = smf.glm('y ~ x', data=df, family=sm.families.Poisson()).fit()
fit_x.summary()

Poisson regression | x

以下の結果が得られました。

	coef	std err	z	Log-likelihood
Intercept	1.2917	0.364	3.552
x	0.0757	0.036	2.125
				-235.39

ここで、std errは標準誤差のことで、つまりは推定したパラメータの標準偏差になります。標準偏差の導出には正規分布を仮定しています。また、zは z 値のことです。

$\lambda$ （平均種子数）は以下のようになります。

\lambda_i = e^{1.2917 + 0.0757x_i}

上記の関数をプロットしてみます。

import numpy as np
import math

sns.scatterplot(x='x', y='y', hue='f', data=df)
x = np.arange(df.x.min(), df.x.max() + 0.01, 0.01)
y = [math.exp(fit_x.params['Intercept'] + fit_x.params['x'] * x_i)  for x_i in x]
plt.plot(x, y)

Predict | x

サイズ $x$ の増加につれて平均種子数が増加するという線が描かれています。

種子数が施肥処理の有無に依存する統計モデル

施肥の有無 $f_i$ を説明変数として組み込んだモデルを検討します。

施肥の有無 $f_i$ はカテゴリ変数となっています。カテゴリ変数は通常ダミー変数に置き換えます。施肥の有無 $f_i$ をダミー変数 $d_i$ に置き換えると以下のようになります。

d_i = \left\{ \begin{array}{ll} 0 & (f_i = \mathrm{C}) \\ 1 & (f_i = \mathrm{T}) \end{array} \right.

モデルは以下のようになります。

\lambda_i = e^{\beta_1 + \beta_3 d_i}

つまり、個体 $i$ が施肥処理なし（ $f_i=C$ ）の場合は

\lambda_i = e^{\beta_i}

となり、個体 $i$ が施肥処理あり（ $f_i=T$ ）の場合は

\lambda_i = e^{\beta_i + \beta_3}

となります。

以下の対数尤度関数を最大化する $\beta_1$ 、 $\beta_2$ を推定します。

\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i}e^{-\lambda_i}}{y_i!}

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

fit_f = smf.glm('y ~ f', data=df, family=sm.families.Poisson()).fit()
fit_f.summary()

Poisson regression | f

以下の結果が得られました。

	coef	std err	z	Log-likelihood
Intercept	2.0516	0.051	40.463
f[T.T]	0.0128	0.071	0.179
				-237.63

つまり、 $\lambda$ （平均種子数）は以下のようになります。

\lambda_i = e^{2.0516} \quad (f_i = \mathrm{C})

\lambda_i = e^{2.0516 + 0.0128} \quad (f_i = \mathrm{T})

施肥処理を加えると平均種子数が少しだけ増えるという結果となりました。

log-likelihoodは-237.63となっており、 $x_i$ を説明変数としたモデルのlog-likelihood -235.39よりも小さくなっています。つまり、モデルのデータに対する当てはまりが悪くなっています。

種子数がサイズと施肥処理の有無に依存する統計モデル

最後に、サイズ $x_i$ と施肥の有無 $f_i$ を説明変数として組み込んだモデルを検討します。

モデルは以下のようになります。

\lambda_i = e^{\beta_1 + \beta_2x_i + \beta_3 d_i}

以下の対数尤度関数を最大化する $\beta_1$ 、 $\beta_2$ 、 $\beta_3$ を推定します。

\log L(\beta_1, \beta_2, \beta_3) = \sum_i \log \frac{\lambda_i^{y_i}e^{-\lambda_i}}{y_i!}

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

fit_full = smf.glm('y ~ x + f', data=df, family=sm.families.Poisson()).fit()
fit_full.summary()