みどり本第3章の例題続きです。
前回は例題のデータを分析して、以下のような特徴が得られたのでした。今回はこういった特徴を考慮しながら統計モデルを作っていきます。
- 種子数
:
- 非負の整数（離散値）
- 区間[2, 15]
- 平均と分散がまあまあ近い
- 体サイズ
:
- 正の実数（連続値）
- 区間[7.19, 12.4]
- 分散は小さい
- 施肥処理の有無
(Tが施肥処理有り):
- 因子（カテゴリ）
- 各カテゴリの要素数は
- （散布図とヒートマップより）サイズ
- （箱ひげ図より）肥料の効果は無さそう
ポアソン回帰で統計モデリング
個体ごとの平均種子数が体サイズ
や施肥処理
に影響されるようなモデルを設計。
モデルを当てはめる基本の流れは第2章と同じ。
- 標本が従うと考えられる真の確率分布（母集団分布）を何かしら仮定
- 尤度関数（仮定した確率分布の確率質量(または密度)関数
を、全ての標本について総積をとったもの）を最大にするパラメータを求める
1. 種子数が個体の体サイズのみに依存する統計モデル
説明変数が数量のみの統計モデル
（施肥処理の有無関係無く、すべての標本のを用いる。）
1. ある個体の種子数がである確率はポアソン分布に従うと仮定
$$
y_i \sim Po(\lambda_i)
$$
$$
p(y_i\mid\lambda_i)=\frac{\lambda_i^{y_i}exp(-\lambda_i)}{y_i !}
$$
式の通りですが、個体の種子数
はポアソン分布に従うと仮定するものの、その分布の形を決めるパラメータ
が個体ごとに異なるので、ここでは個体ごとに異なる形をしたポアソン分布から種子数が発生すると考えています。私は最初ここで躓いたので、念の為…。
今回は、個体ごとに異なる平均
を説明変数
の関数として定義する必要がある。以下のように定義する。
$$
\lambda_i=exp(\beta_1+\beta_2 x_i)
$$
両辺対数をとると
$$
log(\lambda_i)=\beta_1+\beta_2 x_i
$$
(の関数)=(線形予測子)となっているとき、左辺の関数をリンク関数（今回は対数関数なので、対数リンク関数関数）と呼ぶ。
ポアソン回帰の一般化線形モデル（GLM）で対数リンク関数を使う理由：ポアソン分布の平均は非負である必要があるから、上式だとexp(線形予測子)>=0となり都合が良い。また、要因の効果が積で表せるので分かりやすい。
2. 対数尤度関数を最大化するパラメータ値を推定
$$
\log L(\beta_1,\beta_2)=\sum_i\log\frac{\lambda_i^{y_i}exp(-\lambda_i)}{y_i !}
$$
今回の場合は、、
の推定値を求める
複数のパラメータの対数尤度関数を解析的に最適化するのは困難だが、プログラム上では数値的な試行錯誤によって求まるので問題ない。
# http://www.statsmodels.org/devel/examples/notebooks/generated/glm_formula.html import statsmodels.formula.api as smf fit_x = smf.poisson('y ~ x', data=df).fit() fit_x.summary()
Optimization terminated successfully.
Current function value: 2.353863
Iterations 4
| coef | std err | |
|---|---|---|
| Intercept | 1.2917 | 0.364 |
| x | 0.0757 | 0.036 |
表中のcoefの列が各パラメータの最尤推定値。
最尤推定値をモデルに代入し、平均種子数を予測
$$
\lambda = exp(1.2917+0.0757x)
$$
xに対応する予測平均値を描画。
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)
サイズの増加に応じて平均種子数が増加する、観測データに合った"それっぽい"線が引けました。
さて、これで種子数が個体のサイズのみに依存する統計モデルが作れました。次回は、施肥処理の効果に依存するモデルを考えていきます。
コメント