ロジステイック回帰
Yoshihiko Baba
2021/5/28
使用するライブラリ: readxl, tableone
今回は、実際に論文に掲載されたデータを使って、ロジスティック回帰分析を行ってみましょう。
はじめに
1948年にアメリカにおいてフラミンガム研究(Framingham study)のために開発された。 (冠状動脈性疾患に関する大規模なコホート研究)
目的変数 \(y\) は、２値変数 (ex. 生 or 死) のオッズ比
単変量解析 \(\rightarrow\) 関連因子の同定 \(\rightarrow\) 多変量解析
ln OR\(_y\) = \(\alpha\) + ln OR\(_1\) \(\times\) \(x_1\) + ln OR\(_2\) \(x_2\) + \(\cdots\)
対数の底は自然対数 \(e\)。
OR \(_i\)：
\(x_i\) が２値変数の時(ex. 男=0 or 女=1)、\(y\) の男に対する女のオッズ比
\(x_i\) が連続変数の時、\(y\) の \(x_i\) が1大きくなる時のオッズ比
この式が論文中に出てくることはほとんどない。 通常は、表形式で示されている。
| OR | 95%CI | |
|---|---|---|
| 性別 (女性) | 1.16 | 0.81-1.68 |
| 年齢 \(\geq\) | 1.11 | 0.75-1.62 |
| BMI | 1.11 | 0.75-1.62 |
解釈
式を見てわかる通り、一つのアウトカム (\(y\)) を、複数の要因の多変量解析としている。
まず、\(y\) は、アウトカムの発生オッズ比である。
\(x_i\) は、それぞれの説明変数である。 上述の表の場合、以下のような解釈になる。 年齢が１歳上がると、イベント発生のオッズ比が2.3となる。 性別（女性）であることは、男性であることよりもイベント発生のオッズ比が 1.7 となる。ただし、95%CIが 1 を含んでいるので、ここは有意差がない。
\(^{発展}\) 機械学習
機械学習でも２値変数の予測モデルがある。
ニューラルネットワーク
サポートベクトルマシン
決定木
ランダムフォレスト
ロジスティック回帰分析と機械学習の比較を行う研究が増えてきている。
- Maroco J, Silva D, Rodrigues A et al (2011) Data mining methods in the prediction of Dementia: A real-data comparison of the accuracy, sensitivity and specificity of linear discriminant analysis, logistic regression, neural networks, support vector machines, classification trees and random forests. BMC Research Notes, 4(1), 1-14.
データ
library(readxl)
strUrl <- "https://doi.org/10.1371/journal.pone.0243548.s003"
strDestfile <- "journal_pone_0243548.xlsx"
curl::curl_download(strUrl, strDestfile)
dfPone0243548 <- read_excel(strDestfile)
dfPone0243548 <- subset(dfPone0243548, Imputation_Detaset == 7)データを確認すると、ところどころ NA (欠損値) があるものの、比較的正規化されたデータであることが分かります。
論文によると、代入 (imputation) をしたようなので、Imputation_Detaset == 7 だけを抽出しました。
視覚的にも確認してみましょう。
library(GGally, warn.conflicts=FALSE, quietly=TRUE)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2dfPone0243548$frailty <- as.factor(dfPone0243548$frailty)
dfPone0243548$volunteer <- as.factor(dfPone0243548$volunteer)
dfPone0243548$sports <- as.factor(dfPone0243548$sports)
dfPone0243548$neighborhood <- as.factor(dfPone0243548$neighborhood)
dfPone0243548$religious <- as.factor(dfPone0243548$religious)
dfPone0243548$salon <- as.factor(dfPone0243548$salon)
dfPone0243548$gender <- as.factor(dfPone0243548$gender)
dfPone0243548$smoking <- as.factor(dfPone0243548$smoking)
dfPone0243548$medication <- as.factor(dfPone0243548$medication)
dfPone0243548$age <- as.factor(dfPone0243548$age)
dfPone0243548$education <- as.factor(dfPone0243548$education)
dfPone0243548$working <- as.factor(dfPone0243548$working)
dfPone0243548$livingArrengement <- as.factor(dfPone0243548$livingArrengement)
ggpairs(data = dfPone0243548, columns = 9:14, mapping = aes(color = frailty))
関数
- as.factor: base ライブラリの関数。数値型や文字型を因子 (factor) 型にする。
- ggpairs: GGally ライブラリの関数。R の基本 graphics パッケージの pairs を拡張している。
9:14 は、R 特有の表現で、c(9, 10, 11, 12, 13, 14) を意味します。
mapping = aes(color = frailty) の、frailty は dfPone0243548 の列名。frailty = 0, 1, 2 で、それぞれに色を割り当てる。
データの正規化
省略
外れ値 (outlier) と代入 (imputation)
省略
多重共線性と完全分離
省略
表１の作成
library(tableone)
CreateTableOne(data = dfPone0243548,
strata="frailty",
vars=c("TotalSP","volunteer","sports", "neighborhood","religious","salon","gender","age","BMI","smoking","medication","education","working","livingArrengement"),
factorVars=c("volunteer","sports", "neighborhood","religious","salon","gender","age","smoking","medication","education","working","livingArrengement"))## Stratified by frailty
## 0 1 2 p
## n 315 273 28
## TotalSP (mean (SD)) 2.25 (1.49) 1.92 (1.41) 1.68 (1.39) 0.006
## volunteer = 1 (%) 118 (40.3) 89 (36.9) 4 (18.2) 0.109
## sports = 1 (%) 102 (35.1) 50 (21.0) 3 (14.3) 0.001
## neighborhood = 1 (%) 205 (71.2) 139 (58.9) 6 (28.6) <0.001
## religious = 1 (%) 137 (47.2) 106 (45.7) 10 (45.5) 0.935
## salon = 1 (%) 91 (31.1) 64 (27.1) 13 (54.2) 0.022
## gender = 2 (%) 127 (40.3) 97 (35.5) 8 (28.6) 0.292
## age = 2 (%) 168 (53.3) 135 (49.5) 3 (10.7) <0.001
## BMI (mean (SD)) 22.83 (3.15) 23.14 (3.03) 22.50 (3.59) 0.350
## smoking = 2 (%) 297 (94.3) 256 (93.8) 27 (96.4) 0.841
## medication (%) 0.108
## 0 54 (17.1) 66 (24.2) 8 (28.6)
## 1 125 (39.7) 98 (35.9) 13 (46.4)
## 2 136 (43.2) 109 (39.9) 7 (25.0)
## education (%) 0.042
## 1 104 (33.0) 85 (31.1) 5 (17.9)
## 2 120 (38.1) 96 (35.2) 7 (25.0)
## 3 91 (28.9) 92 (33.7) 16 (57.1)
## working = 2 (%) 85 (27.0) 73 (26.7) 3 (10.7) 0.164
## livingArrengement = 2 (%) 260 (82.5) 216 (79.1) 23 (82.1) 0.567
## Stratified by frailty
## test
## n
## TotalSP (mean (SD))
## volunteer = 1 (%)
## sports = 1 (%)
## neighborhood = 1 (%)
## religious = 1 (%)
## salon = 1 (%)
## gender = 2 (%)
## age = 2 (%)
## BMI (mean (SD))
## smoking = 2 (%)
## medication (%)
## 0
## 1
## 2
## education (%)
## 1
## 2
## 3
## working = 2 (%)
## livingArrengement = 2 (%)手法
いくつかの方法があるが、もっとも代表的なのは以下の通り。
- ロジスティック回帰分析に先立ち、説明変数の選定を行うための「単変量解析」を行う。
- p＜0.15 の説明変数をロジスティック回帰分析に組み入れる。
例:
- Seiler S et al (2012) Driving Cessation and Dementia: Results of the Prospective Registry on Dementia in Austria (PRODEM), PLoS ONE, 7(12), e52710. [Open Access]
- Tanskanen M et al (2017) Population‐based analysis of pathological correlates of dementia in the oldest old. Annals of Clinical and Translational Neurology 4(.)3), 154-165.
- Mao HF et al (2016) Developing a referral protocol for community-based occupational therapy services in Taiwan: A logistic regression analysis. PloS one, 11(2), e0148414. [Open Access]
ただし、この論文はこのステップを経ていないようなので、飛ばします。
データ整形
frail (2) を削除し、robust (0) と prefrail (1) だけにする。
dfPone0243548 <- subset(dfPone0243548, frailty != '2')単変量解析
ロジスティック回帰
library("epiDisplay", warn.conflicts=FALSE, quietly=TRUE)
GlmPone0243548M1 <- glm(frailty ~ TotalSP,
family=binomial(logit),
data=dfPone0243548)
summary(GlmPone0243548M1)##
## Call:
## glm(formula = frailty ~ TotalSP, family = binomial(logit), data = dfPone0243548)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2598 -1.1228 -0.9309 1.2330 1.4458
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.19172 0.14534 1.319 0.18712
## TotalSP -0.16071 0.05759 -2.790 0.00526 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 812.14 on 587 degrees of freedom
## Residual deviance: 804.23 on 586 degrees of freedom
## AIC: 808.23
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmPone0243548M1)##
## Logistic regression predicting frailty : 1 vs 0
##
## OR(95%CI) P(Wald's test) P(LR-test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.005 0.005
##
## Log-likelihood = -402.1157
## No. of observations = 588
## AIC value = 808.2313関数
glm: stats ライブラリの関数。Generalized Linear Models (一般線形化モデル, GLM) を作成する。形式が独特で、“目的変数 ~ 説明変数1 + 説明変数2” という形式。多変量解析で用いる lm は二値変数を取らないのに対し、 glm は二値変数に対応している。
summary: base ライブラリの関数。モデルの統計値 (p 値や 95%信頼区間) などを返す。
GlmPone0243548M2 <- glm(frailty ~ TotalSP + gender + age + BMI + smoking + medication + education + working + livingArrengement,
family=binomial(logit),
data=dfPone0243548)
summary(GlmPone0243548M2)##
## Call:
## glm(formula = frailty ~ TotalSP + gender + age + BMI + smoking +
## medication + education + working + livingArrengement, family = binomial(logit),
## data = dfPone0243548)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5543 -1.1069 -0.8887 1.2080 1.6172
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.43495 0.82055 0.530 0.59606
## TotalSP -0.17304 0.05931 -2.918 0.00353 **
## gender2 -0.17534 0.18595 -0.943 0.34572
## age2 -0.14162 0.19190 -0.738 0.46053
## BMI 0.02702 0.02802 0.964 0.33498
## smoking2 -0.19764 0.38038 -0.520 0.60336
## medication1 -0.42838 0.23366 -1.833 0.06676 .
## medication2 -0.40753 0.23700 -1.720 0.08552 .
## education2 -0.07542 0.20539 -0.367 0.71346
## education3 0.05663 0.22021 0.257 0.79704
## working2 0.02042 0.20497 0.100 0.92066
## livingArrengement2 -0.21848 0.21927 -0.996 0.31905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 812.14 on 587 degrees of freedom
## Residual deviance: 794.29 on 576 degrees of freedom
## AIC: 818.29
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmPone0243548M2)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.84 (0.75,0.94) 0.004
##
## gender: 2 vs 1 0.82 (0.58,1.14) 0.84 (0.58,1.21) 0.346
##
## age: 2 vs 1 0.86 (0.62,1.18) 0.87 (0.6,1.26) 0.461
##
## BMI (cont. var.) 1.03 (0.98,1.09) 1.03 (0.97,1.09) 0.335
##
## smoking: 2 vs 1 0.91 (0.46,1.81) 0.82 (0.39,1.73) 0.603
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.65 (0.41,1.03) 0.067
## 2 0.66 (0.42,1.02) 0.67 (0.42,1.06) 0.086
##
## education: ref.=1
## 2 0.98 (0.66,1.45) 0.93 (0.62,1.39) 0.713
## 3 1.24 (0.82,1.86) 1.06 (0.69,1.63) 0.797
##
## working: 2 vs 1 0.99 (0.69,1.42) 1.02 (0.68,1.53) 0.921
##
## livingArrengement: 2 vs 1 0.8 (0.53,1.21) 0.8 (0.52,1.24) 0.319
##
## P(LR-test)
## TotalSP (cont. var.) 0.003
##
## gender: 2 vs 1 0.345
##
## age: 2 vs 1 0.46
##
## BMI (cont. var.) 0.335
##
## smoking: 2 vs 1 0.604
##
## medication: ref.=0 0.145
## 1
## 2
##
## education: ref.=1 0.823
## 2
## 3
##
## working: 2 vs 1 0.921
##
## livingArrengement: 2 vs 1 0.319
##
## Log-likelihood = -397.1467
## No. of observations = 588
## AIC value = 818.2934ref が違うのか、値が論文とは一致しません。 relevel 関数を使って、性別の ref を 2 に変えてみます。
dfPone0243548$gender <- relevel(dfPone0243548$gender, ref=2)
GlmPone0243548M2 <- glm(frailty ~ TotalSP + gender + age + BMI + smoking + medication + education + working + livingArrengement,
family=binomial(logit),
data=dfPone0243548)
summary(GlmPone0243548M2)##
## Call:
## glm(formula = frailty ~ TotalSP + gender + age + BMI + smoking +
## medication + education + working + livingArrengement, family = binomial(logit),
## data = dfPone0243548)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5543 -1.1069 -0.8887 1.2080 1.6172
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.25961 0.82296 0.315 0.75241
## TotalSP -0.17304 0.05931 -2.918 0.00353 **
## gender1 0.17534 0.18595 0.943 0.34572
## age2 -0.14162 0.19190 -0.738 0.46053
## BMI 0.02702 0.02802 0.964 0.33498
## smoking2 -0.19764 0.38038 -0.520 0.60336
## medication1 -0.42838 0.23366 -1.833 0.06676 .
## medication2 -0.40753 0.23700 -1.720 0.08552 .
## education2 -0.07542 0.20539 -0.367 0.71346
## education3 0.05663 0.22021 0.257 0.79704
## working2 0.02042 0.20497 0.100 0.92066
## livingArrengement2 -0.21848 0.21927 -0.996 0.31905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 812.14 on 587 degrees of freedom
## Residual deviance: 794.29 on 576 degrees of freedom
## AIC: 818.29
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmPone0243548M2)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.84 (0.75,0.94) 0.004
##
## gender: 1 vs 2 1.23 (0.88,1.71) 1.19 (0.83,1.72) 0.346
##
## age: 2 vs 1 0.86 (0.62,1.18) 0.87 (0.6,1.26) 0.461
##
## BMI (cont. var.) 1.03 (0.98,1.09) 1.03 (0.97,1.09) 0.335
##
## smoking: 2 vs 1 0.91 (0.46,1.81) 0.82 (0.39,1.73) 0.603
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.65 (0.41,1.03) 0.067
## 2 0.66 (0.42,1.02) 0.67 (0.42,1.06) 0.086
##
## education: ref.=1
## 2 0.98 (0.66,1.45) 0.93 (0.62,1.39) 0.713
## 3 1.24 (0.82,1.86) 1.06 (0.69,1.63) 0.797
##
## working: 2 vs 1 0.99 (0.69,1.42) 1.02 (0.68,1.53) 0.921
##
## livingArrengement: 2 vs 1 0.8 (0.53,1.21) 0.8 (0.52,1.24) 0.319
##
## P(LR-test)
## TotalSP (cont. var.) 0.003
##
## gender: 1 vs 2 0.345
##
## age: 2 vs 1 0.46
##
## BMI (cont. var.) 0.335
##
## smoking: 2 vs 1 0.604
##
## medication: ref.=0 0.145
## 1
## 2
##
## education: ref.=1 0.823
## 2
## 3
##
## working: 2 vs 1 0.921
##
## livingArrengement: 2 vs 1 0.319
##
## Log-likelihood = -397.1467
## No. of observations = 588
## AIC value = 818.2934関数
- relevel: stats ライブラリの関数。因子 (factor) 型で、ref となる変数を設定する。See https://stat.ethz.ch/R-manual/R-devel/library/stats/html/relevel.html
ステップワイズ
GlmPone0243548M2 <- glm(frailty ~ TotalSP + gender + age + BMI + smoking + medication + education + working + livingArrengement,
family=binomial(logit),
data=dfPone0243548)
GlmPone0243548M2Stepwise <- step(GlmPone0243548M2, direction = "both")logistic.display(GlmPone0243548M2Stepwise)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.84 (0.75,0.94) 0.003
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.63 (0.4,0.98) 0.041
## 2 0.66 (0.42,1.02) 0.61 (0.39,0.96) 0.031
##
## P(LR-test)
## TotalSP (cont. var.) 0.003
##
## medication: ref.=0 0.069
## 1
## 2
##
## Log-likelihood = -399.4363
## No. of observations = 588
## AIC value = 806.8725関数
- step: stats ライブラリの関数。Akaike’s Information Criteria (AIC) 方式によるステップワイズ。direction オプションで、“forward”, “backward”, “both” を選択できる。 See https://stat.ethz.ch/R-manual/R-devel/library/stats/html/step.html
社会活動数 (TotalSP) と 薬の種類数 (medication) だけになりました。
GlmPone0243548M2Stepwise <- step(GlmPone0243548M2, direction = "backward")logistic.display(GlmPone0243548M2Stepwise)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.84 (0.75,0.94) 0.003
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.63 (0.4,0.98) 0.041
## 2 0.66 (0.42,1.02) 0.61 (0.39,0.96) 0.031
##
## P(LR-test)
## TotalSP (cont. var.) 0.003
##
## medication: ref.=0 0.069
## 1
## 2
##
## Log-likelihood = -399.4363
## No. of observations = 588
## AIC value = 806.8725“backward” も、社会活動数 (TotalSP) と 薬の種類数 (medication) だけになりました。
GlmPone0243548M2Stepwise <- step(GlmPone0243548M2, direction = "forward")logistic.display(GlmPone0243548M2Stepwise)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.85 (0.76,0.95) 0.84 (0.75,0.94) 0.004
##
## gender: 1 vs 2 1.23 (0.88,1.71) 1.19 (0.83,1.72) 0.346
##
## age: 2 vs 1 0.86 (0.62,1.18) 0.87 (0.6,1.26) 0.461
##
## BMI (cont. var.) 1.03 (0.98,1.09) 1.03 (0.97,1.09) 0.335
##
## smoking: 2 vs 1 0.91 (0.46,1.81) 0.82 (0.39,1.73) 0.603
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.65 (0.41,1.03) 0.067
## 2 0.66 (0.42,1.02) 0.67 (0.42,1.06) 0.086
##
## education: ref.=1
## 2 0.98 (0.66,1.45) 0.93 (0.62,1.39) 0.713
## 3 1.24 (0.82,1.86) 1.06 (0.69,1.63) 0.797
##
## working: 2 vs 1 0.99 (0.69,1.42) 1.02 (0.68,1.53) 0.921
##
## livingArrengement: 2 vs 1 0.8 (0.53,1.21) 0.8 (0.52,1.24) 0.319
##
## P(LR-test)
## TotalSP (cont. var.) 0.003
##
## gender: 1 vs 2 0.345
##
## age: 2 vs 1 0.46
##
## BMI (cont. var.) 0.335
##
## smoking: 2 vs 1 0.604
##
## medication: ref.=0 0.145
## 1
## 2
##
## education: ref.=1 0.823
## 2
## 3
##
## working: 2 vs 1 0.921
##
## livingArrengement: 2 vs 1 0.319
##
## Log-likelihood = -397.1467
## No. of observations = 588
## AIC value = 818.2934“foward” は、前の二つとは大きく異なる結果となりました。
モデルの評価
多重共線性 (multicollinearity)
Model 2 (stepwise せず) の多重共線性を調べる。
library(car, warn.conflicts=FALSE, quietly=TRUE)
vif(GlmPone0243548M2)## GVIF Df GVIF^(1/(2*Df))
## TotalSP 1.046141 1 1.022810
## gender 1.153660 1 1.074086
## age 1.305097 1 1.142409
## BMI 1.070778 1 1.034784
## smoking 1.167293 1 1.080413
## medication 1.142855 2 1.033946
## education 1.143508 2 1.034093
## working 1.173773 1 1.083408
## livingArrengement 1.055319 1 1.027287vif: car パッケージの関数。
なお、VIF 値を求める関数は、多くのパッケージで提供されています。 具体的には、AED (corvif), car, rms, DAAG.
この表の一番右側を見ます。
- Fox J & Monette G (1992) Generalized Collinearity Diagnostics, JASA, 87, 178-183.
VIF 値の解釈は、多変量解析の場合、10以下で適切、10以上は多重共線性により不適切。ただし、ロジスティック回帰では、VIF の解釈はできないとする論文や、2.5 以下が適切とする論文もある。（要確認）
適合度 Hosmer-Lemeshow 検定
省略