線形回帰分析〜その6：モデルの作り方

前回の記事で、regression modelに入れるべき変数の選び方の基本が理解できたと思います。今回は、因果推論におけるモデル（「予測モデル」を作りたい訳ではないので注意してください）の作り方を解説します。

因果推論におけるモデルの作り方、変数の選び方は、幾つか存在します。

DAGを用いる方法

疫学的に最も好まれる方法であり、DAGに基づいてモデルに入れる因子を決定します。そして、そのDAGを描くためには、expert knowledgeが必要になります。

例えば、吸入麻酔が患者の術後せん妄に与える影響について調べたいとします。データ上は、暴露因子（exposure）=1が吸入麻酔薬となり、暴露因子（exposure）=0がプロポフォールなどの静脈麻酔薬、といった割付になります。DAGに基づいて交絡因子（confounder）を考える際、その麻酔方法の選択に影響を与える因子を考えなければなりません。麻酔科医は、患者の年齢、性別、PONVの既往、乗り物酔いの既往、smoking、喘息、COPD、心機能など、様々な因子を元に吸入麻酔薬にするか静脈麻酔薬にするか選択しています。すなわち、これらの因子がconfounderの可能性があり、モデルに入れるべき因子ということになります。DAGにおいてどのような矢印を描くかの判断は、疫学者や統計学者だけではできません。臨床の場で働く医師の知識（expert knowledge）が必要なんですね。

また、過去の研究も参考にする必要があります。過去の研究がせん妄に与える因子を明らかにしている場合は、そのような矢印を含んだDAGを描く必要があります。

そして、DAGによって決定されたconfoundersをモデルに入れ、その他にbiasを引き起こすようなpathwayがDAG上に無いことが確認できれば、モデル作成終了となります。

メリット

メリットとしては、multiple testingを減らせることです。DAGを用いてモデルに入れるべき変数を決めることにより、検定回数を減らし、multiple testing problemを可能な限り減らすことができます。

デメリット

Expert knowledgeにしても、過去の研究を参考にするにしても、最終的に変数を選ぶのは研究を行なっている人間であり、その意味では非常に主観的（subjective）です。特に、統計学者やデータサイエンティストといった、データや数字を相手に仕事をしている人間からすると、この方法には非常に抵抗があるようです。

全ての変数をmodelに入れる方法

自分がデータとして持っている全ての変数をモデルに入れてしまうやり方です。

メリット

誰が行なっても同じ結果になるため、客観的（subjective）な方法です。

デメリット

一方で、前述のように、regressionで解析して得られたp-valueの数だけ検定を行なっていることになるため、multiple testing problemは回避できていません。そして、もし変数の数が多ければ、βのstandard errorが大きくなり、不安定なモデル（overfitting）となってしまいます。

また、前回の記事で説明したように、collinear covariateをモデルに入れてしまうと検出力が低下してしまいます。持っている変数を全て入れてしまうと、その変数の中にcollinear covariateが入ってしまう可能性もあります。

単変量解析でスクリーニング

それぞれの変数とY (outcome)の関係性を解析することで、モデルに入れる変数を決める方法です。例えば、変数とYとの関連を示したp-valueが0.05（or 0.10 or 0.20）を下回った場合のみ、その変数をモデルに入れます。

メリット＆デメリット

Confounderはoutcomeとも関係があるはずなので、単変量解析でp<0.05となるような変数はconfounderである可能性が高いということになります。そのため、この方法ではできるだけ多くのconfounderをモデルに入れることが可能になります。一方で、厳密にはconfounderの定義にp-valueは関係ないため（βの変化、でしたよね？）、confounderの選択を間違う可能性があります。

Collinear covariateはある他の因子（ex. Significant predictor）を通してoutcomeと関係のある因子です。そのため、単変量解析ではp-valueが小さくなる可能性が高く、多くのcollinear covariateを見つけ出すことが可能になります。一方で、collinearityをもつ相方の変数が最終モデルにないと、collinearityそのものが失われてしまうため、本当にモデルから外すべき因子だったのかどうか不明となってしまいます。

Significant predictorはoutcomeと関係性があるはずなので、単変量解析でもp-valueは小さくなることが予想されます。しかし、significant predictorの中には他の変数と合わせて初めてsignificantとなるものがあるため、単変量解析で全てのsignificant predictorを見つけ出すことはできません。

デメリット

単変量解析でスクリーニングを行った全ての検定と、最終モデルに入れた変数の数だけp-valueが発生することになり、multiple testing problemが問題となります。

Automated Selection Algorithm

殆どの統計ソフトでは、自動的に必要な変数を見つけ出しモデルを作成するアルゴリズムを使用することができます。多くはp-valueやAICなどを元に変数を選択していきます。

メリット＆デメリット

客観的な方法である一方で、変数選びに施行された全ての検定数だけp-valueが生まれるため、multiple testing problemが問題となります。また、βの変化を見ていないので全てのconfounderを選べているか怪しいですし、collinear variableが多い場合には、誤って大事なsignificant predictorを除いてしまうなどということも起こり得ます。

例：Backward selection based on p-value

ここでは、p-valueの大きい変数から除いていくbackward eliminationを紹介します。

mod_full<-lm(bwt~low+age+lwt+race+smoke+ptl+ht+ui+ftv,data=df_bw)
mod_null<-lm(bwt~smoke,data=df_bw)
drop1(mod_full,scope=~low+age+lwt+race+ptl+ht+ui+ftv,test="F")

Outputは以下のようになります。

## Single term deletions
## 
## Model:
## bwt ~ low + age + lwt + race + smoke + ptl + ht + ui + ftv
##        Df Sum of Sq      RSS    AIC  F value    Pr(>F)    
## <none>              32721641 2313.7                       
## low     1  38151468 70873109 2457.8 200.5417 < 2.2e-16 ***
## age     1    296441 33018082 2313.4   1.5582 0.2136214    
## lwt     1    317044 33038685 2313.5   1.6665 0.1984566    
## race    2   1541268 34262909 2318.4   4.0508 0.0190956 *  
## ptl     3    498939 33220580 2310.5   0.8742 0.4556695    
## ht      1    353016 33074657 2313.7   1.8556 0.1749129    
## ui      1   2496430 35218072 2325.6  13.1224 0.0003838 ***
## ftv     5    148902 32870544 2304.5   0.1565 0.9778107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-valueの最も大きな変数を除きます。今回は、ftvが最もp-valueが大きいため、ftvを除きましょう。

# Remove the variable with the largest p-value ("ftv")
drop1(update(mod_full,~.-ftv),test="F")

すると、以下のようなoutputが得られます。

## Single term deletions
## 
## Model:
## bwt ~ low + age + lwt + race + smoke + ptl + ht + ui
##        Df Sum of Sq      RSS    AIC  F value    Pr(>F)    
## <none>              32870544 2304.5                       
## low     1  39567230 72437774 2451.9 213.0601 < 2.2e-16 ***
## age     1    313731 33184275 2304.3   1.6894  0.195374    
## lwt     1    336670 33207214 2304.5   1.8129  0.179884    
## race    2   1679583 34550126 2310.0   4.5221  0.012151 *  
## smoke   1   1083145 33953688 2308.7   5.8325  0.016751 *  
## ptl     3    560214 33430757 2301.7   1.0055  0.391684    
## ht      1    350600 33221143 2304.5   1.8879  0.171177    
## ui      1   2624926 35495470 2317.1  14.1346  0.000231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

次は、この中で最もp-valueの大きい変数を除きます。これを繰り返し、全て変数のp-valueが0.05未満となるまで続ける方法です。

例：Forward selection based on AIC

変数のないモデルから次々と加えていくforward selectionという方法もあります。AICに基づいて選択する場合、

# Using package
step(mod_null, scope=list(lower=~smoke,upper=mod_full),direction = "forward", k = 2)

とするだけで簡単にモデルが作成できます。

ハイブリッド法

上記の方法を組み合わせて、モデルを作成することができます。

例えば、DAGを用いて最低限必要な変数を選び、その後に統計ソフトに自動で追加の変数を選んでもらい、最後にβの変化が大きい変数をconfounderとして追加で救い上げる、といった方法です。

まとめ

因果推論におけるモデルの作成方法は、一つだけではありません。そして、どの方法を選ぶかは、因果推論における基本的な概念の理解が必要になります。そして、実はどの方法も一長一短で、最終的に選ぶのはauthorです。ですが、少なくともなぜその方法を選んだのか説明できるようにはなりたいですね。

 

⇒医療従事者に必要な統計学と疫学（目次）へ

References

• John Orav. BST 213: Applied Regression for Clinical Research. Harvard T.H. Chan School of Public Health