投稿日:2018.2.13 最終更新日:2020.07.26
説明変数がが単一かつ連続変数の場合
ロジット関数
(標準)ロジスティック関数ロジット関数の逆関数=
サンプル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")]) #散布図の表示
ちなみに、そのまま単回帰分析を行ってみます.
summary(res <- lm(formula=y~x, data=df) ) abline(res) #lines(range(df$x), res$coef[1] + res$coef[2]*range(df$x)) と同じ直線
h
目的変数は0、1なのですが、上図ではかなりはみ出してしまいます.0より少なくても、1より大きくても、良い予測とは言えません.したがって、この回帰直線は適当なモデルとは言えません.そこで指数を用いたロジスティックモデルを使用します.
ロジスティック変換 logistic transformation
目的変数が2値変数(品種=virginica=1、品種=別の品種=0)の回帰分析をおこなうためには、2値変数に与えられる上限と下限を取り除けば解決します.ロジスティック変換は、0~1の区間制限を取り除くことで、統計的推測に柔軟性を持たせた方法といえます.
品種が"virginica(y=1)"になる確率をとした場合のロジスティック(ロジット)変換 を考えてみます.
ロジット(logit)は、0から1の値をとる確率p に対し、そのオッズの対数から計算される値.
Pのロジスティック変換
グラフで確認してみます
pl <- function(p){ return(log(p/(1-p))) } plot(pl, xlab ="P", ylab = "f(P)") abline(0,0)
ロジットf(p)がどのような値になっても、0<p< 1となっています.
ロジスティック回帰モデル
2値の変数(virginica=1、別の品種=0)が1になる確率のロジットを目的変数にしたものが、ロジスティック回帰分析になります。
ロジスティック回帰モデル
パラメータの推定(最尤推定法) 
尤度関数
このままでは大変なので対数をとります
対数尤度関数
yiは0か1です.例えば例題のy1と y150を見てみます.それぞれy1=1, y150=0となっています.
y1のときは
y150のときは
となります.これを、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
という解になりました
#対数尤度を求めてみます 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")

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