Traffine I/O

日本語

2022-12-23

一般化線形モデル

一般化線形モデル(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は標準誤差のことで、つまりは推定したパラメータの標準偏差になります。標準偏差の導出には正規分布を仮定しています。また、zz 値のことです。

\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()

Poisson regression | f

次の結果が得られました。

coef std err z Log-likelihood
Intercept 1.2631 0.370 3.417
f[T.T] -0.0320 0.074 -0.430
x 0.0801 0.037 2.162
-235.29

つまり、\lambda(平均種子数)は次のようになります。

\lambda_i = e^{2.0516 + 0.0801 x_i} \quad (f_i = \mathrm{C})
\lambda_i = e^{2.0516 + 0.0801 x_i - 0.0320} \quad (f_i = \mathrm{T})

施肥処理を行うと平均種子数が少なくなると解釈されています。具体的には、e^{−0.032} \approx 0.969 より施肥処理を行うと平均種子数が0.969倍になることになります。

log-likelihoodは-235.29となっており、これまでの2つのモデルよりもデータに対する当てはまりが良くなっています。

参考

https://kuboweb.github.io/-kubo/ce/IwanamiBook.html
https://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/49477/4/kubostat2008c.pdf

Ryusei Kakujo

researchgatelinkedingithub

Focusing on data science for mobility

Bench Press 100kg!