Chapter 4 コホート研究
4.1 ロジスティック回帰
ロジスティック回帰(logistic regression) は、1948年にアメリカにおいてフラミンガム心臓研究（Framingham Heart Study）のために開発された統計手法です。 これは、冠状動脈性疾患に関する大規模なコホート研究でした。
一般的には、以下のような流れになります。
- 関連因子 \(x\) から、目的変数 \(y\) が発生するオッズ比を求める。
- 目的変数 \(y\) は、２値変数 (例えば「生」と「死」)。
- 関連因子 \(x\) は、同時に複数の分析ができる（\(x_1\)、\(x_2\) のように添え字をつけて区別する）。
- 関連因子 \(x\) は、連続変数、２値以上の変数をとることができる。
解析の結果は、オッズ比（またはオッズ比の対数）が得られます。では、オッズ比とは何でしょうか。
リスク比、オッズ比、ハザード比は、統計上の重要な指標です。
まず、理解しやすいリスク比から始めます。二つの例で計算します。ある疾患に対して、性別による差があるかどうかを調べてみたとします。
例1: 1100人のうち100人が死亡が死亡しました。生存している1000人のうち500人が男性でした。死亡した100人のうち50人が男性でした。
男性の女性に対するリスク比 = (50/100) / (500/1000) = 1.0
例2: 1100人のうち100人が死亡が死亡しました。生存している1000人のうち500人が男性でした。死亡した100人のうち90人が男性でした。
男性の女性に対するリスク比 = (590/100) / (500/1000) = 1.8
このように、「ある要因の影響がない場合、リスク比は 1.0」となります。また、「正の影響がある場合は 1 より大きくなり、負の影響がある場合は 0 に近づく」という性質があります。
オッズ比も同様に、「ある要因の影響がない場合、リスク比は 1.0」となります。
例1: 男性の女性に対するオッズ比 = (50/50) / (500/500) = 1.0
例2: 男性の女性に対するオッズ比 = (90/10) / (500/500) = 9.0
この例で見るように、影響がある場合、リスク比より１から大きく離れる傾向があります。
関連因子の同定は、統計学的には単変量解析で p < 0.10 またはp < 0.15 のものを選ぶことがあります。このほか、先行文献などから関連があるとされているものを見つけることもあります。
ln OR\(_y\) = \(\alpha\) + ln OR\(_1\) \(\times\) \(x_1\) + ln OR\(_2\) \(x_2\) + \(\cdots\)
対数は、高校数学IIで学びます。対数の底は自然対数 \(e\) です。
OR \(_i\)：
- \(x_i\) が２値変数の時(ex. 男=0 or 女=1)、\(y\) の男に対する女のオッズ比。 なお、３値以上の因数型 (A, B, C, …) の場合、(A, B), (A, C) … と２値変数に変換します。
- \(x_i\) が連続変数の時、\(y\) の \(x_i\) が1大きくなる時のオッズ比。
この式が論文中に出てくることはあまりありません。 通常は、表形式で以下のように示されているでしょう。
| OR | 95%CI | 型 | |
|---|---|---|---|
| 社会参加数 | 0.83 | (0.74-0.96) | 数値 |
| 性別 (女性) | 1.16 | (0.81-1.68) | 因子 |
| 年齢 ≧65 | 1.11 | (0.75-1.62) | 因子 |
| BMI | 1.11 | (0.75-1.62) | 数値 |
| 喫煙 | 1.18 | (0.56-2.50) | 因子 |
上述の表の場合、以下のような解釈になります。 社会参加数が１増えると、イベント発生のオッズ比が0.83となります。 これは、95％信頼区間 (CI) が 1 を含んでおらず、有意差であると言えます。 年齢が１歳上がるとイベント発生のオッズ比が 1.11 となります、性別（女性）は、イベント発生のオッズ比が 1.16 となりますが、95%CI が 1 を含んでいるので、有意差があるとは言えないと判断します。
4.2 論文 Abe T (2020) PLoS One
今回は、実際に論文に掲載されたデータを使って、ロジスティック回帰分析を行ってみましょう。
- Abe T, Okuyama K, Kamada M, Yano S, Toyama Y, Isomura M, … & Miyazaki R (2020). Social participation and physical prefrailty in older Japanese adults: The Shimane CoHRE study. PloS one, 15(12), e0243548. https://doi.org/10.1371/journal.pone.0243548
この論文は、島根県隠岐島において行われた横断研究です。高齢者は身体的・精神的に衰える「フレイル」という状態になりやすいですが、これを放置してしまいますと要介護状態になります。最近では、社会的参加が多い高齢者はフレイルになりにくいという研究もあります。とはいえ「社会参加」という言葉は漠然としていて人によって解釈が異なります。そこで、日本の農村部における社会参加をより詳しく調べてみました。
社会参加としては、ボランティア活動、運動、町会活動、宗教、地域の高齢者サロンに分けました。それぞれの参加頻度を聞き取りました。結果としては、社会参加が多いほどプレフレイルにはなりにくく、種類としては運動がプレフレイルになりにくく、町会活動はフレイルとプレフレイルになりにくいという結果となりました。
ここで用いている統計手法はロジスティック回帰です。社会参加については、OR = 0.83 (95% CI, 0.74–0.94)でした。 (OR = 0.47; 95% CI, 0.31–0.73)。 (OR = 0.57; 95% CI, 0.37–0.86)。
4.3 xlsx データの読み込み
library(readxl)
download.file("https://doi.org/10.1371/journal.pone.0243548.s003", "journal_pone_0243548.xlsx", mode="wb")
dfFrailty <- read_excel("journal_pone_0243548.xlsx")この論文は、サイト上でエクセルファイルを提供しています。 関数でダウンロードし、保存します。 なお、エクセルファイルへの URL は、論文サイト中にあります。
関数
download.file(): インターネット上からダウンロードし、二つ目の引数で指定した名前で保存する。read_excel():readxlライブラリの関数。.xslx ファイルを読み込む。
ライブラリ名が read_xl なのに、関数名が read_excel となっているのは、困ったものです。
データを確認すると、ところどころ NA (欠損値) があるものの、比較的正規化されたデータであることが分かります。
念のため、因子型データを因子として定義します。
dfFrailty$frailty <- as.factor(dfFrailty$frailty)
dfFrailty$volunteer <- as.factor(dfFrailty$volunteer)
dfFrailty$sports <- as.factor(dfFrailty$sports)
dfFrailty$neighborhood <- as.factor(dfFrailty$neighborhood)
dfFrailty$religious <- as.factor(dfFrailty$religious)
dfFrailty$salon <- as.factor(dfFrailty$salon)
dfFrailty$gender <- as.factor(dfFrailty$gender)
dfFrailty$smoking <- as.factor(dfFrailty$smoking)
dfFrailty$medication <- as.factor(dfFrailty$medication)
dfFrailty$age <- as.factor(dfFrailty$age)
dfFrailty$education <- as.factor(dfFrailty$education)
dfFrailty$working <- as.factor(dfFrailty$working)
dfFrailty$livingArrengement <- as.factor(dfFrailty$livingArrengement)同じことをなんども繰り返す時、通常はループ処理を行います。 本書では、極力プログラミング的なことを避けるため、同じ処理でもそのまま記述しています。
関数
as.factor(): 数値型や文字型を因子 (factor) 型にする。
4.3.1 外れ値 (outlier) と代入 (imputation)
論文によると、代入 (imputation) をしたようなので、Imputation_Detaset という列があります。 どうやら、Imputation_Detaset == 1 となっているデータは代入前で、Imputation_Detaset == 10 となっているのが解析対象データのようです。
ここでは、多重代入をやってみましょう。
R では、mice というパッケージを使用します。
関数
mice():miceライブラリの関数。代入に関するオブジェクトを作成する。complete():miceライブラリの関数。代入を行ってデータフレームを返す。
是非、上のコードを実行してみて、dfFrailtyImp0 と dfFrailtyImp0MICE を見比べてみてください。NA がなくなっていることが分かると思います。
代入前
代入後
4.4 表１の作成
library(tableone)
dfFrailty <- dfFrailtyImp10
CreateTableOne(data = dfFrailty,
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.30 (1.48) 1.98 (1.41) 1.79 (1.42) 0.012
## 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.011
## 1 104 (33.0) 76 (27.8) 4 (14.3)
## 2 118 (37.5) 98 (35.9) 7 (25.0)
## 3 93 (29.5) 99 (36.3) 17 (60.7)
## working = 2 (%) 87 (27.6) 74 (27.1) 4 (14.3) 0.308
## livingArrengement = 2 (%) 262 (83.2) 219 (80.2) 24 (85.7) 0.565
## 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 (%)4.5 ロジスティック回帰解析
4.5.1 データ整形
frail (frailty = 2) を削除し、robust (frailty = 0) と prefrail (frailty = 1) だけにします。
dfFrailty <- subset(dfFrailty, frailty != '2')4.5.2 ロジスティック回帰
library("epiDisplay", warn.conflicts=FALSE, quietly=TRUE)
GlmFrailty <- glm(frailty ~ TotalSP,
family=binomial(logit),
data=dfFrailty)
summary(GlmFrailty)##
## Call:
## glm(formula = frailty ~ TotalSP, family = binomial(logit), data = dfFrailty)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2563 -1.1260 -0.9427 1.2297 1.4319
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.18361 0.14802 1.240 0.21481
## TotalSP -0.15288 0.05766 -2.652 0.00801 **
## ---
## 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: 805.01 on 586 degrees of freedom
## AIC: 809.01
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmFrailty)##
## Logistic regression predicting frailty : 1 vs 0
##
## OR(95%CI) P(Wald's test) P(LR-test)
## TotalSP (cont. var.) 0.86 (0.77,0.96) 0.008 0.008
##
## Log-likelihood = -402.5047
## No. of observations = 588
## AIC value = 809.0094関数
glm(): stats ライブラリの関数。Generalized Linear Models (一般線形化モデル, GLM) を作成する。形式が独特で、“目的変数 ~ 説明変数1 + 説明変数2” という形式。多変量解析で用いる lm は二値変数を取らないのに対し、 glm は二値変数に対応している。summary(): base ライブラリの関数。モデルの統計値 (p 値や 95%信頼区間) などを返す。
GlmFrailtyM2 <- glm(frailty ~ TotalSP + gender + age + BMI + smoking + medication + education + working + livingArrengement,
family=binomial(logit),
data=dfFrailty)
summary(GlmFrailtyM2)##
## Call:
## glm(formula = frailty ~ TotalSP + gender + age + BMI + smoking +
## medication + education + working + livingArrengement, family = binomial(logit),
## data = dfFrailty)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4964 -1.0925 -0.8905 1.2177 1.5817
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.29116 0.82077 0.355 0.72278
## TotalSP -0.16273 0.05933 -2.743 0.00609 **
## gender2 -0.17045 0.18704 -0.911 0.36214
## age2 -0.10833 0.19657 -0.551 0.58157
## BMI 0.02519 0.02796 0.901 0.36766
## smoking2 -0.19633 0.38229 -0.514 0.60756
## medication1 -0.42510 0.23318 -1.823 0.06830 .
## medication2 -0.43096 0.23769 -1.813 0.06982 .
## education2 0.08037 0.20768 0.387 0.69877
## education3 0.27247 0.22089 1.234 0.21739
## working2 0.01984 0.20712 0.096 0.92367
## livingArrengement2 -0.17658 0.22227 -0.794 0.42694
## ---
## 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.00 on 576 degrees of freedom
## AIC: 818
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmFrailtyM2)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## TotalSP (cont. var.) 0.86 (0.77,0.96) 0.85 (0.76,0.95) 0.006
##
## gender: 2 vs 1 0.82 (0.58,1.14) 0.84 (0.58,1.22) 0.362
##
## age: 2 vs 1 0.86 (0.62,1.18) 0.9 (0.61,1.32) 0.582
##
## BMI (cont. var.) 1.03 (0.98,1.09) 1.03 (0.97,1.08) 0.368
##
## smoking: 2 vs 1 0.91 (0.46,1.81) 0.82 (0.39,1.74) 0.608
##
## medication: ref.=0
## 1 0.64 (0.41,1) 0.65 (0.41,1.03) 0.068
## 2 0.66 (0.42,1.02) 0.65 (0.41,1.04) 0.07
##
## education: ref.=1
## 2 1.14 (0.76,1.69) 1.08 (0.72,1.63) 0.699
## 3 1.46 (0.97,2.19) 1.31 (0.85,2.02) 0.217
##
## working: 2 vs 1 0.97 (0.68,1.4) 1.02 (0.68,1.53) 0.924
##
## livingArrengement: 2 vs 1 0.82 (0.54,1.25) 0.84 (0.54,1.3) 0.427
##
## P(LR-test)
## TotalSP (cont. var.) 0.006
##
## gender: 2 vs 1 0.362
##
## age: 2 vs 1 0.582
##
## BMI (cont. var.) 0.367
##
## smoking: 2 vs 1 0.608
##
## medication: ref.=0 0.133
## 1
## 2
##
## education: ref.=1 0.447
## 2
## 3
##
## working: 2 vs 1 0.924
##
## livingArrengement: 2 vs 1 0.427
##
## Log-likelihood = -397.0016
## No. of observations = 588
## AIC value = 818.0032参照（性別の場合、「男」か「女」）が違うのか、値が論文とは一致しません。
relevel() 関数を使って、性別の参照を 2 に変えてみます。
dfFrailty$gender <- relevel(dfFrailty$gender, ref=2)
dfFrailty$age <- relevel(dfFrailty$age, ref=2)
dfFrailty$smoking <- relevel(dfFrailty$smoking, ref=2)
dfFrailty$medication <- relevel(dfFrailty$medication, ref=2)
dfFrailty$education <- relevel(dfFrailty$education, ref=3)
dfFrailty$working <- relevel(dfFrailty$working, ref=2)
dfFrailty$livingArrengement <- relevel(dfFrailty$livingArrengement, ref=2)
GlmFrailtyM3 <- glm(frailty ~ TotalSP + gender + age + BMI + smoking + medication + education + working + livingArrengement,
family=binomial(logit),
data=dfFrailty)
summary(GlmFrailtyM3)##
## Call:
## glm(formula = frailty ~ TotalSP + gender + age + BMI + smoking +
## medication + education + working + livingArrengement, family = binomial(logit),
## data = dfFrailty)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4964 -1.0925 -0.8905 1.2177 1.5817
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.493316 0.751555 -0.656 0.51157
## TotalSP -0.162733 0.059325 -2.743 0.00609 **
## gender1 0.170448 0.187039 0.911 0.36214
## age1 0.108332 0.196574 0.551 0.58157
## BMI 0.025186 0.027958 0.901 0.36766
## smoking1 0.196332 0.382292 0.514 0.60756
## medication0 0.425100 0.233181 1.823 0.06830 .
## medication2 -0.005858 0.197126 -0.030 0.97629
## education1 -0.272474 0.220895 -1.234 0.21739
## education2 -0.192105 0.210573 -0.912 0.36161
## working1 -0.019844 0.207124 -0.096 0.92367
## livingArrengement1 0.176581 0.222272 0.794 0.42694
## ---
## 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.00 on 576 degrees of freedom
## AIC: 818
##
## Number of Fisher Scoring iterations: 4logistic.display(GlmFrailtyM3)##
## Logistic regression predicting frailty : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI)
## TotalSP (cont. var.) 0.86 (0.77,0.96) 0.85 (0.76,0.95)
##
## gender: 1 vs 2 1.23 (0.88,1.71) 1.19 (0.82,1.71)
##
## age: 1 vs 2 1.17 (0.84,1.62) 1.11 (0.76,1.64)
##
## BMI (cont. var.) 1.03 (0.98,1.09) 1.03 (0.97,1.08)
##
## smoking: 1 vs 2 1.1 (0.55,2.17) 1.22 (0.58,2.57)
##
## medication: ref.=1
## 0 1.56 (1,2.44) 1.53 (0.97,2.42)
## 2 1.0223 (0.7095,1.473) 0.9942 (0.6756,1.463)
##
## education: ref.=3
## 1 0.69 (0.46,1.03) 0.76 (0.49,1.17)
## 2 0.78 (0.53,1.15) 0.83 (0.55,1.25)
##
## working: 1 vs 2 1.03 (0.71,1.48) 0.98 (0.65,1.47)
##
## livingArrengement: 1 vs 2 1.22 (0.8,1.85) 1.19 (0.77,1.84)
##
## P(Wald's test) P(LR-test)
## TotalSP (cont. var.) 0.006 0.006
##
## gender: 1 vs 2 0.362 0.362
##
## age: 1 vs 2 0.582 0.582
##
## BMI (cont. var.) 0.368 0.367
##
## smoking: 1 vs 2 0.608 0.608
##
## medication: ref.=1 0.133
## 0 0.068
## 2 0.976
##
## education: ref.=3 0.447
## 1 0.217
## 2 0.362
##
## working: 1 vs 2 0.924 0.924
##
## livingArrengement: 1 vs 2 0.427 0.427
##
## Log-likelihood = -397.0016
## No. of observations = 588
## AIC value = 818.0032関数
relevel(): stats ライブラリの関数。因子 (factor) 型で、ref となる変数を設定する。See https://stat.ethz.ch/R-manual/R-devel/library/stats/html/relevel.html
4.6 モデルの評価
4.6.1 多重共線性 (multicollinearity)
Model 3 の多重共線性を調べる。
library(car, warn.conflicts=FALSE, quietly=TRUE)
vif(GlmFrailtyM3)## GVIF Df GVIF^(1/(2*Df))
## TotalSP 1.041654 1 1.020615
## gender 1.166296 1 1.079952
## age 1.368763 1 1.169941
## BMI 1.063765 1 1.031390
## smoking 1.178558 1 1.085614
## medication 1.150842 2 1.035748
## education 1.134117 2 1.031964
## working 1.209635 1 1.099834
## livingArrengement 1.049269 1 1.024338vif: 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 以下が適切とする論文もあります。（要確認）
4.7 サンプルサイズ
単変量のロジスティック回帰のサンプルサイズを計算するパッケージには以下のようなものがあります。
WebPowerpowerMediation
library(WebPower)
wp.logistic(p0 = 0.4, p1 = 0.5, alpha = 0.05, power = 0.80)library(powerMediation)
SSizeLogisticBin(p1 = 0.4, p2 = 0.5, B = 0.5, alpha = 0.05, power = 0.80)多変量解析では、最低限必要なイベント発生サンプル数 = 説明変数の数 × 10というサンプルサイズ計算式があります。
イベント発生率が0.5であれば、これに２をかける必要があります。
Peduzzi P, Concato J, Kemper E, Holford TR, Feinstein AR (1996) A simulation study of the number of events per variable in logistic regression analysis.: J Clin Epidemiol, 49(12):1373-9.