今回は傾向スコアマッチングについてまとめます。傾向スコアは観察研究やランダム化が行われない実験研究で大きなメリットを持つ方法です。マッチングにあたってはSAS14.2より追加されたPSMATCHプロシジャを使用してやっていきます。
傾向スコア自体の概要に関しては以前の記事を参考にどうぞ。
傾向スコアマッチングの概要
傾向スコア（PS; Propensity Score）とは、共変量で条件付けた際の個人が治療を受ける確率です。このスコアを用いたマッチングの利点としては、傾向スコアの算出に用いた観察された共変量Ｘについて両群でのバランスが期待される点です。
傾向スコアマッチングの基本的な流れは、
- 傾向スコアの算出
- ペアの作成
- 効果の推定
といった感じです。傾向スコアの算出方法としては最もメジャーなのはロジスティック回帰ですが、それ以外にも様々な統計、機械学習の手法が提案されています。ただあくまでもNo unmeasured confounders（未観測交絡がない）の仮定の下での議論になるので、その点には注意が必要です。マッチング後の両群での共変量の分布をみるなどといったチェックも必要になります。
また、そもそもマッチングという方法に対しての批判としてあげられるようにマッチングした集団とはいったい何なのかといったことや、結果の解釈には十分な考察が必要です。
DAGsの紹介
傾向スコアのモデルには交絡因子をすべて含むことが、因果効果の推定には必要な家庭となります。しかし実際には選択結果が妥当なものであるかの判断が難しい場合があります。その共変量の選択時に役立つ方法の一つがDAGsです。
DAGs自体については上の記事を参考にしていただければいいかなと思います。
簡単に言うと、DAGsを使いつつ、視覚的に変数同士の関係性を探っていくことでより妥当な共変量の選択につなげていこうということです。個人的な話ですが、自分の周りで医学研究をなさっている方も議論上ではしばしば使っていらっしゃいます。
また交絡因子ではないですが、因果効果の推定のばらつきを低減させるという目的でアウトカムを予測する因子を傾向スコアのモデルに組み込むことも有用です。そういった際にもDAGsは有用かなと思います。
マッチングに関して
マッチングと一言にいっても様々な手法、構成比がありますがここでは簡単に代表的なもののみ紹介します。
構成比
マッチングさせる介入（曝露）群と対照群の比の取り方の代表的なもの
- One-to-one matching
- One-to-Many matching
- Full matching
最もメジャーな構成比は1：1です。これは他のマペアの構成方法と比較して例数が少なくなりやすいというデメリットはありますが、マッチ後のそれぞれの群のサイズが等しくなります。そのためマッチングを利用しての推定における検出力の低下が最も軽微です。
マッチング方法
マッチングを行う際に使用するアルゴリズムとして主なものには、
- 貪欲マッチング（Greedy matching）
- 最適マッチング（Optimal matching）
があります。これら自体もかなり内容が重いので詳細については別記事等で参照していただけると大変助かります。
貪欲マッチングは、ある任意の個体に対して傾向スコアが最も近い非暴露群の個体を逐次に探索する方法ですが、これはさらに次のような方法があります。
- 最近傍法（Nearest neighbor）
介入（曝露）群から任意に選択された1個人に対して、対照群から最も近い傾向スコアの値を持つ被験者をペアとする
- キャリパーの設定
キャリパーについてはそれ単独ではなく、併用した形で用います。最近傍法で問題となるのはいくら傾向スコアの値がサンプルの中で最も近くても、値が離れていもいいのかということです。キャリパーはそのマッチングさせる際の傾向スコアの距離、許容範囲を指す言葉で、解析においては設定することが一般的です。おおよそ傾向スコアの推定値の対数をとった標準偏差の0.2、0.25倍が推奨されています。ちなみにPSMATCHでのデフォルト設定は0.25です。
PSMATCHプロシジャ
まずは今回の解析に当たって使用するSASのPSMATCHプロシジャの基本的なステートメントを見ていきます。
それぞれのオプションについて代表的なものをいくつか
class
分類を行う変数の指定
psmodel
傾向スコアを算出するロジスティック回帰モデルの設定を行うオプションです。
PSMODEL （曝露を示す二値変数）<(Treated='曝露を意味する値')> = （共変量）
match
- CALIPER=
キャリパーの大きさ設定
- DISTANCE=
使用する距離の種類を設定
- METHOD=
マッチング方法の設定
- WEIGHT=
Weightの設定
asses
- ALLCOV
PSMODELステートメントで指定した共変量の分布を評価（二値変数or連続変数） - LPS
傾向スコアのロジット変換の差を評価 - PS
傾向スコアを評価 - VAR =（共変量）
指定した共変量の分布を評価（二値変数or連続変数）PSMODELステートメントで指定されていない変数も可
これらのオプションの指定の後に/（スラッシュ）の後にいくつかのオプションをさらに指定できます。例えば、
- plot=
バランスの評価のプロットを作成
- STDBINVAR =
Standardized Mean Differencesにおける二値変数を標準化の有無
- STDDEV =
Standardized Mean Differencesで利用する標準偏差の種類を指定
- VARINFO
治療群と対照群の変数情報の表示
output
アウトプットステートメント、いつも通りです。
より詳細なオプションについて知りたい場合はこちらからどうぞ。
分析例紹介（SmokingWeight）
PSMATCHプロシジャを使ったATTの推定の傾向スコアマッチングの例として紹介するのは、この動画（Tools for Assessing a Propensity Score Model in SAS/STAT）で例として出されているものです。SASの元記事だと解析コードなどが一部省略している部分されているので、ここでは一つ一つ丁寧に見ていこうと思います。（個人的学習の意もこめて）
使用データ
解析で使用するデータセットは、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例を解析対象とします。本当は選択バイアスの存在とか考慮しないといけないんだろうとは思ってはいます。
1566例のうち、曝露群は403, 非曝露群は1163で、マッチング方法はOptimal Variable Ratio Matchingとして解析を行います。なお傾向スコアの算出にはロジスティック回帰分析を用います。
解析結果
マッチングさせた両群のサイズ、傾向スコアに関してはこんな感じ。
曝露群のサイズに合わせた1:1マッチングが行われています。
また、各説明変数に関して両群での分布は下図。
また、PlOTオプションでは各連続変数の累積分布を作成しますが、それを見るとマッチングさせる前とさせた後での共変量の分布の差が分かります。下記はAge（年齢）という変数についてのプロットですが、マッチングさせる前は喫煙した群のほうが年齢が高い傾向がありましたが、マッチング後はそれが改善しています。（ほかの変数についてはぜひコードを実行してみてください）
Standardized Mean Differences（標準化された平均の差） plotでは、曝露群と、非曝露群での変数の標準化した平均値の差をプロットしています。縦軸が各変数で、横軸が標準化後の値です。
また、マッチングを行うことで両群の各変数に対する分散の比がより1に近づき、改善される傾向がうかがえます。
調整前の単純な体重変化の群間差は下表のとおり約2.54kgです。
次にここまでやってきた傾向スコアマッチングで推定された治療群における平均因果効果（ATT）は3.2792kgと、交絡因子で調整していない場合の2.54kgと比較して大きく増加しています。
今回の分析の補足、拡張についてはまた時間を見つけてまとめようと思います。また、分析よりもモデルに組み込む変数の選択やバイアスの低減のための作業が重要なので、なかなか観察研究は奥が深いなと思っています。
サンプルコード
データセット自体の作成コードはかなりの量になるので、使用データの項目よりSAS Helpのページに飛んで利用してください。
/*データセットの加工*/
data SmokingWeight;
set SmokingWeight;
if nmiss(of _numeric_) > 0 then delete;
id = _N_;
run;
data SmokingNoResp;
set SmokingWeight;
drop Change Weight;
run;
/*マッチング後の共変量のバランス確認*/
ods graphics on;
proc psmatch data=SmokingNoResp;
class Sex Race Education Exercise Activity Quit ;
psmodel Quit(Treated='1') = Sex Age Education Exercise Activity YearsSmoke PerDay;
match method=varratio(kmin=1 kmax=4) distance=lps caliper=0.5;
assess lps var=(age YearsSmoke BaseWeight PerDay) /plots=(CDFPlot BoxPlot StdDiff);
output out(obs=all)= SmokeMatched1 weight=matchattwgt matchid=MID;
run;
/*調整前の群間差*/
proc ttest data= SmokingWeight;
class Quit;
var Change;
run;
/*マッチング後のATTの推定*/
data SmokeAnalysis;
merge SmokeMatched1 SmokingWeight;
by id;
run;
proc ttest data=SmokeAnalysis;
class Quit;
var Change;
weight matchattwgt;
run;
参考
Tools for Assessing a Propensity Score Model in SAS/STAT
https://www.sas.com/content/dam/SAS/support/en/sas-global-forum-proceedings/2019/3056-2019.pdf