【EZR】傾向スコアによる逆数重み推定(IPW)の解析手順【図解】

記事の要約：EZRで傾向スコアによる逆数重み推定(IPW)の解析手順【図解】

この記事で学べること

• EZRによる傾向スコアの計算方法
• 傾向スコアを用いた逆数重み推定法（IPW）とは？
• 傾向スコアを用いた逆数重み推定量の計算方法
• 重みづけをしたあとの背景情報のバランスの確認方法

この記事のまとめ

• 傾向スコアは背景情報から予測された処置をうける確率
• 傾向スコアを用いた逆数重み法を使えばグループ間の背景が揃う
• EZRを使えば簡単に傾向スコアの計算と重みづけ推定ができる
• EZRや因果推論に関する書籍は以下がおすすめ
  • （因果推論をRで学びたい方向け）統計的因果推論の理論と実装 (Wonderful R)
  • （因果推論をPythonで学びたい方向け）つくりながら学ぶ! Pythonによる因果分析: 因果推論・因果探索の実践入門
  • （因果推論の理論も学びたい方向け）調査観察データの統計科学: 因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学) 

傾向スコアとは

傾向スコアに関する基本的な説明はこちらの記事で解説しています．この記事でもざっくりと要点だけを復習しておきます．

1. 比較するグループ間でもともとの背景が違うと比較結果は信頼できなくなる（交絡）
2. グループ間の背景情報のバランスをとるために傾向スコアを用いた手法が利用できる
3. 傾向スコアは「背景情報から予測される処置をうける確率」として計算できる
4. 傾向スコアを用いて交絡を取り除く方法はマッチング・逆数重み法（IPW）などがある

つまりグループ間の比較をフェアにするのであれば，グループ間の状態を揃えてあげる必要があり，そのときに傾向スコアはとても便利だよ，ということです．

こちらの記事では傾向スコアを使ったマッチングの解析手順について説明していました．今回は良く使われる傾向スコアの利用方法2ということで，傾向スコアによる逆数重み推定法（IPW; Inverse Probability Weighted estimator）について紹介します．

傾向スコアを用いた逆数重み推定法（IPW）とは

傾向スコアを用いた逆数重み法とは，その名前の通り，各対象者のデータを傾向スコア（厳密には所属するグループに割り当てられた確率）の逆数で重みづけしてあげる方法です．重みづけするというのは「水増しする」というイメージをもってもらえば分かり易いかと思います．例えば重みが「２」であれば，その対象者のデータは2人ぶんとして考える，といったイメージです．

さて，それでは傾向スコアの逆数で重みづけるとはどういうことでしょうか？図でみてみましょう．まず，以下のように，より具合の悪い人が処置をうけているような例を考えてみます．このようなデータでは，もともとの状態がグループ間で異なるためにフェアな比較ができず，処置群に不利な結果がでてしまいますよね．

それでは傾向スコアを利用した逆数重みを適用してみましょう．基本的なルールとして，以下の操作を覚えておいて下さい．

• 処置群：処置をうけた確率の逆数（=傾向スコア の逆数）で重みづけ
• 非処置群：処置をうけなかった確率の逆数（=1-傾向スコア の逆数）で重みづけ

まずは「処置群」について考えます．処置群では，それぞれの対象者のもつ傾向スコアの逆数をかけた分だけ，人数を水増しします．例えば傾向スコアが1/3=33%の場合，その逆数は3/1なので，人数を3倍してあげればOKです．他の対象者についても同様に計算して，ここでは水増しした人数を赤色で表示しておきます．

次に非処置群について考えてみます．非処置群では「処置をうけない確率」の逆数で水増しします．処置をうけない確率は，1-傾向スコアで計算できます．例えば傾向スコアが1/3=33%の場合，1-1/3=2/3であるため，その逆数は3/2となります．従って人数を3/2倍してあげればOKです．他の対象者についても同様に計算して，ここでは水増しした人数を緑色で表示しておきます．

これで重みづけは完了です！重みづけ前後で，対象者の分布がどのように変わったか確認してみましょう．

どうでしょうか？重みづけをした後は，両端の極端な傾向スコアをもつ対象者以外の人数は均等になりましたね！この状態であれば，重みづけ前の状態よりも，グループ間の比較をフェアに行うことができそうですね．

あとは両端の極端な傾向スコアをもつ対象者についてですが，これはもっと端的にいうと「相手のグループに類似した傾向スコアをもつ人がいない」対象者たちです．傾向スコアマッチングでは，マッチング相手が見つからず，マッチングデータに抽出されなかった人たちですね．この人たちの扱い方については種々方法がありますが，比較するグループ間で重ならない範囲の傾向スコアをもつ対象者は重みづけ前に解析から除外する，あるいは傾向スコアを5%や95%などの値に丸め込む等，この辺りも状況によって色々な手法が提案されています．今回も，そのような両端の対象者を切り捨てて解析を行えば，両グループ間の対象者の分布が一致することが分かりますね．

傾向スコアを利用した逆数重み法の重要なコンセプトは，ここまで理解できればまず大丈夫です．あとは，この重みづけした解析をどうおこなうかですが，ポイントが2つあります．

• グループ間の比較を重みづけして行う（実際にデータが水増しされるわけではない）
• 重みづけした推定量は標準誤差の推定に工夫が必要

この辺りの実際の解析の操作方法については，以下のEZRでの解析手順で説明していきます．ただEZRだけでは，重みの計算まではできても，実際にその重みを適用した解析自体は実行できません．そのためEZRのスクリプト画面に，Rのコードを追記する必要がありますので，ちょっと上級者向けになります．

事前準備：データのダウンロードと読み込み

それでは早速データセットを読みこんで解析していきましょう．

データセットは練習用のCSVデータを準備しましたので，手元にちょうどよいデータがなければ，これらをダウンロードしてご利用ください．

準備１．参考データ（CSV）をダウンロード

練習用CSVデータをダウンロード

1 ファイル 42.61 KB

ダウンロード

※後ほどEZRで読み込みたいデータの場所を指定する必要があるため，デスクトップなど自分が分かりやすい場所に保存しておくことをおすすめします！

準備２．EZRを起動 ⇒ CSVファイルの読み込み

CSV・ExcelファイルのEZRにおける読みこみ方は，以下で詳しく説明していますので，参考にしてみてください．今回はデータセット名はデフォルトのまま Dataset と名前をつけて読みこんでいます．

【EZR】CSV・Excelデータの読み込み【図解】

EZRでデータ解析をしてみたいけど，データの読みこみからつまづいていませんか？ この記事ではEZRを利用したCSV・Excelデータの簡単な読み込み方法を紹介しています． 実際の実行画面をもちいて図解していますので，初心者の方も安心して取り組めます． 早速EZRでのデータの読みこみ方を学びましょう！

robustwife.com

2021.01.16

傾向スコアを用いた逆数重み推定の手順

傾向スコアを用いた逆数重み推定を行う手順は大きく３つに別れます．

1. 傾向スコア＆各群に所属する確率の逆数（重み）を計算する
2. 重み付け前後で背景情報のバランスを確認する
3. 重み付けを考慮しながらグループ間の比較を行う

1．傾向スコア＆各群に所属する確率の逆数（重み）を計算する

まずは傾向スコアを推定してみましょう．このステップは基本的に傾向スコアマッチングの記事で紹介した内容と同じものになりますが，傾向スコアを計算すると同時に，各群に所属する確率の逆数（重み）を計算します．

また傾向スコアマッチングの記事で紹介したように，今回も傾向スコアはロジスティック回帰モデルを使って推定してみましょう．

今回は処置変数である treat_bin (1:処置あり, 0:処置なし) に関する傾向スコア（treat_bin=1となる予測確率）を背景情報である cnf1 と cnf2 から計算してみます．

1. 名義変数の解析
2. 二値変数に対する多変量解析（ロジスティック回帰）
3. 目的変数 欄にカーソルをあてて処置変数 treat_bin をダブルクリック（※1）
4. 説明変数 欄にカーソルをあてて交絡因子 cnf1 と cnf2 をそれぞれダブルクリック（※2）
5. □傾向スコア変数を自動作成する にチェックを入れる
6. □逆数重みづけ(IPW)のデータセットを作成する にチェックを入れる
7. OK

※1．画面を開いた時点ですでに「目的変数」欄にカーソルが当たっているはずです
※2．treat_bin をダブルクリックしたら自動的に「説明変数」欄にカーソルが移動します

これで各対象者について，背景情報 cnf1 と cnf2 から計算された傾向スコアと逆数重みがデータDatasetに保存されているはずです．早速データの中身を以下手順で見てみましょう．

データの右端に「PropensityScore.GLM.1」という新しい変数ができていますね．これがcnf1とcnf2から予測した各研究対象者において処置が行われる確率であり，傾向スコアになります．

また「weight.ATE.GLM.1」という新しい変数も同時に作成されています．これは各グループに所属する確率の逆数となりますので，今回逆数重み法で利用する重みになります．

２．重み付け前後で背景情報のバランスを確認する

それでは重みづけ前後で，背景情報が揃ったかどうかを確認してみましょう．背景の揃い方の確認する方法も色々とありますが，今回は背景表の計算時に Standardized Mean Difference (SMD) を使ってみます．SMDは０に近いほどグループ間の背景情報の差が小さいと考えることができます．

傾向スコアマッチングの場合は，マッチング前のデータセットと，マッチング後のデータセットでそれぞれ背景表を作成し，同時にSMDを計算するだけでOKでした．しかし重みづけした背景表はEZRでは簡単に計算できないので，以下のようなステップを実行する必要があります．

1. Datasetを使って重みづけしない背景表を作成する
2. Rスクリプトを追記して重みづけデータを作成（背景表でのみ利用します）
3. ###サンプルの背景データのサマリー表の出力### 以降を選択ハイライトして「実行」

Datasetを使って重みづけしない背景表を作成する

1. グラフと表
2. サンプルの背景データのサマリー表の出力
3. [群別する変数] treat_bin, [カテゴリー変数] cnf2, [連続変数(正規分布)] cnf1
4. [カテゴリー変数の検定方法] 独立性のカイ二乗検定, [標準化差を表示する] Yes
5. OK

ここで得られる背景表は重みづけ前のものなので，当然両グループでバランスがとれていません．得られた表は後ほど，重みづけしたものと並べて眺めてみることにしましょう．

Rスクリプトを追記して重みづけデータを作成（背景表でのみ利用します）

ここでは重みづけをした背景表およびSMDを計算してみます．EZRには重みづけした背景表を自動で作成することができません．そこで，上記で作成した重みなしの背景表のRスクリプトに少し手を加えて，重みづけした背景表を作成してみたいと思います．

1. ###サンプルの背景データのサマリー表の出力### の直下に以下のRコードをペースト
2. 以下の4か所に修正を加える
   • Dataset → svyDataset （※2か所！）
   • CreateCatTable → svyCreateCatTable
   • CreateContTable → svyCreateContTable

# 追記部分 --------------------------
install.packages("survey") # オンラインで一度だけ実行
library(survey)
svyDataset = svydesign(ids = ~1, data = Dataset, weights = ~weight.ATE.GLM.1)
# -----------------------------------

###サンプルの背景データのサマリー表の出力### 以降を選択ハイライトして「実行」

それでは上記の追記をしたコードを実行するために「###サンプルの背景データのサマリー表の出力###」意向をマウスで選択（青くハイライトされます）して，実行ボタンをクリックしてみましょう．

※ install.package(“survey”) 行の実行はオンライン環境が必要です．ここではパッケージのインストールをしています．インストール元を問われたら（国名などを選ぶ画面がでてきたら），0-Cloudを選んで進めてもらえば問題ありません．

作成した背景表を重みづけ前後で並べてみましょう．重みづけ前はSMDが0.2付近ですが，重みづけ後は非常に小さくなっていますね．実際にそれぞれの背景情報の値をみてみても，重みづけ後は各要素の割合や平均値のバランスが良くなったように見えます．これでグループ間での背景情報の不均衡を取り除くことができました！

３．重み付けを考慮しながらグループ間の比較を行う

重みづけによって背景も揃えることができたので，やっとこれでグループ間の結果変数をフェアに比べることができそうですね．次のステップでは，グループ間で結果変数を比べてみます．この時に注意しなければいけないことが，

1. 処置効果の推定は重みづけを考慮して行う
2. 一般化推定方程式などで標準誤差を推定する

という2点です．（１）については当然ですね．（２）については少し難しくなってくるので読み飛ばしてもらっても大丈夫ですが，IPW推定を行う際に単純な線形回帰モデルなどを利用してしまうと，重みさえ考慮していれば処置効果自体は問題なく推定できるのですが，推定値のばらつき具合を示す標準誤差という指標が過小評価されてしまいます．その点を考慮しつつ解析を行うことができる推定手法として，一般化推定方程式（Generalized Estimating Equation: GEE）があります．

GEEを利用した解析はEZRではできませんので，以下のコードをRスクリプト画面に張り付けて実行する必要があります．以下では，結果変数である outcome_num を処置変数である treat_bin 間で比較する際に，処置確率の逆数である weight.ATE.GLM.1 で重みづけしています．

install.packages("geepack") # 一番最初に使うときだけ実行
library(geepack)

fit = geeglm(
  outcome_num ~ treat_bin, 
  data = Dataset, 
  weights = weight.ATE.GLM.1, 
  id = id,
  family = gaussian(link = "identity"), 
  corstr = "independence"
  )

# 結果の出力
summary(fit)

上記について簡単にポイントを説明しておきます．

• outcome_num ~ treat_bin: 結果変数 ~ 処置変数
• data = Dataset: 利用しているデータセット名を設定
• weights = weigth.ATE.GLM.1: 重み
• id = id: 各対象者を示す連番などの変数（各症例にユニークなものをあてること）
• family = gaussian(link = “identity”): 分布族．解析によって変える．
• corstr = “independence”: 相関行列．近接法で推定していて頑健なので気にしなくて良い．

family（分布族）については解析に利用する結果変数によって変える必要があります．生存時間分戦機を行う場合などは特殊なので，ここでは連続の場合と，2値の場合だけを出しておきます．

解析やデータに応じて，それぞれの変数の名前を書き換える必要があります．もし違うデータや変数を利用する場合には，データ内の変数名とコード内の変数名が完全に一致しているかどうか，しっかりと確認するようにしてみてください．

結果としては，重みづけをした後の処置効果の推定値は 54.439 となり，処置によって結果変数であるoutcome_numは54.498だけ上昇することが分かりました．この結果については，既にcnf1やcnf2による交絡の影響は受けていないと考えることができます．

これでIPWの解析は終了です．EZRの機能だけではカバーしきれないところがあるので，すこし作業が煩雑になりましたが，慣れれば難しいものではありませんので，ぜひご自身のデータでもチャレンジしてみてください．

最後にひとつだけコメントを添えておくと，今回IPWで推定した処置効果は，厳密には「全ての対象者における平均的な処置効果（Average Treatment Effect: ATE）」と呼ばれるものです．本来データには各対象者は処置を受けた・受けていないのどちらか一方のグループに所属していますが，仮想的に「対象者全員が処置をうけた場合・全員が処置をうけなかった場合」を考えた時にどれくらいの効果があったのかを考えているということです．同じ対象者同士でくらべているので，背景も当然そろいますよね．これがIPWを含む因果推論手法における反実仮想というコンセプトになります．

これに対して「実際に処置をうけた対象者における平均的な処置効果（Average Treatment effect on the Treated; ATT）」という考え方もあります．関心があるのは実際に処置をうけた人たちであり，その人たちがもし処置を受けていなかった場合と比べると，どのくらい効果があったのだろうか？ということを考えるということですね．このあたりはまた面白いところなので，また別の記事で紹介します．

今回の内容を通して，以下の論文が非常によくまとまっていて分かり易いので，良ければ参考にしてみてください．

Austin PC, Stuart EA. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. Stat Med. 2015 Dec 10;34(28):3661-79. doi: 10.1002/sim.6607. Epub 2015 Aug 3. PMID: 26238958; PMCID: PMC4626409.

まとめ

今回はEZRを用いた傾向による逆数重み推定法（IPW）の解析手法に関する一連の流れについて紹介しました．EZRを使えば難しい操作なく傾向スコアの推定と重みの計算が可能です．

重みづけを行った後の処置グループ間の背景を比べたり，重みづけをしながらグループ間で結果変数を比較する際には，Rのコードを書き足す必要がありましたが，この記事の解説通り進めてもらえば問題ありません！

またEZRで使える因果推論系の機能も日々アップデートされているので，今後も楽しみですね．

EZR開発者の先生が出されている書籍や，因果推論についてもっと専門的なところまで勉強したいという方にお勧めの書籍も以下に紹介しておきます．

• （因果推論をRで学びたい方向け）統計的因果推論の理論と実装 (Wonderful R)
• （因果推論をPythonで学びたい方向け）つくりながら学ぶ! Pythonによる因果分析: 因果推論・因果探索の実践入門
• （因果推論の理論も学びたい方向け）調査観察データの統計科学: 因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)