はじめに
ロジスティック回帰をツールを使ってやる機会があったので勉強用に
あえてRのglm()をつかわずにやってみようと思い記事にしました．
実行環境
今回の実行環境です．
R.Version()$version.string
[1] "R version 3.6.2 (2019-12-12)"
ライブラリは次のものを使いました．
library(tidyverse)
ロジスティック回帰とは
次のような直線があるとします．
y = ax+b
ここでデータ，$(x_i,y_i)$を示す点とこの直線の誤差 $e_i$ を次のように定義します.
e_i = y_i-(ax_i+b)
この誤差の二乗和を最小にする$a$,$b$を求めた時，直線 $y = ax+b$を，$(x_i,y_i)$ の $y$ の $x$ への回帰直線と呼びます．
しかし，応答変数が不良品率のような確率をとる時この直線のまま考えるのは適切でない場合があります．
（説明変数を大きく/小さくすると応答変数が0または1を超えてしまうため．）
そこで，次のようなロジットモデルが用いられます．
L(p) = ln\frac{p}{1-p}=ax+b
上の式を $p$ について解くと，
p = \frac{1}{1+exp(-ax-b)}
となり$x$の値に関わらず， $p$ が $[0,1]$ からはみ出すことがなくなります．
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)
ロジスティック回帰では，応答変数が $p$ という確率で1となり，確率 $1-p$ で0になる二項分布について考えます．
パラメータの推定には最尤法を用います．
尤度関数 $L(\theta)$ は次のようになります．
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$ も同じなので
対数をとった対数尤度関数を最大化します．
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を用います．
url <- "http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic.txt"
titanic <- read_csv(file = url)
変数として年齢を使おうと思ったのですが，欠損値が多いので
補います．ここでは単純に，利用できるデータの平均をとって代入しています．
titanic$age <- ifelse(is.na(titanic$age),mean(titanic$age,na.rm = T),titanic$age)
まずはRのglm()でロジスティック回帰を実行します．
説明変数にはクラス，性別，年齢を使ってみます．
res <- glm(survived ~ pclass + sex + age,
family = binomial, data = titanic)
summary(res)
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の数値に置き換えます．
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)
# 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
負の対数尤度関数を定義します．
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関数を使います．
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
(Intercept) pclass2nd pclass3rd sexmale age
3.55226110 -1.17077691 -2.43067246 -2.46337730 -0.04223535
同じ推定値になりました．
AICも確認します．
2*res2$value+2*length(res2$par)
1175.728
こちらもglm()関数と同じになっていることが確認できました．
もう少し一般化する
上記のコードでは，説明変数の数が変わるたびに関数を定義しなければならないので
もう少し一般化したコードにします．
行列を使うとうまくいきます．
データの持ち方も少し変更します．
Data_mat<-as.matrix(Data[,-1])
Data_mat<-cbind(1,Data_mat)
Data_mat
## 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としています．
先ほどの関数を行列を使って定義しなおします．
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)
}
最尤法でパラメータを推定します．
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
(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も確認しておきます．
2*res2$value+2*length(res2$par)
[1] 1175.728
無事，同じ値になっていました．
まとめ
ロジスティック回帰を実行できる関数を定義し，glm()関数と同じ答えになるか
確かめました．
自分の手でコードを書いてみると，意外と理解できていないところが
多いことに気づけました．
ご意見や誤りなどあればコメントをお願いします．
参考文献
ロジスティック回帰を理解する(2)_パラメータ推定について
R Language タイタニックデータセットのロジスティック回帰