Stata分析機能例題集
傾向スコア・マッチング（プロペンシティ・スコア・マッチング）

データの準備

早速下記図1のようにコマンドを入力します。

cattaneo2には新生児の体重bweightと、その母親の妊娠期間中の喫煙習慣mbsmokeを始めとする様々なデータが入っています。
傾向スコア分析を実行する前に、少しデータについて考察の為、ヒストグラムを作図しましょう。

以下のグラフが表示されます。（図2）

2つとも似たような分布を示していますが、これを記述統計量で見てみましょう。

meanの項目から新生児の体重には300g程度の差があることが分かります。
ここで、データファイルcattaneo2に含まれる母親の数を確認します。

全体では4500程度の母親がいて、妊娠期間中に喫煙習慣のあった母親は864人であることが分かります。
喫煙習慣のなかった3,778人に比べると、かなり少ないことが分かります。

話を新生児の体重に戻します。
例えば、喫煙習慣の有無によって分けた2つのグループの新生児体重に有意差があるか、対応のない2群のt検定（等分散を仮定）を実行してみましょう。

対立仮説はnonsmoke-smokerが0より大きいとして良いと思います。
平均体重の差は約275gで、新生児体重に有意差ありとなりました。

次に、傾向スコア分析を行う際の問題意識を確認します。
以下のコマンドを実行して、2つのグループの母親について基本的な統計量を調べてみましょう。

喫煙習慣という「処置」が、新生児の体重という「アウトカム」に与える影響について考える場合、新生児の体重に影響を及ぼすであろう母親に関する変数の分布には「差がない状態」が理想的です。
つまり、分布が同じなら交絡を考えなく済みます。
喫煙習慣の無い母親のグループを見ると、既婚者である割合が0.75で、喫煙者の0.47よりも、その大きさが目立ちます。
fbabyは第一子である事を示すダミー変数です。
喫煙習慣のない母親のグループでは第一子である割合が0.45で、喫煙者グループの0.37よりも高いことが分かります。
次に示すマッチングを行うことの意味は、あたかも処置のランダム割付ができ、交絡を考えなくても良い状態になる所にあります。

マッチングの考え方

カテゴリカルなmmarried（既婚者は1）やfbaby（第一子は1）は、相手のグループから同じ値の人を探すのは簡単ですが、mage（年齢）、medu（母親の学歴）などの連続変数が入ってくると、一番近い人を探すのための計算量はかなり大きくなります。
つまり、変数が増えるほど、似たような人を探すのには手間がかかります。
そこで傾向スコア分析では、処置の割り当て（mbsmoke）を被説明変数とするロジットモデルを推定し、そこからゼロになる確率（傾向スコア）を求め、この傾向スコアをマッチングさせるというアプローチを採用します。
このようにすれば、ロジットモデルの説明変数について個別にマッチングを行うよりも、はるかに計算量は少なくなるという考え方です。
そこで、確認の意味でロジットモデルを別個に推定してみましょう。

Stataではモデル推定後、そのモデル式の理論値をpredictコマンドで求めることができます。
実際に次のコマンドを実行して、0を取る確率p_hat0を計算してみましょう。

傾向スコアは基本的にゼロになる確率のことを指しますが、本来の目的は変量の情報を集約することにありますので、1となる確率を利用したり、正規分布に近づけるために、さらに変換する場合もあります。
このようにして作成した傾向スコアの分布を確認してみましょう。

この図はマッチング前の状態ですので、それぞれnonsmoke=3,778, smoker=864から作成したものです。

ATEを求める

傾向スコア分析は相手の群から傾向スコアが一番近い人を探し、アウトカムの差を取って、最後に差の総平均（Average Treatment Effect）を求めるというものです。
ここで、目的の値を一気に求めてみましょう。

先頭の2つの単語、teffects psmatchが一つのコマンドとなっています。
半角スペースを挟んで（bweight）はATEを求めるアウトカム、処置の割り当てを示すロジットモデルは前出のmbsmoke mmarried c.mage##c.mage fbaby meduとします。

傾向スコア分析を利用した結果、喫煙習慣（Treatment）により新生児の体重は210g減るという結果になりました。
対応の無い2群のt検定とは、かなり異なる値になりました。

反実仮想とは

ここからは、StataのPDFマニュアルを用いて傾向スコア分析の反実仮想（counterfactual）の考え方を簡潔に説明します。
最初にAさんが処置を受けた時のアウトカムをY1、処置を受けていない時のそれをY0とします。
実際には、過去に戻って喫煙習慣を変更して、もう一度、子供を出産することは不可能ですが、似た人を探しだしてくることで、あたかもこれを行っているかのごとく考えます。
複数の人について、この（Y1-Y0）の差の平均を求めたものがATEとなります。

これは観測したデータそのものです。
血圧降下剤を服用した人（Treated）と、そうでない人のアウトカム（血圧）の平均は同じように見えます。
つまり、血圧降下剤を服用しても効き目がない、ということになります。
しかし、その2群で体重の分布には差がありそうです。
乱暴な言い方ですが、すごく痩せた人と、すごく太った人では、薬の効果も違うように思えます。

そこで、傾向スコア分析を行って、相手の群から傾向スコアの最も近しい人を仮想の自分として対にします。

図で、中空の点は反実仮想でマッチングさせたデータに相当します。
この時、降圧剤を服用したグループの血圧平均は137で、服用していない群では182となりました。
かなり、大きな差があることが分かりました。

メニュー操作

今回、初めてStataを操作した方のために、メニュー操作とダイアログボックスを使って操作する方法をご紹介しておきます。
メニュー操作で表6の結果を出してみましょう。

メインメニューで統計/処置効果/アウトカム（連続）/傾向スコアマッチングを選択します。

先にコマンド入力した情報に相当するものを図のように入力します。
ロジットモデルの独立（説明）変数についてStata独自の記法がありますので、これを説明しておきます。

とありますが、これは次のようにも書けます。

つまり、mageの一乗項と二乗項を一緒に説明変数として利用する場合は、##を利用することで、少し簡単に書けということです。
さらに、c.mageのcの役割ですが、これは交差項が二乗項を作成する場合、変数が整数型だとStataは自動的に「この変数はカテゴリ変数である」と理解してダミー変数を作成します。
今、利用している母親の年齢mageは整数ではありますが、カテゴリ変数ではありませんので、明示的に連続変数である事を示すc.を変数名の前に付けなければなりません。
このようにStataにはモデルを記述する際に注意すべき事柄がいくつかありますので、それは折に触れてこの連載の中で紹介することにします。
既にお分かりのように、交絡があるときのアウトカムの大きさを評価する場合は、傾向スコア分析を利用すれば、処置効果を正しく推定できます。
ただし、ロジットモデルの定式化には十分、時間をかけてください。
また、「キャリパ」というものを利用することも考慮しなければなりません。

先の章で平均処置効果（ATE）を求めてみました。
実際に用いたコマンドは次のようなものでした。

分析のモチベーションはもちろん、妊娠期間中の喫煙習慣（mbsmoke）が新生児体重（bwight）に与える単独の影響を調べる事です。
そして、交絡要因として母親が既婚（mmarried）であるか、母親の年齢（mage）、妊娠第一期（prenatal1）、第一子（fbaby）を利用しました。
結果として、喫煙習慣のある母親の新生児は平均で約235kg小さいことが分かりました。
今回のコマンドで最後にオプションのgen()を利用しています。
これは後段の分析をより速く行うためのお作法ですので、気にせず進みましょう。

さて、今回のテーマは推定後に必ず実行しなければならない、共変量に関するバランスチェックの機能を体験する事です。
それでは最初にマッチング前の原データについて、共変量の情報を調べてみましょう。
teffects psmatchというコマンドの場合、共変量のバランスチェックには3つのコマンドが用意されていますので、順番に見ていくことにしましょう。

オプションbaselineを利用すると、マッチング前の状態を調べることができます。
二群の平均と分散を比較していますが、もちろん、平均と分散が近いほど分析者にとってうれしい状態です。
しかし、ここでは平均、分散ともにバラツキが見てとれます。

それでは、次にオプションを外してマッチング後の状態について見てみることにします。

マッチング後の平均は元データのそれに比べ0に近くなり、分散は1に近くなっています。
この表作成に用いたStandardized differenceとはどのようなものなのか、気になるところかと思いますが、それはまた、時間のある時にStataのオンラインマニュアル（PDF）を調べてみてください。

次はこれを図（密度分布）で確認しましょう。

マッチングしたmageの分布曲線（図1の右側）がほぼ重畳していることが分かります。
先の表3や、この図のマッチング後のズレが目立っている場合、ロジットモデルの共変量を他のものに入れ替える必要があります。

このチェックを箱ひげ図で行うコマンドもあります。

もちろん、右側の図にあるように、マッチングした時のnonsmokerとsmokerの分布がほぼ同じ状態になっていることが理想的な状態です。

さて、実は手元のサンプルデータは喫煙習慣の有無を示すデータmbsmokeの他に、具体的に何本吸うのかというデータmsmokeも用意されています。
データを順番に確認してみましょう。
最初に喫煙習慣の有無mbsmokeを確認します。

次にその内訳データmsmokeを確認します。

1～5本の人がむしろ少なめで、6～10本、11本以上の人が300人以上います。
コマンドをteffects ipwに変更すると、これらのグループにおける喫煙習慣のATEを調べることもできます。

POmean（Potential-Outcome means）とあるのは、反実仮想の考え方を利用したときの、喫煙習慣の無い人たちの新生児体重の平均（3402g）です。
これを基準にして、ATEは5本以下の場合は-179g、6～10本の場合は-259g、11本以上だと-237gとなっています。
奇妙な結果になっていることが分かります。
すなわち、中途半端（5～10本）に吸うよりも、11本以上吸ったほうが、まだましだ、という結果になっています。

バランスチェックを行ってみましょう。

表7に示すように母親の年齢はマッチングさせても分散に無視できない差があります。
このような結果になったら、割り付けのモデル（ロジットモデル）の定式化を再検討しましょう。
そして、推定後の結果に不合理な結果が出なくなるように、試行錯誤を繰り返しましょう。
これは前出のteffects psmatchコマンドを利用した時でも全く同じです。
これでお分かりいただいたようにATEが有意であるから、それでよし、というものではありません。
共変量のバランスチェックが大事であり、それを無視してしまうと、どのような結果になるのか、ご理解いただけたと思います。
今回はteffects psmatchの他に、teffects ipwというコマンドを利用してATEを推定しました。
StataにはATEを推定するコマンドとして次に示すようなものが用意されています。

前の章では、teffects psmatchだけでしたが、今回はteffects ipwも利用しました。
teffects psmatchはロジットモデルにより求めた傾向スコアを利用してマッチングするものでした。
ipwはどんな理屈でマッチングさせるのか、最後に簡単に触れておきましょう。

図3は前回も示した概念図です。
中空のデータをロジットモデルを用いて算出した傾向スコアによるマッチングで選び出します。
一方、teffects ipwはロジットモデルから得た理論値の逆数を利用して図4のようなデータを作成し、ATEを求めます。
基本的な発想としては反実仮想の考えに基づいていますが、実際の計算で利用しているデータは実データを重み付けしたものです。

よって、ATEの標準誤差の計算方法には特別なものは利用していません。
言い換えると、図3の中空のマッチングデータを利用するが故にteffects psmatchによるATEの標準誤差の計算は難しく、その点においてAbadie and Imbens（2011）はとても貢献度の高いものと言えます。
実際、Stataでも以前からpsmatch2というadoファイルを利用すれば、ATEの計算だけはできましたが、その結果画面をみると標準誤差は空白になっていました（つまり、有意性の検定はできなかった）。

次は Abadie and Imbens（2011）の貢献と題して、彼らの論文のアイデアをご紹介します。

今回は処置効果推定量ATEの推定値とその標準誤差の計算方法を中心にStataの機能をご紹介します。
少しだけ数式を用いて説明する箇所がありますが、他の統計解析ソフトとの違いを示す意味で、重要な箇所ですのでご容赦ください。
傾向スコア分析とは平均処置効果（Average Treatment Effect）を推定することですが、毎回申し上げますように、交絡要因の影響を排除して処置効果の平均値を推定することを目的とします。
Rosebaum and Rubin（1983）はこの目的を達成するために傾向スコアを利用することを提案しました。
共変量Xを個別にマッチングさせるのではなく、傾向スコアだけをマッチングさせることで、いわゆる、“次元の呪い”から解放されます。
例えば、血圧降下剤の薬効を測定する場合、体重や身長などのベースライン特性を個々にマッチングさせるのは大変です。
Rosebaum and Rubin（1983）は共変量を用いたロジットモデルから推定した傾向スコアだけをマッチングさせることで交絡の影響を排除できることを示しました。

これはRubin（1974）の考えた反実仮想に基づくATEの定義式です。

ある人が処置を受けた時のアウトカムをY(1)、そして同じ人が仮に処置なしであった場合のアウトカムをY(0)とします。
処置の割付をWで表現すると、アウトカムと処置の関係は次のようになります。

τを求めるにあたりRosebaum and Rubin（1983）は3つの仮定を設定しています。
傾向スコア分析に興味のある方なら、その中の一つである、“強く無視できる割り当て（strong ignorability）”と呼ばれる仮定を一度は耳にしたことがあるのではないでしょうか。
この仮定と他の計算式を用いることで、次に示す(2)式により観測可能なデータから(1)式が計算できます。

平均処置効果

上記は共変量Xが同じで処置がある時のYの条件付期待値から、処置がない時のYの条件付期待値を引くことで、目的の処置効果を推定できることを示しています。
Abadie and Imbens（2012）は (2)式の具体的な計算式として次式を提案しています。

ここでˆθは傾向スコアを推定するロジット（またはプロビット）モデルのパラメータです。
ポイントはマッチングセットを示すJM(i, ˆθ)です。
Abadie and Imbens（2012）はマッチングセットにおける標本の“復元あり”を前提としており、一度、マッチングに利用した標本はそれ以降、何度でも利用可能な設定を採用しています。
Mはマッチングに利用する標本の数で、通常は1とします。
例えば、M=2の場合、片方の群のAさんに対し、もう一方の群からは傾向スコアが一番近いBさんと二番目に近いCさんのアウトカムを平均し、それをAさんのYiから引きます。
M=1であるにも関わらず、Aさんに対してBさんとCさんの傾向スコアがまったく同じ場合も同じ要領で計算します。

実際にStataで試してみましょう。

オプションgen()を利用することでマッチング情報データセットにシリーズid（1から74）を作成しました。
74はマッチングした時に同点の人が74人いることを示しています。
例えば、データエディタを開いて一行目の情報を見てゆくと、右端のほうに次のような情報があります。

データエディタ1行目の母親は変数mbsmokeの情報から非喫煙者（nonsmoker）であることが分かります。
この母親の傾向スコアにマッチングする喫煙習慣のある母親は7人いることが分かります。
2行目の数字は、各母親のデータエディタ中の行番号です。
次のコマンドを実行して喫煙習慣があることを確認してみましょう。
nはデータエディタ上での行番号を意味します。

この7人のアウトカムYjの平均を求め、Y1から引いた個別の処置効果はいくつでしょうか。
それは次のコマンドで簡単に調べることができます。

teffects psmatchコマンドでATEを推定した後で、predictコマンドをオプションなしで実行すると、teに各人の処置効果を格納します。
データエディタ画面を表示して実際に1行目の母親に対する値を調べてみると、−337.14である事が分かります。

Abadie and Imbens の標準誤差

ここでは推定結果の標準誤差に注目します。
Rosebaum and Rubin（1983）以降、τˆNの分散推定量の計算にブートストラップ法を利用する方法も提案されましたが、それには不都合な点があることが指摘されていました。
実際、Stataの古いバージョン（v12以前）でも傾向スコア分析用のadoファイルpsmatch2を利用すれば、ATEの推定は可能でしたが、このadoファイルでは敢えてATEの標準誤差は推定しませんでした。
試しにpsmatch2を検索し、インストールしてみましょう。

キーワード検索の結果を画面に表示します。
画面中央のやや下側に表示されるPSMATCH2のリンクをクリックし、個別のオンラインヘルプページにジャンプしたら、click here to installのリンクをクリックしてインストールを実行します。
インストールが完了したら、改めてデータを読み込んで、ATEを推定してみましょう。

結果を見れば明らかなように、psmatch2ではS.E.（標準誤差）の項目は算出されません。
このような理由でStata12まではATEの推定は正規なコマンドとしてはサポートしていませんでしたが、Abadie and Imbens（2012）を採用することで、Stata13から正式なコマンドとして内蔵されることになりました。
傾向スコア分析を実行できるソフトウェアはStata以外にもあると思いますが、その標準誤差の計算に、どのような理論を利用しているのか、注意を払う必要がありそうです。
ここでは Abadie and Imbens（2012）に示されているATEと標準誤差の計算式を下記式に示します。

ポイントは式の右辺第二項が分散を小さくする方向に修正を加えている点です。
この中に含まれるcやIまで細かく説明することは避けますが、理屈としては推定した傾向スコアを用いると、その分散は小さくなる方に修正されるということです。
Rosebaum and Rubin（1983）の設定した3つの仮定が成り立ち、もし、我々が標本から真の傾向スコアが入手できて、それでATEの推定量τˆ∗Nが分かるとすれば、統計学の理論を使って次の関係になることが分かります。

しかし、残念ながら、傾向スコアはあくまで推定値ですから上記式のように単純には行きません。
そこで、Abadie and Imbense（2012）は次のような関係になることを示しました。

Rosebaum and Rubin（1983）以来、数多くの提案がなされ、傾向スコアによる解析手法は徐々に進歩しています。
ここでStataに興味のあるユーザのために次の二点を強調しておきたいと思います。

• τˆNの計算において傾向スコアが同点の場合の処理は(3)式を採用している。
  よって、乱数を使って計算に利用する標本を決める方法は利用していない。
• τˆNの標準誤差の計算には(6)式を利用している。
  これによりτˆNの有意性検定が行えるようになった。

回帰調整

このように見てきますと、傾向スコアteffects psmatchコマンドを利用することが必須であるかのように感じあるかもしれませんが、傾向スコア分析といった場合、他にも第二回でご紹介したteffects ipwとteffects ra（Regression Adjustment ： 回帰調整）というものがあります。
この連載の最後にteffects raについて基本的な考え方をご紹介しておきます。
StataのPDFマニュアル[TE]にある次の図でまずはイメージをつかみましょう。

この図4の塗りつぶされている点は実際に観測したデータで、中空で同色のものは潜在的なアウトカムを示します。
二群に対してそれぞれ回帰直線をフィットし、その理論値を全てのデータについて求め、それらの平均を二群で計算します。
最後に両者の差を取り、ATEとします。
まずは、teffects raコマンドで一気に推定を実行し、推定結果をraという名前で保存します。

コマンドの後ろにはアウトカムモデル（被説明変数をbweightとする回帰モデル）を入力します。
prenatal1は妊娠第一期に来院した場合は1、mmarriedは母親が既婚者の場合に1、mageは母親の年齢、そしてfbabyは第一子である場合に1を取るダミー変数です。
最後のmbsmokeは処置を示します。

ATEは−239となっています。
この計算の舞台裏を見てみましょう。
図4にあるように、喫煙習慣の有無により分けた2つのグループを考えて回帰直線をフィットします。
最初に喫煙習慣の無い母親のデータだけを利用してモデルを推定し、次に、すべての母親の潜在的なアウトカムとしてy0を算出します。
そしてy0の平均をsuコマンドで調べ、平均値をスカラーy0mに格納します。

次は喫煙習慣のある母親について同じことを行います。
潜在的なアウトカムを計算するときは全員のものを求めますのでご注意ください。

最後に2つの平均の差を求めます。

teffects raのATEと同じ値（−239）が算出されましたでしょうか。
回帰調整による傾向スコア分析も便利な分析手法ですが、統計モデルを利用する際に分析者が常に注意すべきこととしてモデルの“定式化”という問題があります。
誤った回帰モデルを利用しているとすれば、当然、ATEの推定値にもズレが生じます。
この問題に対応する推定方法としてAIPW（Augmented inverse-probability weighting）という推定手法が用意されていますので、最後にその利用方法をご紹介しましょう。
AIPWによるATEの推定は次の要領で行います。

1. 処置モデル（Treatment model ： デフォルトはロジットモデル）のパラメータを推定し、各標本に対応する確率の逆数を求めます。
   これはステップ3で加重として利用します。
2. 処置群とコントロール群で個別に回帰分析を実行し、潜在的なアウトカムを参加者全員に対して求めます。
3. 最後の平均を計算する時にステップ1で求めた確率の逆数を利用して加重平均を求めます。
   これを実際にStataで実行してみましょう。
   コマンドは teffects aipwになります。

2つの推定結果を比べて見てみましょう。
次ページ図6の中ほどの枠で囲った部分は各群で回帰モデルを推定した情報ですので、まったく同じものになります。
AIPWでは一番下の処置のモデルを利用した推定した確率の逆を利用して加重平均を求めますので、結果として一番上に示すATEが−232となりました。
AIPWの場合、アウトカムモデル（回帰モデル）、または、処置モデル（ロジットモデル）のどちらか一方が正しければ、単独のRAやIPWによる処置効果推定量よりも、より効率的な推定量2を提供するというメリットがあります。
傾向スコア分析（teffects psmatchコマンド）をサポートした理由と、teffects raコマンドの具体的な計算方法、そして二重に堅牢な推定手法である aipwコマンドをご紹介しました。
Stataはもう一つの二重に堅牢な推定手法であるipwra（inverse-probability-weighted regression adjustment）/生存時間をアウトカムとする場合の分析手法もサポートしています。
是非、新しい分析手法にチャレンジしてください。