医療データ奮闘記

公衆衛生大学院に入った内科系専門医が医師として培った現場感と大学院で培った統計の知識を交えながら、医療や疫学や統計に関する素朴な疑問や本音をつらつら書いています。

なんでもかんでも交絡因子を調整してはいけない

統計 R 後輩レクチャー

なんでもかんでもいれてはいけない

医学分野の論文で、「これが代替の説明変数かな、とか、とりあえず多めに調整しておけば良いか、と思って説明変数を多数入れる」事がよく行われている。1

collider bias （日本語では、合流点バイアス？）というが、あまり知られていない。

実際何も考えずに行うと、どうなるか実験する。2

DAG (Directed Acyclic Graph)

f:id:yiliaojingji:20190820112447p:plain
DAG_collider

想定しているモデルの因果ダイアグラムを描いておく

意図としては、交通事故の発生確率(連続変数)と年齢の関係を調べたいが、交絡因子として動脈硬化度(ASOや脳梗塞後などの上流因子として連続変数を想定)を考慮したい。

その際、腎機能に関してのデータがあったので、eGFRと動脈硬化度と年齢と腎機能低下に関するリスクは概ね解っていると仮定する。

腎機能と交通事故にはさすがに直接関係ないけれど、動脈硬化の度合いを他の説明変数で調整できないので、仕方なくeGFRを代替交絡因子として調整しようと考える。というようなイメージ

# サンプルサイズ
n.sample <- 5000
# それぞれの値を標準化したとして連続変数としてランダムに作成
age_v              <- rnorm(n.sample,mean=0,sd=1)
arteriosclerosis_v <- rnorm(n.sample,mean=0,sd=1)
totalrisk_v        <- rnorm(n.sample,mean=0,sd=1)
koutujiko_v <- 0.5*age_v + 1.5*arteriosclerosis_v + rnorm(n.sample,mean=0,sd=1)
egfr_v <- 1*age_v + 0.5*arteriosclerosis_v + 0.5*totalrisk_v + rnorm(n.sample,mean=0,sd=1)
data_v <- data.frame(AGE = age_v, ARTERIOSCLEROSIS = arteriosclerosis_v, TOTALRISK = totalrisk_v
                             ,KOTSUJIKO = koutujiko_v, eGFR = egfr_v)
head(data_v)

seed設定していないので各自確認。

res1 <- lm(KOTSUJIKO ~ AGE, data=data_v)
summary(res1)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE, data = data_v)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.496 -1.203 -0.001  1.219  7.214 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02460    0.02574  -0.955    0.339    
## AGE          0.51232    0.02560  20.010   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.82 on 4998 degrees of freedom
## Multiple R-squared:  0.07417,    Adjusted R-squared:  0.07399 
## F-statistic: 400.4 on 1 and 4998 DF,  p-value: < 2.2e-16

交通事故発生割合と年齢の間の関係は、初めに設定した真実(β =0.5)あたりの結果になる。

eGFRも関係あるかもしれへんやんと思ってしまうと、、

動脈硬化に関して調整しておいた方が良いかと思ったとする。まずは動脈硬化も入れて解析すると、

res2 <- lm(KOTSUJIKO ~ AGE + ARTERIOSCLEROSIS, data=data_v)
summary(res2)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + ARTERIOSCLEROSIS, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4910 -0.6877 -0.0116  0.6815  3.8000 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.007593   0.014210  -0.534    0.593    
## AGE               0.519588   0.014131  36.769   <2e-16 ***
## ARTERIOSCLEROSIS  1.513968   0.014173 106.819   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.004 on 4997 degrees of freedom
## Multiple R-squared:  0.718,  Adjusted R-squared:  0.7179 
## F-statistic:  6362 on 2 and 4997 DF,  p-value: < 2.2e-16

どちらも有意差はある。

もしeGFRを代理で用いると、

res3 <- lm(KOTSUJIKO ~ AGE + eGFR, data=data_v)
summary(res3)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + eGFR, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5586 -1.1474 -0.0204  1.1407  6.1140 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02056    0.02426  -0.847    0.397    
## AGE          0.01741    0.03117   0.558    0.577    
## eGFR         0.49286    0.01965  25.087   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.715 on 4997 degrees of freedom
## Multiple R-squared:  0.1777, Adjusted R-squared:  0.1774 
## F-statistic:   540 on 2 and 4997 DF,  p-value: < 2.2e-16

年齢の効果が完全に消えてしまう。

参考にeGFRのリスクで調整すると、

res4 <- lm(KOTSUJIKO ~ AGE + TOTALRISK, data=data_v)
summary(res4)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + TOTALRISK, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5377 -1.2012  0.0051  1.2126  7.2430 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02476    0.02574  -0.962    0.336    
## AGE          0.51238    0.02560  20.013   <2e-16 ***
## TOTALRISK    0.02504    0.02543   0.985    0.325    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.82 on 4997 degrees of freedom
## Multiple R-squared:  0.07435,    Adjusted R-squared:  0.07398 
## F-statistic: 200.7 on 2 and 4997 DF,  p-value: < 2.2e-16

こっちだとAGEに関しては正しい推定結果が出てくる。

まとめると

解析No 調整内容 交通事故に及ぼす年齢の影響
res1 単変量(年齢のみ) 有意差あり
res2 動脈硬化度で調整 有意差あり
res3 eGFRで調整 有意差なし
res4 eGFRのリスクで調整 有意差あり

個人的な意見としては、

DAG自体にも恣意性が入るので、

  • DAGは自分の頭の中で書いて、

  • この説明変数を除いた解析も入れておいたら文句を言われまい

という理論武装をしておく(modelを複数作る)のが良いと思う。

  1. 時としてreviewerに求められたりする↩

  2. http://takehiko-i-hayashi.hatenablog.com/entry/20130919/1379560929 を参考にさせて頂いた↩