前回のIPWの記事の続きとして今回はSASでのIPW法での解析例について紹介をします。前回の記事はこちら。
個人的な意見ですが、IPW推定はプロシジャを利用してしまえば比較的簡単に実行が可能です。むしろそれ以上に時間をかけて考慮すべきなのは、変数間の関係、つまり調整すべき変数（交絡因子）が一体何であるかということかと思います。また、もちろんそれ以外にも様々な仮定が成り立っているかや、他のバイアスの存在、バイアスの方向性に関する議論などもあり、単純に解析を回せばいいわけではないなと特に最近痛感しています。
CAUSALTRTプロシジャ
CAUSALTRTプロシジャでは、次の3つのいずれかを実行可能です。
- 回帰モデル
- IPW
- AIPW
基本的な構文はこんな感じ。
PROC CAUSALTRT <options>;
BOOTSTRAP < options > ;
BY variables;
CLASS variables <(options)> …<variable<(options)>> </ global-options>;
FREQ variable;
MODEL outcome <(variable-options)><= <effects>> </ model-options>;
OUTPUT <OUT=SAS-data-set> <keyword=name …keyword=name>;
PSMODEL treatment <(variable-options)><= effects </ psmodel-options>>;
この赤く染めたところが回帰モデルに関する部分、青く染めた部分が傾向スコアに関する部分です。それぞれの手法の必要最低限の構文を紹介すると次のような感じになります。なお今回はL1, L2, L3がすべての交絡因子であると仮定します。
- 回帰モデル
modelステートメントで交絡因子を指定
psmodelステートメントで交絡因子を指定
- AIPW
model, psmodelステートメントで交絡因子を指定
オプションに関して
IPW推定でよく用いるオプションのみここでは紹介します。さらに詳細に知りたい場合には、SASのHelpより適宜確認していただければいいかなと思います。
- proc causaltrtステートメント
このステートメントではデータセットの指定（data=データセット名）や推定方法に関するオプションの設定などを行います。代表的なものには次のようなものがあります。- ATT
推定対象をATEではなく、ATT（介入群における平均因果効果）に変更。この際にはmethodはIPWRかREGADJに変更する必要あり。 - COVDIFFPS
傾向スコアモデルに組み込んだ共変量の重み付け、重みなしの標準化された平均差、分散比の計算と表示 - method=推定方法
メジャーなものには以下のようなものがあります。- AIPW
- IPW
- IPWR
- IPWS
- ATT
- psmodelステートメント
psmodelステートメントでは、傾向スコアのモデル式に関するオプションを指定します。基本的な書き方は、
psmodel 目的変数（介入を示す変数）＝説明変数（モデルに組み込む変数）
となります。なおこの傾向スコアの算出に当たってはCAUSALTRTプロシジャではロジスティック回帰を用います。
SAS公式のcausaltrtプロシジャに関してはこちらからどうぞ
分析例紹介
使用データ説明
解析で使用するデータセットは、Hernan＆Robinsによって2018年に行われたNHANES I Epidemiologic Follow-Up Studyのサブセットデータになります。
このデータセットに関するコードはSAS Help Center内で公開がされているので、必要な方は適宜そちらを参照ください。
この研究における目的は禁煙が体重に与える因果効果を測定することです。ベースライン時と10年後に健康診断が行われ、下記のような医学的、行動に関する項目について調査が行われています。
- Activity
日常生活における活動量（0, 1, 2）
- Age
1971年時点での年齢
- BaseWeight
1971年時点での年齢（kg）
- Change
ベースラインと10年後の体重変化（kg）
- Education
教育レベル（0, 1, 2, 3, 4, 5）
- Exercise
定期的に行う運動量（0, 1, 2）
- PerDay
1971年に1日に吸ったたばこの本数
- Quit
ベースラインと10年後の調査期間中の喫煙の有無（喫煙した場合に1、それ以外は0）
- Race
人種（白人は0 それ以外は1）
- Sex
性別（男性は0、女性は1）
- Weight
1981年時点の体重
- YearsSmoke
個人の喫煙歴
下表の通り、今回の解析でも用いるのは簡略化のために全1746オブザベーションのうち、欠損値を含まない1566オブザベーションです。この内訳としては、介入群のサンプルサイズが403、対照群のサイズが1163となっています。
解析内容
傾向スコアの算出に当たって、共変量として設定するのは以下のものです。
- Sex
- Age
- Education
- Exercise
- Activity
- YearsSmoke
- PerDay
なお曝露の有無を示す変数はQuitです。
上記の設定で傾向スコアの算出を行なった場合の両群の傾向スコアの分布はこんな感じです。それほど極端な値を含んでいるわけではないので、今回は特にサンプルを除かず、このまま解析に進んでいきます。
次に確認するのは、傾向スコアで条件づけた場合の共変量の分布です。
下図は傾向スコアで重み付けをする前の解析集団と、重み付けをした疑似集団（pseudo population）における共変量（Age, YearSmoke）のカーネル密度です。多少の誤差はあるにしても、重み付けをしたほうが共変量の分布が似たものになっていることがわかるかと思います。なおこのplotの作成は、plots=(PSDist pscovden(effects(Age YearsSmoke)))の部分で行っていますが、ここで指定できる共変量は連続値に限られます。
それ以外の共変量の分布も確認するために以下の標準化された差、分散比に関する表を出力し、その内容を確認します。標準化された差の値が0.1未満（分散比であれば1に近づく）となれば、一般的には両群でバランスが取れていると判断していいわけですが、重み付けした後はどの変数も十分バランスされていることがわかる表からみて取れるかと思います。この表は、proc causaltrtステートメントでcovdiffpsと書くことで計算、表示することができます。
上記の結果から、今回解析にあたって想定した傾向スコアモデルは妥当なものであるとし、IPW法を用いて平均因果効果（ATE）を推定しに行きます。
推定結果は下表の通り。
曝露群における推定値が4.9850、非曝露群における推定値が1.7954であることから、ATEは3.1896と推定がされました。
つまり、元のセッティングを振り返ってみると、喫煙による10年間の体重増加はおよそ3.19kgであると推測されるわけです。
傾向スコアマッチングとの比較
次に以前の記事で紹介した傾向スコアマッチングとの比較を行います。使用データは全く同じで、想定している共変量も同一です。記事についてはこちらから。
なお傾向スコアマッチングでは推定対象がATTとなっているので、傾向スコアを用いた逆確率重み付けに際しても、推定対象をATTとします。この場合proc causaltrtステートメントでは、ATT、method=IPWRと指定する必要があります。というわけでそれぞれの方法による出力結果を見てみます。傾向スコア、共変量調整前のATTの算出コードについては上記の過去記事に記載があるのでそちらを参照ください。
- 調整前
- 傾向スコアマッチング
- IPWR法
曝露群における平均因果効果（ATT）の値は、共変量で調整前は2.5406kg、傾向スコアマッチングでは3.2792kg、IPWR法では3.2756kgと推定がされています。傾向スコアで調整を行った後ろ二つの推定結果は非常に似た値（約3.28kg）となっており、未調整の結果（約2.54kg）と比較してその値が大きく増加していることがわかるかと思います。
コード
データセットの作成に関しては、使用データの項目をみてください。
proc causaltrt data=smokingweight covdiffps method=IPWS;
class Sex Race Education Exercise Activity Quit /desc;
psmodel Quit = Sex Age Education Exercise Activity YearsSmoke PerDay
/plots=(PSDist pscovden(effects(Age YearsSmoke)) );
model Change;
run;
/*IPWR法によるATTの推定*/
proc causaltrt data=smokingweight covdiffps method=IPWR ATT;
class Sex Race Education Exercise Activity Quit /desc;
psmodel Quit = Sex Age Education Exercise Activity YearsSmoke PerDay
/plots=(PSDist pscovden(effects(Age YearsSmoke)) );
model Change;
run;
参考
- SASによる因果推論 CAUSALTRTプロシジャの紹介
- SAS Help Center The CAUSALTRT Procedure