ロジスティック回帰分析の基礎 - 統計学備忘録（R言語のメモ）

投稿日:2018.2.13 最終更新日:2020.07.26

説明変数がが単一かつ連続変数の場合

ロジット関数 $\ \$ $y=b0+b1*x = log\frac{P}{1-P} \ \ (0\verb|＜|P\verb|＜|1)\ \ \ →\ \ \ -∞\verb|＜| f(P)\verb|＜|∞$

(標準)ロジスティック関数 $\ \$ ロジット関数の逆関数= $P = \frac{exp(y)}{1+exp(y)}$

サンプルirisより

品種"virginica＝1"、"別の品種＝0"という２値を目的変数、”Sepal.Length(がく片の長さ)(cm)”を説明変数とします．ロジスティック回帰で求めることは、”Sepal.Length(がく片の長さ)(cm)”から予測される、品種が"virginica"になる確率です．

class(iris$Species)  #属性を確認
is <- as.numeric(iris$Species)   #  iris$Speciesが要因になっているので数字に変換（as.integerでも同じ）
y <- gsub(1, 0, is)  # is の1を0に変換
y <- gsub(2, 0, y)  # is の2を0に変換
y <- gsub(3, 1, y)  # is の3を1に変換
y <- as.numeric(y)#数字に変換：これを忘れると後でエラーが出ます
x <- iris$Sepal.Length
df <- data.frame(y, 1-y, x)
colnames(df)<- c("y", "1-y", "x")
edit(df)  #データフレームの表示
plot(df[,c("x","y")])  #散布図の表示

f:id:yoshida931:20200725222848p:plain:w250
f:id:yoshida931:20200726103742p:plain:w300
ちなみに、そのまま単回帰分析を行ってみます．

summary(res <- lm(formula=y~x, data=df) )
abline(res) 
#lines(range(df$x), res$coef[1] + res$coef[2]*range(df$x)) と同じ直線

f:id:yoshida931:20180208184829p:plain:w300
h 目的変数は0、1なのですが、上図ではかなりはみ出してしまいます．0より少なくても、1より大きくても、良い予測とは言えません．したがって、この回帰直線は適当なモデルとは言えません．そこで指数を用いたロジスティックモデルを使用します．

ロジスティック変換 logistic transformation

目的変数が2値変数（品種=virginica＝1、品種=別の品種＝0）の回帰分析をおこなうためには、2値変数に与えられる上限と下限を取り除けば解決します．ロジスティック変換は、0～1の区間制限を取り除くことで、統計的推測に柔軟性を持たせた方法といえます．

品種が"virginica（y=1）"になる確率を $P$ とした場合のロジスティック（ロジット）変換を考えてみます．

$オッズ=\frac{P(y =1 | xi)}{P(y = 0 | xi)}=\frac{p}{1-p}$

ロジット（logit）は、0から1の値をとる確率p に対し、そのオッズ $\frac{p}{1-p}$ の対数から計算される値．

Ｐのロジスティック変換 $\ \ logit(p)=log\frac{P}{1-P}$

$f(P)=log\frac{P}{1-P} \ \ (0\verb|＜|P\verb|＜|1)\ \ \ →\ \ \ -∞\verb|＜| f(P)\verb|＜|∞$

グラフで確認してみます

pl <- function(p){ return(log(p/(1-p))) }
plot(pl, xlab ="P", ylab = "f(P)")
abline(0,0)

ロジットf(p)がどのような値になっても、0＜ｐ＜ 1となっています．
$f(P) = b0+b1*x$
f:id:yoshida931:20180213133241p:plain:w300

ロジスティック回帰モデル

2値の変数(virginica＝1、別の品種＝0)が1になる確率のロジットを目的変数にしたものが、ロジスティック回帰分析になります。

ロジスティック回帰モデル
$f(P)=log\frac{P}{1-P}=b0 + b1 * x$

$P=\frac{exp(b0+b1 * x)}{1+exp(b0+b1 * x)}$

パラメータの推定（最尤推定法） $P(yi=1|b0,b1)$

尤度関数
$L(θ)= \prod_{i=0}^{150} (\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{yi})*(1-\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{1-yi}$

このままでは大変なので対数をとります
対数尤度関数 $L(θ)= \sum_{i=0}^{150} log(\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{yi}log(1-\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{1-yi}$ 　

yiは0か1です．例えば例題のy1と y150を見てみます．それぞれy1=1, y150=0となっています．
y1のときは
$f\theta (y_1)=log(1-\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{1-yi}$
y150のときは
$f\theta (y_{150})=log(\frac{exp(b0+b1x)}{1+exp(b0+b1x)})^{yi}$
となります．これを、1～150全て掛け合わせることになります．

ここからNewton法により最尤法を解いていくことになります．
qiita.com

Rを使ったロジスティック回帰分析

# glm関数で、応答変数の分布は二項分布なので引数をbinomialにすればOKのはず・・・
(ans <- glm(y~x,family=binomial))
#次も同じ式です
 (ans <- glm(cbind(y,1-y)~x,family=binomial))

Error in weights * y : non-numeric argument to binary operator
というエラーが出たので、例題で使用したｘｙの属性を調べます

typeof(x);typeof(y)
[1] "double"
[1] "character"

ｘは数で、ｙが文字になっていました．ということでｙを数に置きなおして再トライ

y <- as.numeric(y)
(ans <- glm(y~x,family=binomial))

Call:  glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
    -16.320        2.592  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3   AIC: 121.3

$b0=-16.32, b1=2.596$ という解になりました

#対数尤度を求めてみます
x1 <- x[101:150]
x2 <- x[1:100]
e1 <- exp(-(-16.3198+2.592*x1))
e2 <- exp(-(-16.3198+2.592*x2))
(ly <- sum(log(1/(1+e1)))+sum(log(1-1/(1+e2))))

[1] -58.67273

視覚的にイメージしておきます．

par(mfrow = c(1,2))
plot(x,y,yaxt="n",
     ylab = "目的変数(0,1)",
     xlab = "Sepal.Length(がく片の長さ)(cm)")     #まずは散布図
name<-c("0","1")
axis(side=2,at=c(0,1),labels=name)　
abline(lm(y~x))    #回帰直線
p <- function(x){ return(1/(1+exp(-(-16.320+2.592*x)))) } #ｐの関数を作ります
plot(p,-16.320+2.592*x,ylim=c(0,1),xlim = c(3.3,9.3),
     ylab = "virginicaである確率",
     xlab="b0+b1*x")

f:id:yoshida931:20180517150938p:plain:w700

Ｒを使用して最尤法で推定したパラメータの有意性の確認．

summary(ans)

Call:
glm(formula = y ~ x, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9870  -0.5520  -0.2614   0.5832   2.7001  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -16.3198     2.6581  -6.140 8.27e-10 ***
x             2.5921     0.4316   6.006 1.90e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 190.95  on 149  degrees of freedom
Residual deviance: 117.35  on 148  degrees of freedom
AIC: 121.35     #　-2*対数尤度+2*パラメータ数

Number of Fisher Scoring iterations: 5

標準誤差の算出については、またいつか・・・

オッズについて

ロジスティック変換より

がく片の長さxcmのオッズ＝ $\frac{P}{1-P} = exp(-16.3198+2.5921x)$

がく片の長さxcmから1cm増えた場合のオッズ＝ $\frac{P'}{1-P'} = exp(-16.3198+2.5921(x+1))$

オッズ比= $\frac{P(y=1 | x =xi +1)÷P(y=0 | x =xi +1)}{P(y=1 | x =xi )÷P(y=0 | x =xi)}$

$\frac{\frac{P'}{1-P'}}{ \frac{P}{1-P}}= exp(-16.3198+2.5921(x+1)-(-16.3198+2.5921x))=exp(2.5921)=13.35646$

ｘが1増えた場合の対数オッズ比は2.59、オッズ比は13.36になります．

 summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900

ｘには大きな変動がないので、ｘが0.5増えた場合のオッズ比を求めてみます

$\frac{\frac{P'}{1-P'}}{ \frac{P}{1-P}}=exp(2.5921×0.5)=3.654832$

つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます．オッズ比にはｘが含まれていないので、どの幅(xcm)のに対しても同じことが言えます．
さらに切片を解釈するためにｘから６を除した変数x6＝ｘ－６を追加してみます． f:id:yoshida931:20200726162414p:plain:w250

(ans <- glm(cbind(y,1-y)~x6,family=binomial,data=df))

Call:  glm(formula = cbind(y, 1 - y) ~ x6, family = binomial, data = df)

Coefficients:
(Intercept)           x6  
    -0.7675       2.5921  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3        AIC: 121.3

$P(y=1 | x=6) = exp(-0.7675+0) / (1 +exp(-0.7675+0) ) = 0.3170202$
がく片の長さが６cmの場合に、品種が"virginica"である確率は約30%と言える．