Rでロジスティック回帰をglm()なしでやってみる

# はじめに
ロジスティック回帰をツールを使ってやる機会があったので勉強用に
あえてRの`glm()`をつかわずにやってみようと思い記事にしました．

## 実行環境
今回の実行環境です．

```{r}
R.Version()$version.string
```

```{r:Console}
[1] "R version 3.6.2 (2019-12-12)"
```

ライブラリは次のものを使いました．

```{r}
library(tidyverse)
```

# ロジスティック回帰とは
次のような直線があるとします．

```math
y = ax+b
```

ここでデータ，$(x_i,y_i)$を示す点とこの直線の誤差 $e_i$ を次のように定義します.

```math
e_i = y_i-(ax_i+b)
```
この誤差の二乗和を最小にする$a$,$b$を求めた時，直線 $y = ax+b$を，$(x_i,y_i)$ の $y$ の $x$ への回帰直線と呼びます．
しかし，応答変数が不良品率のような確率をとる時この直線のまま考えるのは適切でない場合があります．
（説明変数を大きく/小さくすると応答変数が0または1を超えてしまうため．）

そこで，次のようなロジットモデルが用いられます．

```math
L(p) = ln\frac{p}{1-p}=ax+b
```
上の式を $p$ について解くと，

```math
p = \frac{1}{1+exp(-ax-b)}
```
となり$x$の値に関わらず， $p$ が $[0,1]$ からはみ出すことがなくなります．
Rで次のように関数を定義し図示して確認してみます．

```{r}
 logistic<-function(x)1/(1+exp(-x))
 ggplot() + 
  geom_line(data=data.frame(x=seq(-5,5,by=0.1),y=logistic(seq(-5,5,by=0.1))),
            aes(x=x,y=y),size=1.2) +
  theme(axis.title = element_text(size = 12),
        axis.text = element_text(size = 12)) +
  geom_hline(yintercept = 0,lty=2) +
  geom_hline(yintercept = 1,lty=2)
```
 ![Rplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/532006/ea06e926-ff85-ca99-7994-69543e57513a.png)

ロジスティック回帰では，応答変数が $p$ という確率で1となり，確率 $1-p$ で0になる二項分布について考えます．
パラメータの推定には最尤法を用います．
尤度関数 $L(\theta)$ は次のようになります．

```math
L(\theta|x) = \prod_{i=1}^n \left(\frac{1}{1+exp(-ax_i-b)}\right)^{y_i}\left(1-\frac{1}{1+exp(-ax_i-b)}\right)^{1-y{_i}}
```
推定の際は，$L(\theta)$ を最大にする$\theta$と $lnL(\theta)$ を最大にする $\theta$ も同じなので
対数をとった対数尤度関数を最大化します．

```math
lnL(\theta|x) = \sum_{i=1}^{n} y_iln(\frac{1}{1+exp(-ax_i-b)})+(1-y_i)ln(1-\frac{1}{1+exp(-ax_i-b)})
```
実際にロジスティック回帰をやってみます．

## データ
サンプル用のデータとして，有名な`titanic`を用います．

```{r}
url <- "http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic.txt"
titanic <- read_csv(file = url)
```

変数として年齢を使おうと思ったのですが，欠損値が多いので
補います．ここでは単純に，利用できるデータの平均をとって代入しています．

```{r}
titanic$age <- ifelse(is.na(titanic$age),mean(titanic$age,na.rm = T),titanic$age)
```
まずはRの`glm()`でロジスティック回帰を実行します．
説明変数にはクラス，性別，年齢を使ってみます．

```{r}
res <- glm(survived ~ pclass + sex + age,
                     family = binomial, data = titanic)
summary(res)
```

```{r:Console}
Call:
glm(formula = survived ~ pclass + sex + age, family = binomial, 
    data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6452  -0.6641  -0.3679   0.6123   2.5615  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.552261   0.342188  10.381  < 2e-16 ***
pclass2nd   -1.170777   0.211559  -5.534 3.13e-08 ***
pclass3rd   -2.430672   0.195157 -12.455  < 2e-16 ***
sexmale     -2.463377   0.154587 -15.935  < 2e-16 ***
age         -0.042235   0.007415  -5.696 1.23e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1686.8  on 1312  degrees of freedom
Residual deviance: 1165.7  on 1308  degrees of freedom
AIC: 1175.7

Number of Fisher Scoring iterations: 5
```

今回は，この推定値に自力でたどり着けるようにします．

まずはデータの前処理を行います．
カテゴリー変数を`1`，`0`の数値に置き換えます．

```{r}
Data <- titanic %>% 
  mutate(.,pclass_2nd = ifelse(pclass == "2nd" ,1,0),
           pclass_3rd = ifelse(pclass == "3rd",1,0),
           sex_male = ifelse(sex == "male" ,1,0)) %>% 
  select(.,survived,pclass_2nd,pclass_3rd,sex_male,age)
head(Data,n=5)
```

```{r:Console}
# A tibble: 5 x 5
  survived pclass_2nd pclass_3rd sex_male    age
     <dbl>      <dbl>      <dbl>    <dbl>  <dbl>
1        1          0          0        0 29    
2        0          0          0        0  2    
3        0          0          0        1 30    
4        0          0          0        0 25    
5        1          0          0        1  0.917
```

負の対数尤度関数を定義します．

```{r:Negative_Log_likelihood_function.R}

logistic_reg<-function(x1,x2,x3,x4,b0,b1,b2,b3,b4){
  a<-b0+x1*b1+x2*b2+x3*b3+x4*b4
  logistic(x=a)
}
 
nll_logistic<-function(para,x1,x2,x3,x4,y){
  b0<-para[1]
  b1<-para[2]
  b2<-para[3]
  b3<-para[4]
  b4<-para[5]
  
  tmp<-log(ifelse(y==1,logistic_reg(x1=x1,x2=x2,x3=x3,x4=x4,b0=b0,b1=b1,b2=b2,b3=b3,b4=b4),
                  1-logistic_reg(x1=x1,x2=x2,x3=x3,x4=x4,b0=b0,b1=b1,b2=b2,b3=b3,b4=b4)))
  -1*sum(tmp)
}
```
パラメータ推定にはRの`optim`関数を使います．

```{r}
inits<-res$coefficients
res2<-optim(inits,nll_logistic,x1=Data$pclass_2nd,x2=Data$pclass_3rd,
            x3=Data$sex_male,x4=Data$age, y=Data$survived,
            method = "BFGS",hessian = T) 
res2$par  
```
```{r:Console}
(Intercept)   pclass2nd   pclass3rd     sexmale         age 
 3.55226110 -1.17077691 -2.43067246 -2.46337730 -0.04223535 
```
同じ推定値になりました．
AICも確認します．

```{r:AIC}
2*res2$value+2*length(res2$par)
```
```{r:Console}
1175.728
```

こちらも`glm()`関数と同じになっていることが確認できました．

## もう少し一般化する
上記のコードでは，説明変数の数が変わるたびに関数を定義しなければならないので
もう少し一般化したコードにします．
行列を使うとうまくいきます．

データの持ち方も少し変更します．

```{r}
Data_mat<-as.matrix(Data[,-1])
Data_mat<-cbind(1,Data_mat)
Data_mat
```

```{r:Console}
##           pclass_2nd pclass_3rd sex_male      age
##    [1,] 1          0          0        0 29.00000
##    [2,] 1          0          0        0  2.00000
##    [3,] 1          0          0        1 30.00000
##    [4,] 1          0          0        0 25.00000
##    [5,] 1          0          0        1  0.91670
##    [6,] 1          0          0        1 47.00000
##    [7,] 1          0          0        0 63.00000
##    [8,] 1          0          0        1 39.00000
##    [9,] 1          0          0        0 58.00000
##   [10,] 1          0          0        1 71.00000
##   [11,] 1          0          0        1 47.00000
##   [12,] 1          0          0        0 19.00000
##   [13,] 1          0          0        0 31.19418
##   [14,] 1          0          0        1 31.19418
##   [15,] 1          0          0        1 31.19418
##   [16,] 1          0          0        0 50.00000
##   [17,] 1          0          0        1 24.00000
##   [18,] 1          0          0        1 36.00000
##  [ reached getOption("max.print") -- 1295 行を無視しました ]
```

一列目は切片用にすべて1としています．
先ほどの関数を行列を使って定義しなおします．

```{r:Negative_Log_likelihood_function.R} 
logistic_reg<-function(x,beta){
   logistic(x=x%*%beta)
 }
  
#負の対数尤度関数
nll_logistic<-function(para,y,x){
   tmp<-log(ifelse(y==1,logistic_reg(x=x,beta=para),
                   1-logistic_reg(x=x,beta=para)))
   -1*sum(tmp)
 }
```
最尤法でパラメータを推定します．

```{r}
 inits<-res$coefficients
 res2<-optim(inits,nll_logistic,y=Data$survived,x=Data_mat,
             method = "BFGS",hessian = T) 

 res_mat<-data.matrix(matrix(NA,2,length(res2$par),
                    dimnames = list(c("estimates","std_error"),
                                    names(res2$par))))
 res_mat[1,]<-res2$par
 res_mat[2,]<-sqrt(diag(solve(res2$hessian)))
 res_mat
```
```{r:Console}
          (Intercept)  pclass2nd  pclass3rd   sexmale          age
estimates   3.5522611 -1.1707769 -2.4306725 -2.463377 -0.042235350
std_error   0.3421664  0.2115554  0.1951519  0.154585  0.007415443
```
推定値や標準誤差も一致していました．
最後にAICも確認しておきます．

```{r:AIC}
 2*res2$value+2*length(res2$par) 
```
```{r:Console}
[1] 1175.728
```
無事，同じ値になっていました．

# まとめ
ロジスティック回帰を実行できる関数を定義し，`glm()`関数と同じ答えになるか
確かめました．

自分の手でコードを書いてみると，意外と理解できていないところが
多いことに気づけました．

ご意見や誤りなどあればコメントをお願いします．

# 参考文献

[ロジスティック回帰を理解する(2)_パラメータ推定について](https://qiita.com/g-k/items/7f9eb73d088d35b7ab63)
[R Language タイタニックデータセットのロジスティック回帰](https://riptutorial.com/ja/r/example/9800/%E3%82%BF%E3%82%A4%E3%82%BF%E3%83%8B%E3%83%83%E3%82%AF%E3%83%87%E3%83%BC%E3%82%BF%E3%82%BB%E3%83%83%E3%83%88%E3%81%AE%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E5%9B%9E%E5%B8%B0)