理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

ロジスティック回帰分析の基礎

投稿日: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"という2値を目的変数、”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}の対数から計算される値.

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<p< 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
というエラーが出たので、例題で使用したxyの属性を調べます

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

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

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)))) } #pの関数を作ります
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

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

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

xが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 

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

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

つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます.オッズ比にはxが含まれていないので、どの幅(xcm)のに対しても同じことが言えます.
さらに切片を解釈するためにxから6を除した変数x6=x-6を追加してみます. 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
がく片の長さが6cmの場合に、品種が"virginica"である確率は約30%と言える.