確率変数$X$はパラメータ$\theta$の確率分布$f$に“従う”: $X \sim f(\theta) $
e.g., ある植物が作る種の数$X$は平均値$\lambda$のポアソン分布に従う:
これを一般化線形モデル(GLM)として見ることもできる。
個体$i$の種子数$y_i$は平均値$\lambda_i$のポアソン分布に従う。
平均値$\lambda_i$は他のデータによらず$\beta_0$で一定。
種子数をY軸にして、式を2つに分けただけ…?
説明変数を含むモデルを見ればご利益が分かるかも。
個体$i$の種子数$y_i$は平均値$\lambda_i$のポアソン分布に従う。
平均値の対数$\log(\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}\]
気温も湿度も高いほどビールが売れる、とか
今度は確率分布とリンク関数を変えてみよう。
何かの成否に対する何かの因子の影響、とか
客10人中$y_i$人がビールを注文。
その日$i$の気温$x_i$によって割合が変化。
\[\begin{split} y_i &\sim \text{Binomial}(n,~p_i) \\ \text{logit}(p_i) &= \beta_0 + \beta_1 x_i \\ p_i &= \frac 1 {1 + e^{-(\beta_0 + \beta_1 x_i)}} \end{split}\]
ロジスティック関数↑
何かの成否に対する何かの因子の影響、とか
風が吹けば桶屋が儲かる。
\[\begin{split} y_i &\sim \text{Bernoulli}(p_i) \\ &= \text{Binomial}(1,~p_i) \\ \text{logit}(p_i) &= \beta_0 + \beta_1 x_i \\ p_i &= \frac 1 {1 + e^{-(\beta_0 + \beta_1 x_i)}} \end{split}\]
ロジスティック関数↑
\[\begin{split} y_i &\sim \mathcal{N}(\mu_i,~\sigma^2) \\ \text{identity}(\mu_i) &= \beta_0 + \beta_1 x_i \end{split}\]
最小二乗法の直線あてはめと結果的に同じになる。
単回帰・重回帰と言ったとき一般線形モデルを前提とする人もいる。
質的な説明変数を持つ正規分布・恒等リンクのGLM、と解釈可能。
指示変数 (0 or 1) に変換してから重回帰する。
天気 | → | $x_1$ ☀️ 晴れ | $x_2$ ☔️ 雨 |
---|---|---|---|
☁️ くもり | 0 | 0 | |
☀️ 晴れ | 1 | 0 | |
☔️ 雨 | 0 | 1 |
\[\begin{split} y_i &= \mathcal{N}(\mu_i,\sigma^2) \\ \mu_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} \end{split}\]
くもり☁️ $\beta_0$ を基準に、晴れの効果☀️ $\beta_1$ と雨の効果☔️ $\beta_2$ が求まる。
GLMなら確率分布・リンク関数を変えてもっと柔軟にモデリングできる。
質的変数と量的変数を両方含むGLM、と解釈可能。
正規分布・等分散・恒等リンクなどが仮定される。
天気 | → | $x_1$ ☀️ 晴れ | $x_2$ ☔️ 雨 |
---|---|---|---|
☁️ くもり | 0 | 0 | |
☀️ 晴れ | 1 | 0 | |
☔️ 雨 | 0 | 1 |
\[\begin{split} y_i &= \mathcal{N}(\mu_i,\sigma^2) \\ \mu_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} \end{split}\]
GLMなら確率分布・リンク関数を変えてもっと柔軟にモデリングできる。
確率分布・リンク関数を変えて柔軟にモデリングできる。
特定の組み合わせには名前がある。
名前 | 確率分布 | リンク関数 | 説明変数 |
---|---|---|---|
ポアソン回帰 | ポアソン分布 | log | |
ロジスティック回帰 | 二項分布 | logit | |
一般線形回帰 | 正規分布 | 恒等 | |
分散分析 | 正規分布 | 恒等 | 質的変数 |
共分散分析 | 正規分布 | 恒等 | 質的変数+量的変数 |
確率分布については前章を参照。
リンク関数をもう少しだけ掘り下げたい。
統計モデリングにおいて「まっすぐ以外も表現できる」意味
ほかに probit
, inverse
, sqrt
, etc.
どう選ぶ?
客観的な指標もほしい。
モデルの尤もらしさといえば…
あるモデル$M$の下でそのデータ$D$が観察される確率:
$\text{Prob}(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$
この場合は直線回帰よりもポアソン回帰が良さそう:
この調子で、より尤度の高いモデルを探していけばいいだろうか?
帰無モデル: 説明変数なし。切片のみ。
飽和モデル: データ点の数 ≤ パラメータの数。“データ読み上げ”的モデル
ある植物が作る種の数 $y$ は個体のサイズ $x$ に応じて増える。
観察時に着てた服の色 $x_2$ を追加すると尤度が上がる……?
\[\begin{split} \text{AIC} = -2 (\log L^* - k) = -2 \log L^* + 2k \end{split}\]
ある植物が作る種の数 $y$ は個体のサイズ $x$ に応じて増える。
観察時に着てた服の色 $x_2$ を追加したモデルはAICが増加。
「正しい」ものを選べるわけではない。
予測・理解に useful なものを何らかの基準で選ぶだけ。
All models are wrong, but some are useful. — George E. P. Box
ある説明変数の効果が、別の説明変数によって異なる。
e.g., ビール売上の温度依存性が天気によって異なる。
天気 | $x_1$ |
---|---|
☀️ 晴れ | 1 |
☔️ 雨 | 0 |
\[\begin{split} y_i &= \mathcal{N}(\mu_i,\sigma^2) \\ \mu_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_{1,2} x_{1i} x_{2i} \end{split}\]
雨の日は $x_{1i} = 0$ のため $\beta_0,~\beta_2$ の項だけ。
晴れの日はそれに加えて $\beta_1,~\beta_{1,2}$ の項も。
解釈が一気に難しくなるのでむやみに使わない。
以降のコードで共通のライブラリ読み込み
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
palette = {"Adelie": "#ff6600", "Gentoo": "#c35bcc", "Chinstrap": "#007174"}
# Python
import statsmodels.api as sm
penguins = sm.datasets.get_rdataset('penguins', 'palmerpenguins').data
print(penguins)
# R
library(palmerpenguins)
print(penguins)
species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex year
<fct> <fct> <dbl> <dbl> <int> <int> <fct> <int>
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
どうやら、重いペンギンほど翼長も長い。
sns.relplot(x='body_mass_g', y='flipper_length_mm', data=penguins)
結構たくさん出力される。上がモデルの評価、下が推定結果。
formula = 'flipper_length_mm ~ body_mass_g'
model1 = smf.glm(formula, data=penguins)
results1 = model1.fit()
print(results1.summary())
Dep. Variable: flipper_length_mm No. Observations: 342
Model: GLM Df Residuals: 340
Model Family: Gaussian Df Model: 1
Link Function: identity Scale: 47.795
Method: IRLS Log-Likelihood: -1145.5
Date: Tue, 29 Jun 2021 Deviance: 16250.
Time: 10:32:40 Pearson chi2: 1.63e+04
No. Iterations: 3
Covariance Type: nonrobust
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 136.7296 1.997 68.473 0.000 132.816 140.643
body_mass_g 0.0153 0.000 32.722 0.000 0.014 0.016
$y = 136.7 + 0.0153 x$
y_pred = results1.predict(penguins)
grid = sns.relplot(x='body_mass_g', y='flipper_length_mm', data=penguins)
grid.map(sns.lineplot, x=penguins['body_mass_g'], y=y_pred)
重いペンギンほど翼長も長い。翼長は種によっても違うかも。
sns.relplot(x='body_mass_g', y='flipper_length_mm', hue='species', data=penguins, palette=palette)
Adelieを基準に、ChinstrapとGentooはそれより長め。
体重の効果は単回帰のときより小さい。
formula = 'flipper_length_mm ~ body_mass_g + species'
model2 = smf.glm(formula, data=penguins)
results2 = model2.fit()
print(results2.summary())
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------------
Intercept 158.8603 2.387 66.564 0.000 154.183 163.538
species[T.Chinstrap] 5.5974 0.788 7.101 0.000 4.053 7.142
species[T.Gentoo] 15.6775 1.091 14.374 0.000 13.540 17.815
body_mass_g 0.0084 0.001 13.255 0.000 0.007 0.010
pen_pred = penguins.assign(pred=results2.predict(penguins))
grid = sns.relplot(x='body_mass_g', y='flipper_length_mm', hue='species', data=pen_pred, palette=palette)
grid.map(sns.lineplot, x='body_mass_g', y='pred', hue='species', data=pen_pred, palette=palette)
傾きも種によって違うかも。交互作用を入れてみたい。
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.summary())
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------------------------
Intercept 165.2448 3.551 46.536 0.000 158.285 172.204
species[T.Chinstrap] -13.8639 7.301 -1.899 0.058 -28.174 0.446
species[T.Gentoo] 6.0594 6.051 1.001 0.317 -5.800 17.919
body_mass_g 0.0067 0.001 7.011 0.000 0.005 0.009
body_mass_g:species[T.Chinstrap] 0.0052 0.002 2.683 0.007 0.001 0.009
body_mass_g:species[T.Gentoo] 0.0024 0.001 1.746 0.081 -0.000 0.005
pen_pred = penguins.assign(pred=results3.predict(penguins))
grid = sns.relplot(x='body_mass_g', y='flipper_length_mm', hue='species', data=pen_pred, palette=palette)
grid.map(sns.lineplot, x='body_mass_g', y='pred', hue='species', data=pen_pred, palette=palette)
AICで選ぶなら交互作用入り重回帰のが良さそう。
results1.aic
results2.aic
results3.aic
🔰クチバシの長さと深さで同じ解析をやってみよう。
🔰余裕があったら性別や年なども説明変数に入れてみよう。
何も指定しない場合は正規分布・恒等リンクだった:
formula = 'flipper_length_mm ~ body_mass_g'
model = smf.glm(formula, data=penguins)
#### results.summary()
# Model Family: Gaussian
# Link Function: identity
こう書いたのと同じ:
formula = 'flipper_length_mm ~ body_mass_g'
link = sm.families.links.identity
family = sm.families.Gaussian(link=link)
model = smf.glm(formula, data=penguins, family=family)
正規分布・恒等リンクじゃないものだとなお良し。
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
Slackで報告していただけると嬉しいです。
植物100個体から8個ずつ種子を取って植えたら全体で半分ちょい発芽。
親1個体あたりの生存数はn=8の二項分布になるはずだけど、
極端な値(全部死亡、全部生存)が多かった。個体差?
各個体の生存率$p_i$をそのままパラメータにすると過剰適合。
「パラメータ数 ≥ サンプルサイズ」の“データ読み上げ”モデル。
i.e., この個体は4個生き残って生存率0.5だね。次の個体は2個体だから……
個体の生存能力をもっと少ないパラメータで表現できないか?
各個体の生存率$p_i$が能力値$z_i$のシグモイド関数で決まると仮定。
その能力値は全個体共通の正規分布に従うと仮定:
$z_i \sim \mathcal{N}(\hat z, \sigma)$
パラメータ2つで済む: 平均 $\hat z$, ばらつき $\sigma$ 。
前者は標本平均 $\hat p$ から求まるとして、後者どうする?
普通の二項分布は個体差無し $\sigma = 0$ を仮定してるのと同じ。
正規分布と二項分布の混ぜ合わせ……?
パラメータp(を決めるz)ごとに二項分布を作って、重み付けして足したもの。
固定効果(fixed effects) のみ扱っていたGLMを拡張して、
変量効果(random effect) を混合したモデル。
「混合分布を使うモデル」という意味ではないらしい。
\[\begin{split} y_i &\sim \text{Binomial}(n,~p_i) \\ \text{logit}(p_i) &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \ldots + z_{1i} + \ldots \\ z_{1i} &\sim \mathcal{N}(\mu_1,~\sigma_1) \end{split}\]
e.g.,
個体$i$の種子生存率$p_i$は、
(固定効果) 体サイズ$x_{1i}$と日当たり$x_{2i}$に依存し、
(変量効果) よくわからん個体差$z_{1i}$と植木鉢差$z_{2i}$もある。
推定したパラメータを予測に使うなら固定効果
データに擬似反復が含まれるとき。
ぜんぶ独立のつもりで解析すると推定が偏ったり誤ったり。
植木鉢 | 個体/植木鉢 | 種子/個体 | 疑似反復 | 推定不可 |
---|---|---|---|---|
100個 | 1個体ずつ | 1個ずつ | – | 個体差・鉢差 |
25個 | 1個体ずつ | 4個ずつ | 個体 | 鉢差 |
20個 | 5個体ずつ | 1個ずつ | 植木鉢 | 個体差 |
5個 | 5個体ずつ | 4個ずつ | 植木鉢・個体 | – |
疑似反復あり
→ 観測できなかった個体差・場所差(変量効果)を推定可能
→ そのぶんを差し引いて固定効果を推定したい
→ ここでGLMMの練習はせず、階層ベイズモデルに進む。