統計モデリング概論 DSHC 2023

1. 導入
2. 直線回帰、確率分布、擬似乱数生成
3. 尤度、最尤推定
4. 一般化線形モデル(GLM)
5. 個体差、一般化線形混合モデル(GLMM)
6. ベイズの定理、事後分布、MCMC
7. StanでGLM
8. 階層ベイズモデル(HBM)

ちょっとずつ線形モデルを発展させていく

線形モデル LM (単純な直線あてはめ)

    ↓ いろんな確率分布を扱いたい

一般化線形モデル GLM

    ↓ 個体差などの変量効果を扱いたい

一般化線形混合モデル GLMM

    ↓ もっと自由なモデリングを！

階層ベイズモデル HBM

データ解析のための統計モデリング入門 久保拓弥 2012 より改変

確率分布に長い時間を割いたけど、元はと言えば回帰したいのでした。

ここまでに見た統計モデル

確率変数$X$はパラメータ$\theta$の確率分布$f$に“従う”:  $X \sim f(\theta) $

e.g., ある植物が作る種の数$X$は平均値$\lambda$のポアソン分布に従う:

これを一般化線形モデル(GLM)として見ることもできる→

一般化線形モデル(GLM)として記述してみる

個体$i$の種子数$y_i$は平均値$\lambda_i$のポアソン分布に従う。
平均値$\lambda_i$は他のデータによらず$\beta_0$で一定。

種子数をY軸にして、式を2つに分けただけ…?
説明変数を含むモデルを見ればご利益が分かるかも。

説明変数が1つある一般化線形モデル

個体$i$の種子数$y_i$は平均値$\lambda_i$のポアソン分布に従う。
平均値の対数$\log(\textcolor{#3366ff}{\lambda_i})$はその個体の大きさ$x_i$に比例する。

この場合は単回帰。説明変数が複数あると重回帰。

複数の説明変数を同時に扱う重回帰

\[\begin{split} y_i &\sim \text{Poisson}(\lambda_i) \\ \log(\lambda_i) &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \ldots \end{split}\]

気温も湿度も高いほどビールが売れる架空データ:

ほかの確率分布とリンク関数を使う例を見てみよう。

ロジスティック回帰

• 確率分布: 二項分布
• リンク関数: $\operatorname{logit}(p) = \log \frac {p} {1 - p}$

何かの成否に対する何かの因子の影響、とか

ロジスティック回帰 (狭義)

• 確率分布: ベルヌーイ分布 ($n = 1$ の二項分布)
• リンク関数: $\operatorname{logit}(p) = \log \frac {p} {1 - p}$

何かの成否に対する何かの因子の影響、とか

一般線形モデル (“化”無し) はGLMの一種

• 確率分布: 正規分布
• リンク関数: 恒等関数(なにもせずそのまま)

最小二乗法の直線あてはめと結果的に同じになる。

単回帰・重回帰と言ったとき一般線形モデルを前提とする人もいる。

分散分析 (Analysis of variance, ANOVA) as GLM

質的な説明変数を持つ正規分布・恒等リンクのGLM、と解釈可能。
指示変数 (0 or 1) に変換してから重回帰する。

くもり☁️ $\beta_0$ を基準に、晴れの効果☀️ $\beta_1$ と雨の効果☔️ $\beta_2$ が求まる。

GLMなら確率分布・リンク関数を変えてもっと柔軟にモデリングできる。

共分散分析 (Analysis of covariance, ANCOVA) as GLM

質的変数と量的変数を両方含むGLM、と解釈可能。
正規分布・等分散・恒等リンクなどが仮定される。

GLMなら確率分布・リンク関数を変えてもっと柔軟にモデリングできる。

交互作用

ある説明変数の効果が、別の説明変数によって異なる。
e.g., ビール売上の温度依存性が天気によって異なる。

解釈が一気に難しくなるのでむやみに使わない。

一般化線形モデル(GLM)ふりかえり

確率分布・リンク関数を変えて柔軟にモデリングできる。
特定の組み合わせには名前がある。

確率分布については前章を参照。
リンク関数をもう少しだけ掘り下げたい。

リンク関数

統計モデリングにおいて「まっすぐ以外も表現できる」意味

$\text{identity}(\mu_i)$

$\mu_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \ldots$

説明変数の効果が足し算的に働く。

$\log(\lambda_i)$

$\lambda_i = e^{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \ldots} = e^{\beta_0} \times e^{\beta_1 x_{1i}} \times e^{\beta_2 x_{2i}} \times \ldots$

説明変数の効果が掛け算的に働く。
e.g., $\Delta x_1$ 増えると $e^{\beta_1 \Delta x_{1}}$ 倍になる

$\operatorname{logit}(p_i)$

$p_i = \frac 1 {1 + e^{-(\beta_0 + \beta_1 x_i + \ldots)}} $ (ロジスティック関数)

説明変数の効果が頭打ちになる。
e.g., $\lim_{x \to -\infty} p = 0;~\lim_{x \to \infty} p = 1$

ほかに probit, inverse, sqrt, etc.

statsmodelsにおけるGLMのやりかた

smf.glm の使い方は直線回帰のOLSとほぼ同じ

import statsmodels.api as sm
import statsmodels.formula.api as smf
formula = "weight ~ height"
model = smf.glm(formula, data=r.df_weight)
result = model.fit()
print(model.family)

<statsmodels.genmod.families.family.Gaussian object at 0x1270537d0>

何も指定しない場合は正規分布・恒等リンク。
family= オプションで 確率分布 と リンク関数 を明示的に指定:

identity = sm.families.links.Identity()
gaussian = sm.families.Gaussian(link=identity)
model = smf.glm(formula, data=r.df_weight, family=gaussian)

🔰 とにかくGLMを使ってみる練習

🔰 4-glm.ipynb をJupyterで開き、順々に実行してみよう。
ここまでに登場した回帰分析のPythonコードが書いてあります。

とりあえず当てはめと作図だけ。
結果の解釈やモデルの評価はこの後。

☕️ 休憩 + 質疑応答

データはひとつ、モデルはたくさん

どう選ぶ？

1. メカニズム的に納得できるものを選ぶ
   • ポアソン過程のカウントならポアソン分布、間隔ならガンマ分布
   • n回中k回のように割合的なカウントなら二項分布
2. データを可視化してみて、それっぽい形・性質のものを選ぶ
   • 左右対称のひと山ならとりあえず正規分布
   • 負の値を取らないならガンマ分布
   • 直線的か、指数関数的か、頭打ちか、などなど

客観的な指標もほしい。
モデルの尤もらしさといえば…

尤ゆう度 (likelihood)

あるモデル$M$の下でそのデータ$D$が観察される確率:
$\Pr(D \mid M)$

データ$D$を固定し、モデル$M$の関数とみなしたものが尤度関数:
$L(M \mid D)$

モデルの構造も固定してパラメータ$\theta$だけ動かす場合はこう書く:
$L(\theta \mid D)$ or $L(\theta)$

対数尤度 $\log L$ の形にしたほうがいろいろ便利。

各モデルで最適なパラメータを探して、比較:
$\log L^* (M_1) \text{ vs. } \log L^* (M_2) \text{ vs. } \log L^* (M_3) \ldots$

result = model.fit()
result.llf  # log likelihood

たしかに尤度はあてはまりの良さを表してそう

この場合は直線回帰よりもポアソン回帰が良さそう:

この調子で、より尤度の高いモデルを探していけばいいだろうか？

あてはまりが良ければいいってもんでもない

過剰適合 / 過学習 / overfitting

パラメータを増やせば現データへの適合度・尤度を高くできるが、
予測・理解の役には立たなくなる。

帰無モデル: 説明変数なし。切片のみ。
飽和モデル: データ点の数 ≤ パラメータの数。“データ読み上げ”的モデル

無駄な説明変数を加えても尤度は上がる

ある植物が作る種の数 $y$ は個体のサイズ $x$ に応じて増える。
観察時に着てた服の色 $x_2$ を追加すると尤度が上がる……?

AIC: 赤池情報量基準

\[\begin{split} \text{AIC} = -2 (\log L^* - k) = -2 \log L^* + 2k \end{split}\]

• AICが小さいほど予測精度の良いモデル。
  • 尤度は上げたい。
  • パラメータ数 $k$ が増えるとペナルティ。
• どのデータに対する当てはまりを目指すかという観点
  • 「手元のデータ」に対する対数尤度は $\log L^*$
    
  • 「真のメカニズムから出てくる未来のデータ」に対する
    平均対数尤度の推定量は $(\log L^* - k)$
    (Kullback–Leibler情報量を使って導出するらしい)

result = model.fit()
result.aic

無駄な説明変数の追加でAIC増加

ある植物が作る種の数 $y$ は個体のサイズ $x$ に応じて増える。
観察時に着てた服の色 $x_2$ を追加したモデルはAICが増加。

ほかの情報量基準

• $\text{BIC} = -2 \log L^* + k \log n$
  • パラメータ数 $k$ でペナルティを付けるのはAICと同じ。
  • データの観測数 $n$ に依存する点でAICと異なる。
    感覚としては「AICはデータサイズによるペナルティが無い」
  • (周辺尤度の最大化という観点で導出するらしい)
• WAIC, WBIC
  • AIC, BICを一般化し、広く使えるようにしたもの。
  • 理想的な条件ではそれぞれAIC, BICとほぼ同じ。
    そうじゃない場合(現実的には常に)こちらが優位。
  • WAICは予測の良さ、WBICは真のモデルへの近さ、を表す。

モデル選択の心構え

「正しい」ものを選べるわけではない。
予測・理解に useful なものを何らかの基準で選ぶだけ。

All models are wrong, but some are useful. — George E. P. Box

現実的な注意点・悩みどころ

• 多重共線性(multicollinearity):
  • 説明変数同士が強い相関関係にある
• 変数変換:
  • 気安くやるべきじゃないけど、対数変換などしばしば有用
  • 割り算した値は危険
• 交互作用を入れると解釈が難しくなる。

一般化線形モデル座学まとめ

• 何はともあれ散布図を描く
• 適切な確率分布・リンク関数・説明変数を考える
• パラメータを最尤推定する
• 尤度は「手元のデータへのあてはまり」
• モデルを比較するときは情報量基準を参考にする

☕️ 休憩 + 質疑応答

penguinsデータセット

https://allisonhorst.github.io/palmerpenguins/

import statsmodels.api as sm
penguins = sm.datasets.get_rdataset("penguins", "palmerpenguins", True).data
print(penguins)

penguinsデータセット

https://allisonhorst.github.io/palmerpenguins/

      species    island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g    sex year
  1    Adelie Torgersen           39.1          18.7               181        3750   male 2007
  2    Adelie Torgersen           39.5          17.4               186        3800 female 2007
  3    Adelie Torgersen           40.3          18.0               195        3250 female 2007
  4    Adelie Torgersen             NA            NA                NA          NA     NA 2007
 --                                                                                           
341 Chinstrap     Dream           43.5          18.1               202        3400 female 2009
342 Chinstrap     Dream           49.6          18.2               193        3775   male 2009
343 Chinstrap     Dream           50.8          19.0               210        4100   male 2009
344 Chinstrap     Dream           50.2          18.7               198        3775 female 2009

🔰 penguinsでGLMの練習

🔰 4-glm.ipynb に戻り、次の課題を解いてみよう。
(次ページ以降に解答。まずは自力で。)

1. body_mass_g を横軸、 flipper_length_mm を縦軸に、まず作図。
2. 単回帰して、切片と傾きを求める。そして作図。
3. species で色分けして作図。
4. species も説明変数に加えて重回帰し、切片と傾きを求める。そして作図。
5. 余裕があれば、クチバシの長さと深さを縦横軸にして同様の解析。

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

単回帰の練習: 1. まず作図

どうやら、重いペンギンほど翼長も長い。

grid = sns.FacetGrid(penguins)
grid.map(sns.scatterplot, "body_mass_g", "flipper_length_mm")

単回帰の練習: 2. モデル作成、フィッティング

とりあえずデフォルトの正規分布・恒等リンク。 $y = 136.7 + 0.0153 x$

formula = "flipper_length_mm ~ body_mass_g"
model1 = smf.glm(formula, data=penguins)
results1 = model1.fit()
print(results1.params)

Intercept      136.729559
body_mass_g      0.015276
dtype: float64

print(results1.llf)

-1145.5175473095946

print(results1.aic)

2295.0350946191893

単回帰の練習: 3. フィッティング結果を作図

結果とデータから予測値を作って回帰線を引く。

pen_pred = penguins.assign(pred=results1.predict(penguins))
grid = sns.FacetGrid(pen_pred)
grid.map(sns.scatterplot, "body_mass_g", "flipper_length_mm")
grid.map(sns.lineplot, "body_mass_g", "pred")

重回帰の練習: 1. まず作図

種によって色分けしてみると、傾向の違いが見える。

palette = {"Adelie": "#ff6600", "Gentoo": "#c35bcc", "Chinstrap": "#007174"}
grid = sns.FacetGrid(penguins, hue="species", palette=palette)
grid.map(sns.scatterplot, "body_mass_g", "flipper_length_mm")

重回帰の練習: 2. モデル作成、フィッティング

Adelieを基準に、ChinstrapとGentooはそれより長め。
体重の効果は単回帰のとき(0.0153)より小さい。

formula = "flipper_length_mm ~ body_mass_g + species"
model2 = smf.glm(formula, data=penguins)
results2 = model2.fit()
print(results2.params)

Intercept               158.860261
species[T.Chinstrap]      5.597440
species[T.Gentoo]        15.677470
body_mass_g               0.008402
dtype: float64

print(results2.llf)

-1059.7183131897373

print(results2.aic)

2127.4366263794745

重回帰の練習: 3. フィッティング結果を作図

pen_pred = penguins.assign(pred=results2.predict(penguins))
grid = sns.FacetGrid(pen_pred, hue="species", palette=palette)
grid.map(sns.scatterplot, "body_mass_g", "flipper_length_mm")
grid.map(sns.lineplot, "body_mass_g", "pred")

傾きも種によって違うかも。交互作用を入れてみたい。

交互作用の練習: モデル作成、フィッティング

Adelieを基準に、Chinstrapの傾きが結構違う。
切片の違いは解釈しにくくなった。

formula = "flipper_length_mm ~ body_mass_g + species + body_mass_g:species"
model3 = smf.glm(formula, data=penguins)
results3 = model3.fit()
print(results3.params)

Intercept                           165.244813
species[T.Chinstrap]                -13.863939
species[T.Gentoo]                     6.059376
body_mass_g                           0.006677
body_mass_g:species[T.Chinstrap]      0.005228
body_mass_g:species[T.Gentoo]         0.002362
dtype: float64

print(results3.llf)

-1055.7107640450004

print(results3.aic)

2123.4215280900007

交互作用の練習: フィッティング結果を作図

pen_pred = penguins.assign(pred=results3.predict(penguins))
grid = sns.FacetGrid(pen_pred, hue="species", palette=palette)
grid.map(sns.scatterplot, "body_mass_g", "flipper_length_mm")
grid.map(sns.lineplot, "body_mass_g", "pred")

ここまでの3つのモデルでどれがいいか？

AICで選ぶなら交互作用入り重回帰が良さそう。

results1.aic
results2.aic
results3.aic

余裕があったら追加の練習

🔰クチバシの長さと深さで同じ解析をやってみよう。

🔰 手元のデータ、公共データなどでGLMしてみよう

正規分布・恒等リンクじゃないものだとなお良し。

Pythonパッケージから読み込めるものを探すのもあり。
e.g., sm.datasets.get_rdataset(item, package)

import seaborn as sns
sns.get_dataset_names()
titanic = sns.load_dataset("titanic")

import statsmodels.api as sm
iris = sm.datasets.get_rdataset("iris").data
diamonds = sm.datasets.get_rdataset("diamonds", "ggplot2").data

一般化線形モデル(GLM)まとめ

• 何はともあれ作図して俯瞰
• GLMは統計モデリングの考え方の根幹
  • 確率分布・リンク関数・説明変数
  • 尤度・最尤法によるパラメータ推定
  • 情報量基準などによるモデル選択

参考文献

• データ解析のための統計モデリング入門 久保拓弥 2012
• StanとRでベイズ統計モデリング 松浦健太郎 2016
• RとStanではじめる ベイズ統計モデリングによるデータ分析入門 馬場真哉 2019
• データ分析のための数理モデル入門 江崎貴裕 2020
• 分析者のためのデータ解釈学入門 江崎貴裕 2020
• 統計学を哲学する 大塚淳 2020