はじめに
今回もロジスティック回帰(logistic regression)を扱います。
ただし、少し変わっているので(個人的にはどうやんだろう)ってなった問題です。
yが二項分布にしたがっているときのロジスティック回帰分析です。
普段は2値のロジスティック回帰だと思いますが、二項分布にしたがっているので分類という感じではなく、回帰のような使い方をしています。多クラス分類の話でもないです。
一般的なロジスティック回帰はこちら
https://zenn.dev/totopironote/articles/48e39f9b7fef2b
設定
用量反応実験
マウスがn=10匹の群が11個あり、薬の投与量xが少ない順に1,2,…,11と群に番号がついている。
y_i = 群iのマウスの死亡数 , 死亡率 p_i= y_i /10 とする。
| x |
y |
p |
| 1 |
0 |
0.0 |
| 2 |
0 |
0.0 |
| 3 |
0 |
0.0 |
| 4 |
3 |
0.3 |
| 5 |
6 |
0.6 |
| 6 |
6 |
0.6 |
| 7 |
5 |
0.5 |
| 8 |
9 |
0.9 |
| 9 |
9 |
0.9 |
| 10 |
10 |
1.0 |
| 11 |
10 |
1.0 |
yを二項分布でモデル化すると
y_i \hspace{2pt} i.i.d. \hspace{2pt}Bi(n_i,\pi_i) , i=1,2,...,N
N=11であり、n_iはすべて10である。
より\pi_i はp_iによって不偏推定される。
ロジットパラメータ\lambdaが以下の一次関数でかけるとする。
\lambda_i = \log\left \{\frac{\pi_i}{1-\pi_i}\right\}=\beta_0+\beta_1x_i
xは投与量である。
\piについて書き直すと
\pi_i = (1+e^{-\lambda})^{-1}
最尤推定
最尤推定量(\hat{\beta_0},\hat{\beta_1})を求める。
まず、二項分布の尤度関数Lは
L = \Pi_{i=1}^N\begin{pmatrix}n_i\\
y_i\end{pmatrix}\pi^y_i(1-\pi_i)^{n_i-y_i}
マイナス対数尤度lは
\begin{align*}
l = -\log L &\propto -\sum_{i=1}^N\{y_i\log \pi_i + (n_i-y_i)\log (1-\pi_i)\} \\
&= -\sum_{i=1}^N\left\{y_i\log\left (\frac{\pi_i}{1-\pi_i}\right) +n_i \log(1-\pi_i)\right\}\\
&= \sum_{i=1}^N \{-y_i(\beta_0+\beta_1x_i)+n_i\log (1+e^{\beta_0+\beta_1x_i})\}
\end{align*}
\beta_0,\beta_1で1回微分すると
\begin{align*}
\frac{\partial l}{\partial \beta_0}
&=\sum_{i=1}^N \{\frac{n_i\exp(\beta_0+\beta_1x_i)}{1+\exp(\beta_0+\beta_1x_i)}-y_i\}\\
&=\sum_{i=1}^N(n_i\pi_i-y_i)=0
\end{align*}
\begin{align*}
\frac{\partial l}{\partial \beta_1}
&=\sum_{i=1}^N \{\frac{n_ix_i\exp(\beta_0+\beta_1x_i)}{1+\exp(\beta_0+\beta_1x_i)}-x_iy_i\}\\
&=\sum_{i=1}^N\{x_i(n_i\pi_i-y_i)\}=0
\end{align*}
\pi_iの微分を考えると
\frac{\partial \pi_i}{\partial \beta_0}
= \pi_i(1-\pi_i) , \frac{\partial \pi_i}{\partial \beta_1}
= x_i\pi_i(1-\pi_i)
よって2回微分はx_0=1としている。j,k=0,1に対して
\frac{\partial^2 l}{\partial \beta_j\partial \beta_k}=\sum_{i=1}^Nx_{ij}x_{ik}n_i\pi_i(1-\pi_i)
以上の計算を綺麗にベクトル表記に直すと
\utilde{\beta}=\begin{pmatrix}
\beta_0\\
\beta_1
\end{pmatrix}, X=\begin{pmatrix}
1 &x_{1} \\
1 &x_{2}\\
\vdots & \vdots \\
1 &x_{N}
\end{pmatrix},\utilde{y}=\begin{pmatrix}
y_1\\
\vdots\\
y_N
\end{pmatrix}
\utilde{\pi} = \begin{pmatrix}
\pi_1\\
\vdots\\
\pi_N
\end{pmatrix}, \log\frac{\utilde{\pi}}{1_N-\utilde{\pi}}=\begin{pmatrix}
\log\frac{\pi_1}{1-\pi_1}\\
\vdots\\
\log\frac{\pi_N}{1-\pi_N}
\end{pmatrix}=X\utilde{\beta},\utilde{\mu}=\begin{pmatrix}
n_1\pi_1\\
\vdots\\
n_N\pi_N
\end{pmatrix}
これらを用いると、
マイナス対数尤度は
l = \log(1+\exp(X\utilde{\beta}))-\utilde{y}^TX\utilde{\beta}
1回微分は
\dot{l} = X^T(\utilde{\mu}-\utilde{y})=0
2回微分は
\ddot{l} = X^TVX,
V=diag(v_1,\cdots,v_N),v_i=n_i\pi_i(1-\pi_i)
よってニュートンラフソン法より
\begin{align*}
\utilde{\beta}^{new}&=\utilde{\beta}^{old}-\ddot{l}^{-1}\dot{l}\\
&=\utilde{\beta}^{old}-(X^TVX)^{-1}X^T(\utilde{\mu}-\utilde{y})
\end{align*}
回帰曲線の図示
まず、回帰曲線は最尤推定量(\hat{\beta_0},\hat{\beta_1})を代入したものである。
\hat{\pi}(x)=(1+e^{-(\hat{\beta}_0+\hat{\beta}_1x)})^{-1}
pythonでsklearnを用いて図示する方法を説明する。
2値のロジスティック回帰の文脈に乗っけて使う。
そのためpの値が0.5より大きいものを1,0.5以下を0としている。
x = range(1, 12)
y = [0, 0, 0, 3, 6, 6, 5, 9, 9, 10, 10]
y_ = [0 if i/10 <= 0.5 else 1 for i in y]
from sklearn.linear_model import LogisticRegression
model = LogisticRegression()
X = [[i] for i in x]
model.fit(X, y_)
new_x = np.linspace(0, 12, 100)
new_X = np.array([[i] for i in new_x])
predicted_probabilities = model.predict_proba(new_X)
predicted_p = predicted_probabilities[:, 1]
p =[i/10 for i in y]
plt.scatter(x,p)
plt.xlabel('投与量')
plt.ylabel('死亡率')
plt.plot(new_x, predicted_p, color='r', label='Predicted')
plt.legend()
plt.show()
標準誤差
単純な標準誤差を求める。
sd = [\hat{\pi}(x)(1-\hat{\pi}(x))/10]^{1/2}
| x |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
| sd |
0.058424 |
0.081335 |
0.109124 |
0.137118 |
0.155646 |
0.155277 |
0.136224 |
0.108096 |
0.08043 |
0.057721 |
0.040653 |
まとめ
教科書と若干数値がずれたのが気になるけど、ほぼ2値のロジスティック回帰と同じような感じになった。
参考文献
- Bradley Efron,Trevor Hastie 「大規模計算時代の統計推論」
Discussion