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

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

since2016 ときどきTEXのメモ

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

これまでの記事を整理・更新したり、新たな記事を掲載します。 y2pt.com

y2pt.com

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

ロジット関数 \ \  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")])  #散布図の表示



ちなみに、そのまま単回帰分析を行ってみます.

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)"になる確率を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

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

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") 

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を追加してみます.

(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%と言える.

回帰分析の結果の解釈

これまでの記事を整理・更新したり、新たな記事を掲載します。 y2pt.com

y2pt.com

単回帰分析

dat <- LifeCycleSavings      

# sr(個人貯蓄の合計), pop15(15歳未満の人口の数値%), pop75(75歳未満の人口の数値%), dpi(可処分所得、自由に使える収入), ddpi(dpiの成長率)

reg1 <- lm(sr ~ pop15, data=dat)    
summary(reg1) 

次のような結果が出力されます

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 17.49660    2.27972   7.675 6.85e-10 ***
pop15       -0.22302    0.06291  -3.545 0.000887 ***

結果の解釈について
Intercept=17.4966
「15歳未満の人口が0%の場合には個人貯蓄合計が平均17.49660になる」と解釈できるのですが、15歳未満の人口が0%というのはあまり現実的な解釈ではありませんので、次のように書き換えます.

reg2 <- lm(sr ~ I(pop15-10), data=dat)
summary(reg2)

結果の切片は、15歳未満の人口が10%の場合の平均的な個人貯蓄合計を意味することになります。切片が変わるだけすのでpop15の係数の推定値は変化していません。

              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   15.26642    1.67802   9.098 5.09e-12 ***
I(pop15 - 10) -0.22302    0.06291  -3.545 0.000887 ***

結果の解釈について
Intercept=15.26642
15歳未満の人口が10%の場合には個人貯蓄の合計が約15.3になり、さらに15歳未満の人口が1%増加するごとに個人貯蓄の合計が平均して0.22減少することが推定されます。
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

Residual standard error: 4.03 on 48 degrees of freedom
Multiple R-squared:  0.2075,    Adjusted R-squared:  0.191 
F-statistic: 12.57 on 1 and 48 DF,  p-value: 0.0008866

自由度 n=50, p=1 (単回帰なのでp=1となる)
Residual standard error(残差標準誤差)…でも、残差の「標準偏差」を意味している
  \sqrt{\frac{SSE}{自由度}}
Multiple R-squared (決定係数)
  R2= \frac{SSR}{SST}\

決定係数は説明変数の数が増えるほど1に近づくので自由度で調整
Adjusted R-squared(自由度調整済み決定係数)
  AR2=1- \frac{\frac{SSE}{n-p-1}}{\frac{SST}{n-1}}

F-statistic: 13.92 on 2 and 12 DF, p-value: 0.0007471
  MSR <- SSR/p
  MSE <- SSE/(n-p-1)
  F <- MSR/MSE    p-value=1-pf(F, p, n-p-1) ---> F分布から算出

分散分析表から求めることも可能

anova(reg2)

              Df Sum Sq Mean Sq F value    Pr(>F)    
I(pop15 - 10)  1 204.12  204.12  12.569 0.0008866 ***
Residuals     48 779.51   16.24          

SSR =  204.12     
MSR =  204.12 / 1 =  204.12   
SSE =  779.51
MSE =  779.51 / (50-1-1) = 16.24

回帰式から直接求める方法
注) reg1の係数を使用します

b0 <- reg1$coeff[1] #切片
b1 <- reg1$coeff[2] #傾き
yhat <- b0 + b1*dat$pop15 #推定値(srの推定値)
heikin <- mean(dat$sr) #=推定値の平均
SSE <- sum((dat$sr - yhat)^2)
SSR <- sum((yhat - heikin)^2)
SST <- sum((dat$sr - heikin)^2)

リニューアル準備中のお知らせ

2022-9-26更新

現在のアクセス件数は24万3千件25万3千件です。
必要としていただいている皆様のために、只今リニューアルの準備中です。
しばらく(たぶん半年くらい・・・)投稿や更新が滞りますので予めご了承ください。

これまでの記事を整理・更新したり、新たな記事を掲載します。 y2pt.com

y2pt.com

最小二乗法と最尤推定法

投稿日:2022.3.9

irisを使用した単回帰分析
まずは最小二乗法

dat = iris
fit1 <- lm(Sepal.Length ~ Sepal.Width, data = dat)
summary(fit1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.5262     0.4789   13.63   <2e-16 ***
Sepal.Width  -0.2234     0.1551   -1.44    0.152    
---
Residual standard error: 0.8251 on 148 degrees of freedom
Multiple R-squared:  0.01382,   Adjusted R-squared:  0.007159 
F-statistic: 2.074 on 1 and 148 DF,  p-value: 0.1519

残差の標準誤差や決定係数などが算出されます。
またF値やp値は次に示すような分散分析の結果になります。

 anova(fit1)

             Df  Sum Sq Mean Sq F value Pr(>F)
Sepal.Width   1   1.412 1.41224  2.0744 0.1519
Residuals   148 100.756 0.68078

参考

yoshida931.hatenablog.com

次に最尤推定法による単回帰分析(lmがglmに変わるだけ・・・です)。もちろん答えは同じです。

fit2 <- glm(Sepal.Length ~ Sepal.Width, data = dat)
summary(fit2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.5262     0.4789   13.63   <2e-16 ***
Sepal.Width  -0.2234     0.1551   -1.44    0.152    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.6807844)

    Null deviance: 102.17  on 149  degrees of freedom
Residual deviance: 100.76  on 148  degrees of freedom
AIC: 371.99

ただし対数尤度を利用した解析結果ですのでdeviance、AICが算出されます。
確率分布は正規分布で、リンク関数は恒等リンクがデフォルトなので、本当は以下のように書かれています(算出結果は全く同じ)。

fit3 <- glm(Sepal.Length ~ Sepal.Width, family = gaussian (link = identity), data = dat)
summary(fit3)

family = gaussianのリンク関数は、 identity, log, inverseが使用できるようです。
The gaussian family accepts the links (as names) identity, log and inverse
ちなみにリンク関数を変えたら、(当たり前ですが)異なる結果となります。

fit4 <- glm(Sepal.Length ~ Sepal.Width, family = gaussian(link=log), data = dat)
summary(fit4)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.87900    0.08200  22.914   <2e-16 ***
Sepal.Width -0.03723    0.02668  -1.396    0.165    

(Dispersion parameter for gaussian family taken to be 0.6810537)
Null deviance: 102.17  on 149  degrees of freedom
Residual deviance: 100.80  on 148  degrees of freedom
AIC: 372.05