データ解析のための統計モデリング入門 7章 教科書中の図をPythonで再現してみた

はじめに

個体差などを考慮してデータをモデリングするためには、一般化線形モデルを拡張して、**一般化線形混合モデル（Generalized Linear Mixed Model, 以降GLMMと略します）**などを使う必要性が生じます。
GLMMは、例えばデータ解析のための統計モデリング入門の７章に説明があります。GLMMの考え方、およびRを使ったパラメータ推定方法はこの教科書にわかりやすく書いてあります。

ただし、私は一度読んだだけでは、混合モデルの分布形状やパラメータ推定後のモデルの予測値のプロットをすんなりとは理解できませんでした。この記事では、Pythonを使って、数値積分を実施して混合モデルの分布形状の図を再現する方法を記載します。なお、最尤法によるGLMMのパラメータ推定はRでしかできないようですので、パラメータ推定結果はRの結果をそのままPythonに転記して使うことにします。

データとGLMによるパラメータ推定

データは教科書と同じ、架空植物の各個体から8個の種子を取ってきて、いくつが生存しているかを確認したものを使います。データは教科書のサポートページにあります。「kubobook_2012.zip」を解凍すると、「glmm」というフォルダに「data.csv」が格納されており、これが今回扱うデータとなります。データにはxとして葉数、yとして生存している種子数(0-8)が入力されています。これをpandasのデータフレームとして読み込みます。

# ライブラリの読み込み
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# GLMのためインポート
import statsmodels.api as sm

## 二項計算のためインポート
import scipy.special

df = pd.read_csv("./kubobook_2012/glmm/data.csv")

#後で使うため、x,yごとのデータ件数も求めておく
df_count = pd.DataFrame(df.groupby(["x", "y"]).size(), columns=["counts"]).reset_index()

まずは説明変数である葉数と目的変数である生存種子数の関係を可視化します。ただし、このデータは同一の葉数と生存種子数に複数のデータがあります。散布図をそのまま描いたのでは重なってしまい、どこにデータが多く分布しているかわかりません。
解決方法の一つはalphaを使って１つの点に透明度を持たせて、データが多く重複している箇所は色が濃くするようにすることです。

plt.scatter(df.x, df.y, alpha=0.5)
plt.colorbar()
plt.xlabel("葉数")
plt.ylabel("生存種子数")
plt.show()

もう一つの解決方法は、seabornのstripplotを使い、jitter=Trueを指定してデータ点をランダムにずらして、重ならないように表示することです。

sns.stripplot(x="x", y="y", data=df, jitter=True, color="white", edgecolor="black", linewidth=1)
plt.xlabel("葉数")
plt.ylabel("生存種子数")
plt.show()

GLMによるモデリング

生存種子数は二項分布に従い、葉数 $x_i$ が与えられると生存確率 $q_i$ は以下のロジスティック関数で与えられるとします。
$ q_i = \frac{1}{1 + \exp(-\beta_1 - \beta_2 * x_i)} $

Pythonでパラメータ $\beta_1$, $\beta_2$を推定するには、statsmodelsのGLMを使います。このGLMで二項分布を扱うときは、目的変数は（生存数、死亡数）の2つを渡す必要がありますので、死亡数を計算しておきます。また、切片を評価するためには、定数(const）の列も追加しておく必要があり、statsmodelsのadd_constant()を使います。

# GLMには（生存数、死亡数））で渡す必要があるので、死亡数N-niを計算する
df[["N-y"]] = df.N - df.y
df_y = df[["y", "N-y"]]

# 切片 beta_1を評価するため、定数の列を追加する
df_x = df[["x"]]
df_x = sm.add_constant(data=df_x)

binomial_model = sm.GLM(df_y, df_x, family=sm.families.Binomial()) 

result = binomial_model.fit()
display(result.summary())

以下のように推定されたパラメータなどの結果がでます。

予測結果とデータを比較します。予測は predict()メソッドで実行できます。 また、この教科書では、$\beta_1 = -4$, $\beta_2 = 1$ を真の値としてデータを作成しています。真の分布も同時に比較します。

# 葉数が与えられたとき、真の分布を返す関数
def true_plot(x):
    beta1 = -4
    beta2 = 1
    q = 1/(1 + np.exp(-beta1 -beta2*x))
    return N*q

# 葉数の範囲
rng_x = [x for x in range(2,7)]
# [定数項、葉数]データ
var_x = [ [1,x] for x in rng_x ]

# データのプロット
plt.scatter(df.x, df.y, alpha=0.5)
# stripplotではうまくいかない（理由は調査中）
# sns.stripplot(x=df.x, y=df.y, data=df, jitter=True, color="white", edgecolor="black", linewidth=1)

#GLM予測の結果
plt.plot( rng_x, result.predict(var_x)*N, color="black", label= "GLMの予測")
#真の分布
plt.plot( rng_x, [true_plot(x) for x in range(2,7)], color="green", linestyle="dashed", label="真の分布")

plt.colorbar()
plt.xlabel("葉数")
plt.ylabel("生存種子数")
plt.legend()
plt.show()

GLMの予測結果は真の分布と比較して、直線になっていてデータの特徴を柔軟に捉えられていないようです。

個体差と過分散

上のように、GLMではデータのばらつきをうまくとらえられず、真の分布よりも直線的になっていました。
葉数$x=4$のときについて、もう少し詳しく見てみます。葉数$x=4$のときの生存数分布について、データと推定されたGLMを比較します。GLMの分布を描くとき、パラメータの値は.params.constなどを使って取得します。

x = 4
q = 1/(1+np.exp(-result.params.const - result.params.x * x))

# 葉数x=4のサンプル数20を取得
Nsample_x4 = df.query("x==4").shape[0]

# x=4のときGLMが予測する生存数分布
list_pred_x4 = []

for ni in range(0,9):
    list_pred_x4.append( scipy.special.binom(N, ni) * q**ni * (1-q)**(N-ni) * Nsample_x4 )
    
list_pred_x4

できたグラフはこちら。

GLMの予測は二項分布の形状をしていますが、データはもっとばらついて、生存種子数が2、および7,8あたりにピークがあります。葉数以外にも、個体差があって生存種子数にばらつきが生じていて、データでは分散が大きくなっているようです。これを過分散といいます。GLMでは葉数が決まれば種子が生存する確率も決まると仮定していましたが、この過分散があるデータではこの仮定が成り立ちません。データへのあてはめをよくするには、モデルを拡張する必要があります。

一般化線形混合モデル（GLMM）の分布

ここで一般化線形混合モデルの登場となります。上の生存確率を表すロジスティック関数に、個体差を表す$r_i$を加えます。

$ q_i = \frac{1}{1 + \exp(-\beta_1 - \beta_2 * x_i - r_i)} $

$r_i$はサンプル数（ここでは100）ありますが、100個のデータから$\beta_1$、$\beta_2$と$r_i (i=1-100)$を推定するのは、最尤法の枠組みではナンセンスです。$r_i$も統計モデリングで扱うとこにします。デーiタが得られていない個体差の分布はわからないのですが、ここでは扱いやすい正規分布を仮定します。

$p(r_i|s) = \frac{1}{\sqrt{2\pi}s} \exp(-\frac{r_i^2}{2s^2})$

$\beta_1$, $\beta_2$, $s$が与えられた時、生存種子数$y_i$を得る確率は、

$p(y_i|\beta_1, \beta_2, s) = \int p(y_i|\beta_1, \beta_2, r_i)p(r_i|s) dr_i$
で与えらえれます。ここで、$p(y_i|\beta_1, \beta_2, r_i)$は二項分布$_N C_n q_i^{n_i} (1-q_i)^{N-n_i}$です。
二項分布と正規分布の積を積分することで、確率が得られます。

GLMMの分布を確認する

個体差をモデル化して二項分布と混合することで、生存種子数の分布がどのように変化するのかを確認します。$\beta_1$と$\beta_2$は真の値-4と1を使います。個体差のばらつきの大きさを表す標準偏差$s$を変化させて、どのように変化するか数値積分して確認します。数値積分はscipyのintegrate.quad()を使います。

#種子数など、パラメータの真の値などの設定
N=8
#ni = 4
beta1 = -4 #定数
beta2 = 1 #葉数の係数
xi = 4 #端数 ここでは4を扱う

s=3 #個体差の正規分布の標準偏差

# 葉数と個体に依存する二項分布：GLMMでriを積分するため、葉数は別で決めておき、個体差riの関数として扱う
def Binom(ri):
    qi = 1/( 1 + np.exp(-beta1 -beta2*xi -ri)) 
    #print(qi)
    #print(beta1, beta2)
    return scipy.special.binom(N, ni) * qi**ni * (1-qi)**(N-ni)

# 個体差が従う正規分布
def Gauss(ri):
    #print(s)
    return np.exp(-ri**2/(2*s**2)) /np.sqrt((2*np.pi)) /s

# 二項分布と個体差の正規分布の積：これをriで積分したものが生存種子数の分布を与える
def Prob(ri):
    return Binom(ri)*Gauss(ri)

## 真のパラメータを使って、個体差の標準偏差sを変えたときの分布を描く
# 標準偏差sの範囲
list_s = [0.1, 0.5, 1, 1.5, 2, 3, 5]

list_result_all = []

for s in list_s:

    list_result = []

    # 標準偏差sを与えた時の、生存種子数分布
    for ni in range(0,9):

        list_result.append( integrate.quad(Prob, -100, 100)[0]*Nsample_x4 )
    
    list_result_all.append(list_result)

# 生存種子数の分布を描く
for _i, _s in enumerate(list_s):
    plt.plot(list_result_all[_i], label="s=" + str(_s))
    plt.ylim(-0.5, 6.05)

plt.xlabel("生存種子数")
plt.ylabel("観測された個体数")    
plt.legend()
plt.show()

コードが少し長くなりましたが、結果を確認してみましょう。

$s=0.1$の時は二項分布に近い形をしています。これから$s$を大きくしていき$s=1.5$になると、ほぼ一様な分布形状になります。さらに大きくしていき、$s=3$を超えると両端が大きく中央の生存種子数=4あたりで最少となる、二項分布とは逆の形状をするようになります。標準偏差$s$の値を適切に推定できれば、上のデータへの当てはまりも良くなりそうです。

GLMM推定結果の確認

GLMMを最尤法で推定するパッケージはpythonではないようです。Rを使った推定は教科書に記載されていますので割愛します。（私が推定した値は教科書の値とは一致していません。ただし、こちらでも私と同じ値になっています。教科書の著者が使ったデータと、私たちが入手したデータが少し違うのでしょうか？解明できていませんが、影響は軽微なのでこのまま進めます。）

推定されたパラメータを使って、モデルの予測を確認します。教科書の脚注には「図7.10に示されているような予測をするためには、推定結果を組み合わせた数値積分が必要です。第10章の10.3.2を参照してください」と記載されています。数値積分ならPythonでもできますので、上のGLMMの分布を描くやり方で、推定されたパラメータを使えば予測はできるはずです。

# Rで推定されたパラメータ
beta1 = -4.190
beta2 = 1.005

s = 2.408

list_estimated = []

for ni in range(0,9):

    list_estimated.append( integrate.quad(Prob, -100, 100)[0]*Nsample_x4 )
    
list_estimated

plt.plot( df_count.query("x==4").y, df_count.query("x==4").counts, marker="o", label="データ", linewidth=0) 
plt.plot(list(range(0,9)), list_estimated, label="GLMMの予測")
plt.ylim(-0.5,6.5)
plt.xlabel("生存種子数")
plt.ylabel("観測された個体数")
plt.legend()

教科書の図7.10をほぼ再現した結果を描けました。推定された分散$s=2.408$は真の値よりは小さいものの、GLMと比較するとデータへのあてはめはだいぶ良くなりました。

葉数を変えた時のGLMMの予測と真の分布の比較結果も求めると下図のようになりました。GLMと比較すれば、だいぶデータへの当てはまりがよくなっています。

最後に

この記事では、一般化線形混合モデルの確率分布の振る舞いを調べるため、pythonを使って数値積分を実行して混合モデルの確率分布を求めました。実際に自分で$s$を変化させて確率分布がどのように変わるかを確認すると、だいぶ理解が深まりました。
一応検算などはしたつもりですが、もしこの記事に誤りなどがありましたらお知らせいただきたいです。

Pythonでは簡単にGLMMのパラメータを最尤法で推定できないのが残念ですが、Rを使えということなのでしょうか。
Pythonでも、ベイズ的な扱いはできそうなので、教科書の８章以降も読んで、Pythonで扱えるようにしていきます。