これまでの記事を整理・更新したり、新たな記事を掲載します。 y2pt.com
ロジスティック回帰分析の基礎
これまでの記事を整理・更新したり、新たな記事を掲載します。 y2pt.com
説明変数がが単一かつ連続変数の場合
ロジット関数
(標準)ロジスティック関数ロジット関数の逆関数=
サンプル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")]) #散布図の表示
ちなみに、そのまま単回帰分析を行ってみます．
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 に対し、そのオッズの対数から計算される値．
Ｐのロジスティック変換
グラフで確認してみます
pl <- function(p){ return(log(p/(1-p))) } plot(pl, xlab ="P", ylab = "f(P)") abline(0,0)
ロジットf(p)がどのような値になっても、0＜ｐ＜ 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
というエラーが出たので、例題で使用したｘｙの属性を調べます
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
という解になりました
#対数尤度を求めてみます 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")
Ｒを使用して最尤法で推定したパラメータの有意性の確認．
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増えた場合のオッズ ＝
オッズ比=
ｘが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増えた場合のオッズ比を求めてみます
つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます．オッズ比にはｘが含まれていないので、どの幅(xcm)のに対しても同じことが言えます．
さらに切片を解釈するためにｘから６を除した変数x6＝ｘ－６を追加してみます．
(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
がく片の長さが６cmの場合に、品種が"virginica"である確率は約30%と言える．
回帰分析の結果の解釈
これまでの記事を整理・更新したり、新たな記事を掲載します。 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歳未満の人口が０％の場合には個人貯蓄合計が平均17.49660になる」と解釈できるのですが、15歳未満の人口が０％というのはあまり現実的な解釈ではありませんので、次のように書き換えます．
reg2 <- lm(sr ~ I(pop15-10), data=dat) summary(reg2)
結果の切片は、15歳未満の人口が1０％の場合の平均的な個人貯蓄合計を意味することになります。切片が変わるだけすので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歳未満の人口が１％増加するごとに個人貯蓄の合計が平均して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（残差標準誤差）…でも、残差の「標準偏差」を意味している
Multiple R-squared （決定係数）
R2=
決定係数は説明変数の数が増えるほど1に近づくので自由度で調整
Adjusted R-squared（自由度調整済み決定係数）
AR2=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.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
参考
次に最尤推定法による単回帰分析（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
t分布からの信頼区間
下記のサイトに移転しております