ロジスティック回帰分析をやってみた(1)

自分のためのメモです。

#ロジスティック回帰分析
まず初めに通常のロジスティック回帰分析の説明をします。
##オッズ

通常というのはアウトカム（見つけたい結果）が2値の場合だ。例えば、「１：病気にかかった」、「０：病気にかからなかった」のようなものだ。
ここでオッズとはなんなのかを考える。
語源はギャンブルから来たようである。
オッズというのは競馬のオッズで、「その馬券が的中した場合、何倍の配当がつくのか」を「倍率」として表したものである。
ある事象の起こる確率 p と起こらない確率 1 − p との比 p / (1 − p) のことをオッズと定義した。
ある事象の起こる確率 p が 1/2 を超えることはオッズが 1 を超えることに等しく、$0 ≤ p < 1$ の範囲で確率とオッズは1対1に対応し、確率とオッズは同じものの別表現になっている。
例えば、オッズが5.0倍の単勝を100円分買ったとして、その馬券が的中したら

5.0×100円=500円

の払戻し金を受け取れるのだ。外れたら0円である。

```math
odds=\frac{p}{1-p}
```
odds=5.0なら、p/(1-p)=5.0なので


```math
\begin{align}
5.0&=\frac{p}{1-p}\\
p&=5.0\times (1-p)\\
p&=5.0-5.0p\\
6p&=5.0\\
p&=\frac{5.0}{6}
\end{align}
```

または、「１：病気にかかる」率を20%とする。

「０：病気にかからなかった」率は１−ｐで80％である。そこでこのオッズは
$odds=\frac{0.2}{1-0.2}=\frac{0.2}{0.8}=0.25$
となります。
さて、オッズ比（OR）は以下のように定義されます。

```math
\begin{align}
odds\_ ratio(OR) &= \frac{\frac{p_1}{(1-p_1)}}{\frac{p_0}{(1-p_0)}}\\
 & = \frac{p_0 (1-p_1)}{p_1 (1-p_0)}
\end{align}
```
なぜ、オッズ比を説明するかはこれから説明します。
ロジスティック回帰は係数のオッズ比を求めることになります。
係数が定数項$\beta_0$と変数$x$に対する$\beta_1$の場合で説明します。

```math
\begin{align}
p &=\frac{exp(\beta_0+\beta_1x)}{1+exp(\beta_0+\beta_1x)}\\
odds &=\frac{p}{1-p}\\
log\left(\frac{p}{1-p}\right) &=\beta_0+\beta_1x
\end{align}
```

ここで$x=1$のときと、$x=0$に分けて考えてみます。


```math
x=0の場合\\
p=\frac{exp(\beta_0)}{1+exp(\beta_0)}\\
L1=log\left(\frac{p_0}{1-p_0}\right)=\beta_0
```

とL1は$\beta_0$が対数オッズになります。対数オッズのことをlogitと呼びます。

```math
x=1の場合\\
p=\frac{exp(\beta_0+\beta_1)}{1+exp(\beta_0+\beta_1)}\\
L2=log\left(\frac{p_1}{1-p_1}\right)=\beta_0+\beta_1
```

そこで$L2-L1$を考えると、$\beta_1$のみ残ります。
これはオッズ比になっています。

```math
\begin{align}
L2-L1&=log\left(\frac{p_1}{1-p_1}\right)-log\left(\frac{p_0}{1-p_0}\right)\\
\frac{log\left(\frac{p_1}{1-p_1}\right)}{log\left(\frac{p_0}{1-p_0}\right)}&=L2-L1\\
log\left(\frac{\frac{p_1}{1-p_1}}{\frac{p_0}{1-p_0}}\right)&=\beta_0+\beta_1-\beta_0\\
&=\beta_1
\end{align}
```
すなわち、$x$の係数$\beta_1$はオッズ比になっています。
ロジスティック回帰はオッズ比を求めることになりますが、変数が多いときには内部の共分散で調整した値を計算します。

統計ソフトRではglmという関数を使います。

例として簡単なデータを示します。

```R
head(ldata)
  admit gre  gpa rank rankx
1     0 380 3.61    3     0
2     1 660 3.67    3     0
3     1 800 4.00    1     1
4     1 640 3.19    4     0
5     0 520 2.93    4     0
6     1 760 3.00    2     1
> 
## 
tab1<-xtabs(~admit+rankx, data=ldata)
> tab1
     rankx
admit   0   1
    0 148 125
    1  40  87
```

1項目だけのロジスティック回帰を行います。

```R
out1<-glm(admit~rankx,family="binomial",data=ldata)
> summary(out1)
```
結果は以下のようになります。

```R
>summary(out1)

Call:
glm(formula = admit ~ rankx, family = "binomial", data = ldata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0279  -1.0279  -0.6917   1.3347   1.7593  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3083     0.1782  -7.342 2.11e-13 ***
rankx         0.9459     0.2264   4.178 2.94e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 481.66  on 398  degrees of freedom
AIC: 485.66

Number of Fisher Scoring iterations: 4

```
Estimatesは対数のオッズ比$(log(\beta_0), log(\beta_1))$ですので、指数倍します。

```R
OR<-exp(out1$coefficients)
> OR
(Intercept)       rankx 
  0.2702703   2.5752000 

```
次にオッズ比の95%信頼区間を求めます。

```R
CI<-confint(out1)
> CI
                 2.5 %     97.5 %
(Intercept) -1.6703179 -0.9697971
rankx        0.5078528  1.3968586
> 
```
95%CIはまだ、対数なので、指数化して戻します。

```R
exp(CI)
                2.5 %    97.5 %
(Intercept) 0.1881872 0.3791599
rankx       1.6617193 4.0424808
```

以上、次の回ではもう少し、掘り下げます。