30. 重回帰分析 13（傾向スコア）

• このセクションで使っている R packages

library(broom)
library(devtools)
library(lmtest)
library(margins)
library(MatchIt)
library(sandwich) 
library(tidyverse)
library(patchwork)
library(stargazer)

1. 傾向スコアの定義

• 傾向スコア \(e_i(X_i)\)
  観測された共変量 \(X_i = (X_{1i}, X_{2i},..., X_{kj})\) で条件付けた、処置される（\(D=1\)）確率

\[e_i(X_i) = Pr(D_i = 1|X_i)\] ただし \[(0≦e_i(X_i)≦1)\]

• 傾向スコアとは、\(X\) で条件付けられた処置 \(D\) が 1 になる確率

• 傾向スコアは確率　→　 0 から 1 の間の値をとる (Rosenbaum and Rubin, 1983)

• 極端に小さい値（0.1 など）や極端に大きな値（0.9 など）は使いにくいので、分析のためには望ましくない

2. なぜ傾向スコアが必要なのか？

• 同じ個体に関して「処置を行い」かつ「処置を行わない」ことはできない
  

• 実験研究では、複数の個体を用意し、それらの個体を無作為に「処置群」と「統制群」に割り付けて比較可能な集団を作り、平均処置効果 (ATE) を推定できる
  

• しかし観察研究では、処置の割り付けが無作為ではない
  

• 「処置群」と「統制群」をそのまま比較しても因果効果を適切に推定できない
  

• そこで、処置を受ける確率（= 傾向スコア）が同じ個体同士をペアで比較すればいいのではないかと考える

• 傾向スコアで条件付けしてマッチングし、観測される共変量が同じ個体をペアにして、異なる個体同士を同じ個体として取り扱う
  →　セレクションバイアスを除外して因果効果を推定できる
  

• ランダム化比較試験 (RCT) などの実験研究 では因果推論ができるのに、調査・観察データではできない理由
  →　調査・観察データでは「処置群」と「統制群」が交換できないが、ランダム化比較試験 (RCT) では「処置群」と「統制群」が交換できるから

分析方法	処置を受けることを決める人	「処置群」と「統制群」交換可能性	因果推論の可否
ランダム化比較試験 (RCT)	分析者	可能　　→　セレクションバイアスなし	因果推論ができる
調査・観察データの単純比較	被験者	不可能　→　セレクションバイアスあり	因果推論ができない

• 因果効果は「2 つの潜在的アウトカム間の差」として定義される
  
• しかし、これらの 2 つの潜在的アウトカムは、いずれか一方しか観察できない
  →　個々の個体についての因果効果を推定することはできない
  →　しかし、ランダム化比較試験では、観測しなかった半分の集団と観測したもう半分の集団は、処置を施されたかどうかという違い以外に大きな差がないと考える　　
• 例えば調査対象である有限母集団が 100人（男性50人、女性50人）だとする
  
• ここから無作為に 50 人ずつ抽出することを考える
• 50人ずつ交換可能な処置群と統制群に人々を割り振って平均処置効果を比較する
  →　仮に100人全員が処置を受け受けた場合、あるいは受けなかった場合の結果をかなりの程度推論できる

2.1 調査・観察データの場合

• 100人の有限母集団の時 　　
• 最悪の場合・・・男性50人、女性50人ずつ割り振られてしまうことがあり得る
  
• 「ランダムな割り付けは、自然には実現されない」（ローゼンバウム 2021）

→　男女比にばらつきが出る　→　性別による差が生じてしまう
→　この状況は交換可能ではない（＝グループを入れ替えられない）
→　平均因果効果 (ATE) を推定できない

• 「自然に任せていたのでは公平な比較は望むべくもない」（ローゼンバウム 2021）

2.2 ランダム化比較試験 (RCT) の場合：

• 事前にセレクションバイアスが生じないような割り振りができる
  → 例：性別でブロッキングする
  →　それぞれのブロック内でランダムに半数を処置群と統制群に割り付ける
  →　ブロックした要素（この場合は性別）に関して、処置群と統制群の違いがなくなる

• ランダム化比較試験 (RCT) では被験者が分析を受けるかどうかを決めるのは分析を実施する「研究者」
  →　「処置群」と「統制群」が「均等」である
  →　「処置群」と「統制群」が交換可能である
  →　共変量の分布が確率的に同じになる　→　交絡を除去できる
  →　セレクションバイアスがない
  →　「処置群」と「統制群」の結果変数の平均値を群間比較する
  →　因果推論ができる

• 調査・観察データでは分析を受けるかどうかを決めるのは「被験者」（= セルフセレクション）
  →　「処置群」と「統制群」が「均等」でない
  →　「処置群」と「統制群」が交換できない
  →　共変量の分布が確率的に異なる　→　交絡を除去できない
  →　セレクションバイアスがある
  →　因果推論はできない

• 調査・観察データを使って因果推論を可能にするための方法
  →　「処置群」と「統制群」が交換可能になるように調整する
  →　具体的には、重回帰分析でコントロール変数（= 制御変数）をモデルに投入する

3. 重回帰分析によるバイアスの除去

• 「回帰分析」「分散分析」「共分散分析」の違い

分析の種類	英語表記	説明変数の種類	説明変数の数
単回帰分析	simple regression analysis	連続変数	1つ
分散分析	ANOVA: analysis of variance	ダミー変数	1つ
重回帰分析 (1)	multiple-regression analysis	連続変数	2つ以上
重回帰分析 (2)	multiple-regression analysis	連続変数 + ダミー変数	2つ以上
共分散分析	ANCOVA: analysis of covariance	連続変数 + ダミー変数	2つ以上

• 「共分散分析」は「重回帰分析」のひとつと分類できる
  
• 「重回帰分析 (2)」と「共分散分析」は同じ意味
  
• 統計的因果推論を実行するモデルは処置変数（ダミー変数）を含むので、共分散分析を使う
  
• 統計的因果推論では「重回帰分析」より「共分散分析」と呼ぶのがより正確

3.1 共分散分析の利点

• 共変量 \(X\) の値で条件付ける
  →　観測された共変量 \(X\) のみがセレクションの原因なら
  （＝「強い意味での無視可能性」が成り立つなら）
  （＝「共変量 \(X\) の値が同じ個体同士では処置の割付確率が同じ」なら）
  →　処置群と統制群で共変量 \(X\) の値が同じ個体を 2 群間で比較
  →　共変量 \(X\) が特定の値をとる場合（= 共変量 \(X\) に条件づけられた場合）の処置効果 (\(ATE\)) が推定できる
  
• 重回帰分析（＝共分散分析）では ATE は推定できるが 「処置群の平均処置効果 (ATT) 」は推定できない
  
• 重回帰分析で処置変数に付いている回帰係数（= 回帰直線の傾きのこと）は、共変量 \(X\) ごとに条件付けた処置効果のこと

\[E[Y_i(1)|D_i=1, X_i = x] - E[Y_i(0)|D_i = 0, X_i = x]\\ = E[Y_i(1) |X_i = x] - E[Y_i(0)|X_i = x]\\ =E[Y_i(1) - Y_i(0)|X_i = x]\]

• その条件付けに対して、さらにもう一度期待値を取る
  →　平均処置効果 (\(ATE\)) になる

\[E[[E[Y_i(1) - Y_i(0)|X_i = x]] = E[Y_i(1) - E[Y_i(0)] = ATE\]

• \(ATE\) を計算するためにはセレクションの原因となる共変量全てをコントロールする必要がある

• 共変量 \(X\) の数が多い場合:
  

• 共変量の数が増えても、セレクションの原因となる共変量全ての値に条件付けをすれば、重回帰分析によって平均処置効果 (\(ATE\)) を求めることができる

• 観測された共変量 \(X\) のみがセレクションの原因なら
  （＝強い意味での無視可能性が成り立つなら）
  →　処置群と統制群で共変量 \(X\) の値が同じ個体を 2 群間で比較
  →　共変量 \(X\) が特定の値をとる場合（= 共変量 \(X\) に条件づけられた場合）の処置効果 (\(ATE\)) が推定できる

\[E[Y_i(1)|D_i=1, X_{1j}, X_{2j},...,X_{kj}] - E[Y_i(0)|D_i = 0, X_{1j}, X_{2j},...,X_{kj}]\\ =E[Y_i(1) - Y_i(0)|X_{1j}, X_{2j},...., X_{3j}]\]

• その条件付けに対して、さらにもう一度期待値を取る
  →　平均処置効果 (\(ATE\)) になる

\[E[[E[Y_i(1) - Y_i(0)|X_{1j}, X_{2j},...,X_{kj}]] = E[Y_i(1) - E[Y_i(0)] = ATE\]

3.2 共変量が多い時の問題点

• サンプルサイズ (\(N\)) が小さく、しかも共変量が多いと・・・
  →　比較可能なサンプルが存在しない可能性が出てくる
  →　重回帰分析によって共変量をコントロールできない
  →　重回帰分析による推定が難しくなる

・例: 共変量が 5 個のモデルを考える => \(X_1, X_2, X_3, X_4, X_5\)

• 全ての共変量の値が 0 と 1 の二値変数だとする
  
• 共変量の値のパターンは \(2^5 = 32\) 通り
  →　重回帰分析によって共変量をコントロールするためには・・・
  →　32 の共変量の組み合わせパターンに対して、処置群 (\(D = 1\)) と統制群 (\(D = 0\))を比較するためには、処置群と統制群それぞれに少なくても一つずつ固体が必要
  →　観測数は \(2 x 32 = 64\) 以上必要
  →　そうでなければ、重回帰分析によって共変量をコントロールできない

・例: 共変量が 5 個のモデルを考える => \(X_1, X_2, X_3, X_4, X_5\)

• 全ての共変量の値が 0 〜 100 の連続変数だとする
  →　重回帰分析によって共変量をコントロールするためには・・・
  →　全ての共変量に関して、処置群(\(D = 1\)) と統制群 (\(D = 0\)) を比較するためには、膨大な量の観測数が必要
  →　共変量と結果変数の関係を正しく定式化できない可能性が高い

3.2.1 N が十分大きい場合（2009年総選挙データ）

• ここで、日本の2009年の総選挙データを使って、次のようなモデルを考えてみる
• 「候補者が得た得票率 (voteshare)」を「候補者が現職か否か (inc)」で回帰してみる
  
• 知りたいのは「候補者が現職ならより得票するのかどうか」ということ
  
• ここでは処置変数である「候補者が現職か否か (inc)」以外に、次の 4 つの共変量をモデルに含めてみる
  
• ここで使っているモデルは（処置変数も含めて）共変量が 5 個
  →　共変量の値のパターンは \(2^5 = 32\) 通り
  →　重回帰分析によって共変量をコントロールするためには・・・
  →　32 の共変量の組み合わせパターンに対して、処置群 (\(D = 1\)) と統制群 (\(D = 0\))を比較するためには、処置群と統制群それぞれに少なくても一つずつ固体が必要
  
• 観測数は \(2 x 32 = 64\) 以上必要
  →　そうでなければ、重回帰分析によって共変量をコントロールできない
  
• ここでのサンプル数 (N) は 1124
  →　重回帰分析を行うには十分なサンプル数

変数名	変数の役割	説明
voteshare	結果変数	候補者が得た得票率 (%)
inc	処置変数	候補者が現職（小選挙区当選者）なら 1、そうでなければ 0
more_exp	共変数	候補者が費やす選挙費用が平均以上なら 1, 平均以下なら 0
male	共変量	候補者が男性なら 1、女性なら 0
ldp	共変量	候補者が自民党所属なら 1、それ以外なら 0
seshu_dummy	共変量	候補者が世襲なら 1、それ以外なら 0

• 衆議院議員総選挙の得票データ hr96-21.csv をダウンロード
  

• 1996年から2017年までの総選挙データ(hr96-21.csv) を読み込む

df <- read_csv("data/hr96-21.csv", na = ".")

• 変数を確認

names(df)

 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"

• 二値変数を作成してデータフレームに追加する

df <- df %>% 
  mutate(male = ifelse(gender == "male", 1, 0), # 男性なら 1、女性なら 0 の変数 male
         inc = ifelse(status == 1, 1, 0),　　 　 # 現職なら 1、それ以外なら 0 の変数 inc
         ldp = ifelse(seito == "自民", 1, 0),    # 自民なら 1、それ以外なら 0 の変数 ldp
         more_exp = ifelse(exp > 7551393, 1, 0))# 選挙費用が平均以上なら 1、平均以下なら 0 の変数 more_exp

• 必要な変数だけを抜き出す

df1 <- df %>% 
  dplyr::filter(year == 2009) %>%  # 2009年総選挙データだけを抜き出す  
  dplyr::select(voteshare, inc, more_exp, male, ldp, seshu_dummy) # 必要な変数を指定

• 重回帰分析を実行

model_1 <- lm(voteshare ~ inc + more_exp + male + ldp + seshu_dummy, data = df1,
               family = binomial(link = "logit"))

• 重回帰分析の結果を表示

stargazer(model_1, type = "html")

	Dependent variable:
	voteshare
inc	21.506***
	(1.496)
more_exp	16.836***
	(1.363)
male	4.734***
	(1.347)
ldp	-8.642***
	(1.537)
seshu_dummy	4.920***
	(1.693)
Constant	10.923***
	(1.244)
Observations	1,124
R2	0.471
Adjusted R2	0.469
Residual Std. Error	16.174 (df = 1118)
F Statistic	199.450*** (df = 5; 1118)
Note:	p<0.1; p<0.05; p<0.01

• 重回帰分析の結果を可視化して表示

theme_set(theme_bw(base_family = "HiraKakuProN-W3"))

coef_sjplot1 <- sjPlot::plot_model(model_1, 
                          show.values = T,       #推定値の係数を表示する
                          show.p = T,            #統計的有意性を*で表示する  
                          vline.color = "black", #推定値=0の縦線の色を黒に指定  
                          axis.title = "推定値", title = "N = 1124 (2009年総選挙)",
                          digits = 1) # 推定値の表示を小数点1桁に指定
coef_sjplot1

• 処置変数 (inc) を含め、全ての共変量が得票率に影響を与えており、統計的に有意だと推定

3.2.2 N が小さい場合（2009年総選挙データ）

• ここで使っているモデルは共変量が 5 個
  →　共変量の値のパターンは \(2^5 = 32\) 通り
  →　重回帰分析によって共変量をコントロールするためには・・・
  →　32 の共変量の組み合わせパターンに対して、処置群 (\(D = 1\)) と統制群 (\(D = 0\))を比較するためには、処置群と統制群それぞれに少なくても一つずつ固体が必要
  
• 観測数は \(2 x 32 = 64\) 以上必要
  →　そうでなければ、重回帰分析によって共変量をコントロールできない
  
• サンプルサイズが 64 以下の場合をシミュレーションして、重回帰分析を行ってみる
  
• model_1 で使ったサンプル数は N = 1124
• ここでは N = 1124 から無作為に 60 だけ抽出して小さいサンプルサイズで重回帰分析を行ってみよう
  
• ここでのサンプル数 (N) は 60
  →　重回帰分析を行うには不十分なサンプル数

set.seed(5) # 何度シミュレーションしても同一の結果を得るためにシードを指定  

• N = 1124 から無作為に 60 だけ抽出して d60 と名前を付ける

df60 <- sample_n(df1, 60)

• 無作為抽出したデータのサンプル数 (60) と変数の数 (6) を確認

dim(df60)

[1] 60  6

• 重回帰分析を実行

model_2 <- lm(voteshare ~ inc + more_exp + male + ldp + seshu_dummy, data = df60,
               family = binomial(link = "logit"))

• 重回帰分析の結果を表示

stargazer(model_2, type = "html")

	Dependent variable:
	voteshare
inc	24.365***
	(6.410)
more_exp	21.554***
	(4.942)
male	2.744
	(4.872)
ldp	-6.650
	(7.040)
seshu_dummy	1.309
	(7.138)
Constant	7.623*
	(4.484)
Observations	60
R2	0.626
Adjusted R2	0.591
Residual Std. Error	14.382 (df = 54)
F Statistic	18.042*** (df = 5; 54)
Note:	p<0.1; p<0.05; p<0.01

• 重回帰分析の結果を可視化して表示

theme_set(theme_bw(base_family = "HiraKakuProN-W3"))

coef_sjplot2 <- sjPlot::plot_model(model_2, 
                          show.values = T,       #推定値の係数を表示する
                          show.p = T,            #統計的有意性を*で表示する  
                          vline.color = "black", #推定値=0の縦線の色を黒に指定  
                          axis.title = "推定値", title = "N = 60 (2009年総選挙)",
                          digits = 1) # 推定値の表示を小数点1桁に指定
coef_sjplot2

• サンプル数が十分大きい model_1 と比較すると、サンプル数が 60 の model_2 の分析結果は大きく異なることがわかる

coef_sjplot1 + coef_sjplot2

• model_1 では処置変数である現職 (inc) は得票率を 21.5%ポイント上げるという結果
  

• model_2 では現職だと得票率を 24.4 %ポイント上げると過大推定
  

• model_1 ではより選挙費用を使っている候補者 (more_exp)は得票率を 16.8%ポイント上げるという結果

• model_2 ではより選挙費用を使っている候補者 (more_exp)は得票率を 21.6%ポイント上げると過大推定
  

• model_1 では男性候補者 (male) は得票率を 4.7%ポイント上げるという結果

• model_2 では性別は得票率と無関係という結果(統計的に有意ではない）
  

• model_1 では自民党候補者 (ldp) は得票率を 8.6%ポイント下げるという結果

• model_2 では自民党候補者は得票率と無関係という結果(統計的に有意ではない）
  

• model_1 では世襲候補者 (seshu_dummy) は得票率を 4.9%ポイント上げるという結果

• model_2 では世襲候補者と得票率は無関係という結果(統計的に有意ではない）

• この結果はほんの一例
  

• set.seed(5) と設定せずにシミュレーションすると、実行する度に大きく異なる結果がでる
  → 結果が安定しない
  → 共変量が多く、サンプル数が少ないと重回帰分析による正しい推定ができない

4. セレクションバイアスの問題点

• RCT の場合・・・処置を受ける確率は実験する人が決める
  →　その確率に基づいてランダムに「処置群」と「統制群」を決める
  →　セレクションバイアスが生じない
  →　因果効果を推定できる

• 調査・観察データの場合・・・処置を受けるかどうかは被験者が決める
  →　処置を受ける確率が、共変量の値に依存して決まる
  →　共変量の値によって、処置を受けやすい個体と、処置を受けにくい個体が発生してしまう
  →　処置を受けやすい個体のグループと、処置を受けにくい個体のグループが「交換不可能」
  →　セレクションバイアスが生じる
  →　因果効果を推定できない

・例：「通院（病院に行ったこと）」と「健康状態」の関係を調査・観察データを使って調べる

• データからわかる事実：
• 通院しなかった人の方が健康（健康状態の平均値： 3.21 < 3.93）
  →　常識と反する結果！
  
• その理由: 「通院するかどうか」ということは本人が選べるから
  →　セルフセレクション (self-selection)
  →　「病院に行く」のは「もともと健康状態が悪い人」
  →　もともと健康状態の良い人はそもそも「病院には行かない」

結論・ 調査・観察データを使って分析すると「病院に行くと不健康になる」という誤った推論をしてしまう可能性がある

• ここでは、通院と健康状態の事例におけるセレクションバイアスについて考える

4.1 割り振り（RCT の場合）

• RCT の場合・・・処置を受ける確率は実験する人が決める
  →　その確率に基づいてランダムに「処置群」と「統制群」を決める
  →　個体が処置を受ける確率 \(p\) は、処置群と統制群で等しい
  
• ここでいう「等しい」とは必ずしも \(p = 0.5\) とは限らない
  
• どの個体も「処置群」（あるいは「統制群」）に割り振られる確率が「等しい」という意味
  = 個体の共変量の値の違い（もともとの健康状態の値の違い）が処置群と統制群の間で存在しない
  = 処置群と統制群は交換可能
  
• どの個体も処置群に割り振られる確率は等しく \(p\) だったが、たまたまある個体は処置群に、ある個体は統制群に割り振られたということ
  
• 例えば、処置群の最上位左端の「健康」という個体は、処置群に割り振られる確率は \(p\) であり、たまたま処置群に割り振られた
  
• 例えば、統制群の最上位右側の「健康」という個体は、処置群に割り振られる確率は \(p\) だったが、たまたま統制群に割り振られた
  
• 処置を受ける（通院する）確率が、共変量（もともとの健康状態）の値に依存しない
  →　共変量（もともとの健康状態）の値によって、処置を受けやすい（通院する）個体と、処置を受けにくい（通院しない）個体が発生しない
  →　処置を受けやすい個体のグループと、処置を受けにくい個体のグループが「交換可能」
  →　セレクションバイアスが生じない
  →　因果効果を推定できる

4.2 割り振り（調査・観察データの場合）

• 調査・観察データの場合・・・処置を受けるかどうかは被験者が決める
  →　個体が処置を受ける（通院する）確率が、共変量（もともとの健康状態）の値に依存して決まる
  →　個体が処置を受ける（通院する）確率が、処置群と統制群の間で異なる
  
• 個体ごとに共変量の値（もともとの健康状態の値）が違う
  →　その共変量の値によって処置を受ける確率が決まる
  →　共変量（もともとの健康状態）の値によって、処置を受けやすい（通院する）個体と、処置を受けにくい（通院しない）個体が発生してしまう
  →　処置を受けやすい個体のグループと、処置を受けにくい個体のグループが「交換不可能」
  →　セレクションバイアスが生じる
  →　因果効果を推定できない

共変量の数が少ない場合にはどのような分析が可能なのか？

• ここでは共変量は 2 つしかない： \(X_1\) と \(X_2\)
  \(X_1　∈ \{男性、女性\}\)
  \(X_2　∈ \{健康、不健康\}\)

• 共変量によって処置を受ける（通院する）確率が変わることが分かっている
  →　共変量の値（もともとの健康状態）に条件付けて分析すれば良い
  
• 共変量の組み合わせの値ごとにグループを分ける
  →　同じ色と同じ文字同士（共変量の値が同じ個体同士）を処置群と統制群の間で比較する
  
• 例えば「健康な女性（健康）」は、平均的に処置群と統制群に等しい頻度で割り振られる必要がある（下の図では統制群には健康な人が三人おり、処置群には健康が一人しかいないが、この頻度が同じになるよう割り振る）
  → 2 つのグループの処置効果は比較できる

比較対象	内容
健康と健康	: 健康な女性同士の比較
健康と健康	: 健康な男性同士の比較
不健康と不健康	: 不健康な女性同士の比較
不健康と不健康	: 不健康な男性同士の比較

• 共変量の値に条件付けたそれぞれの平均処置効果（ここでは 4 つ）を推定
  →　全体の期待値（＝平均値）を計算する
  →　全体の処置効果を推定できる　　

共変量の数が多い場合にはどのような分析が可能なのか？

• 共変量の数が数10個、数100個の場合、どのように考えればいいのか？
  →　組み合わせのパターンが非常に増える
  →　各ブロックごとに比較しようとしても、比較する対象が存在しないかもしれない
  （下の図の統制群の右下のセルがそれに該当する）
  →　共変量が多く、観察数が少ないと比較する対象が存在しない確率が増える
  →　比較対象がなければ、そのブロックに関する処置効果の推定ができない
  
• ここでの問題 = 「処置群と統制群で「処置を受ける確率が違うこと」
  →　共変量の値をセルごとに全て揃えなくても、処置を受ける確率が同じもの同士を比較すればいいのではないか 　　

• 事例①：処置群の右上の枠にいる健康（健康な男性）は、もともとの健康状態の値が健康なのだから、確率的に考えると本来なら統制群（通院しない）に多くいるタイプ
  →　しかし、たまたま処置群に割り振られた
  →　この健康（健康な男性）は統制群にいる他の個体と似ているはず
  →　この健康（健康な男性）は統制群の比較対象として（統制群の代表として）上手く使えるのではないか、と考える
  

• 事例②：統制群の左下の枠にいる不健康（不健康な女性）は、もともとの健康状態の値が悪いのだから、確率的に考えると本来なら処置群（通院する）に多くいるタイプ
  →　たまたま統制群に割り振られた
  →　この不健康（不健康な女性）は処置群にいる他の個体と似ているはず
  →　この不健康（不健康な女性）は処置群の比較対象として（処置群の代表として）上手く使えるのではないか、と考える

• 変量の値をセルごとに全て揃えなくても、処置を受ける確率が同じもの同士を比較すればいいのではないかという発想を利用するのが傾向スコア

5. 傾向スコア

• 「傾向スコア」に関しては、高橋将宣著『統計的因果推論の理論と実装』(2022) の Chapter 10 の内容が最新で解説が大変わかりやすい
  
• ここでは、基本的に同書 Chapter 10 で使われているデータを使い、同書に沿って「傾向スコア」について解説する
  
• 傾向スコアをより良く理解するために処置群と統制群の「共変量のバランス」に着目し、人工的に生成したデータを使っての解説を試みる
  
• 「傾向スコア」をマスターしたい人は同書を熟読することを強くおすすめします

5.1 傾向スコアの定義

• 傾向スコア (propensity score) とは
  ・バランシングスコア (blancing score) と呼ばれる関数の一種
  ・バランシングスコアとは、観測される共変量 \(X_i\) の関数 \(b(X_i)\) が与えられた時の \(X_i\) の条件付き分布が、処置群と統制群において同じになる関数のこと
  
• 「観測される共変量 \(X\) で条件付けた処置 \(D\) が 1 になる）確率 」
  = 予測された共変量 \(X_i = (X_{1i}, X_{2i},..., X_{kj})\) で条件付けた、処置される（\(D=1\)）確率
  
• 処置群に割り当てられやすい個体と統制群に割り当てられやすい個体を区別してスコア化したもの

傾向スコア \(e_i(X_i)\) は次の式で表すことができる

\[e_i(X_i) = Pr(D_i = 1|X_i)\] ただし \[(0≦e_i(X_i)≦1)\]

• 傾向スコアは確率　→　 0 から 1 の間の値をとる

• 極端に小さい値（0.1 など）や極端に大きな値（0.9 など）は使いにくいので、分析のためには望ましくない

• 処置を受けやすい個体同士、あるいは処置を受けにくい個体同士を比較する
  →傾向スコアが同じ(＝近い)もの同士を処置群と統制群として比較
  →平均処置効果 (ATE) が推定できる

5.2 共変量のバランス

• ここで使うのは次の 5 つの変数

y0t	潜在的結果変数: \(Y_i(0)\)	観察不可能（理論上の変数）
y1t	潜在的結果変数:\(Y_i(1)\)	観察不可能（理論上の変数）
t1	「t1 = 1: 処置あり」, 「t1 = 0: 処置なし」	観察可能
y	結果変数	観察可能
x1	共変量	観察可能

• ここでは、高橋将宣著『統計的因果推論の理論と実装』(2022) のサポートページで公開されている「本書で使用するデータ（zip）」をダウンロードし data10a.csv を使う

• データを読み取る

df10a <- read_csv("data/takahashi/data10a.csv") # パスは皆さんのパソコン環境に応じて変更して下さい

• 変数を確認する

DT::datatable(df10a)

• 通常のデータ分析では2 つの潜在的結果変数 (\(Y_i(0)\) と (\(Y_i(1)\)) は観察できない
  
• なぜなら、例えば、コロナウィルスのワクチン「接種を受けない = \(Y_i(0)\)」と同時に「接種を受ける = \(Y_i(1)\)」ことはできないから
  
• しかし、ここではドラえもんにお願いして y0t と y1t という2 つの観察不可能な変数が観察できると想定
  →　「因果効果の真値 (ATE)」 と「処置群の平均処置効果 (ATT) の真値」を計算できる

平均処置効果 (ATE) の真値

ATE <- mean(df10a$y1t) - mean(df10a$y0t)
ATE

[1] 9.95

処置群の平均処置効果 (ATT) の真値

ATT <- mean(df10a$y1t[df10a$t1 == 1]) - mean(df10a$y0t[df10a$t1 == 1])

ATT

[1] 9.090909

人工的に設定したATT と ATE の真値 (data.10a.csv) ・ATE の真値 = 9.95
・ATT の真値 = 9.09

• 5 つの変数のうち、私たちが観察可能なのは 3 つだけ: t1, y, x1
• 私たちが観察できるこれらのの変数を使って平均処置効果をナイーブに推定してみる

ナイーブな推定値

tre1 <- mean(df10a$y[df10a$t1 == 1]) # 観察できる y3 （処置を受けた人 t = 1）の平均値  
tre1

[1] 83.72727

con1 <- mean(df10a$y[df10a$t1 == 0]) # 観察できる y3 （処置を受けない人 t = 0）の平均値  
con1

[1] 66.33333

naive <-  tre1 - con1
naive

[1] 17.39394

• 17.4 という推定値が得られたが、これは平均処置効果の真値 (ATE = 9.95) をかなり過大評価している
  
• 観測値の散布図を描いて、このことを可視化してみる

df10a %>% 
  ggplot(aes(x1, y)) + 
  geom_point(aes(color = as.factor(t1))) +
  geom_hline(yintercept = con1, color = "red", linetype = 2) +
  geom_hline(yintercept = tre1, color = "blue", linetype = 2) +
annotate("label", 
           label = "ナイーブな推定値\n= 17.4", 
           x = 69, y = 75,
           size = 3, 
           colour = "purple", 
           family = "HiraginoSans-W3")

処置群 (t = 1) と統制群 (t = 0) での x1 の分布　　

df10a$t1 <- as.factor(df10a$t1)

df10a |> ggplot(aes(x = t1, y = x1)) +
  geom_jitter(aes(color = t1), 
              position = position_jitterdodge(jitter.width = 0.2,
                                              jitter.height = 0),
              show.legend = FALSE) +
  geom_boxplot(aes(fill = t1),
               alpha = 0.5) +
  labs(x = "0:統制群、1:処置群", y = "y", fill = "",
       caption = "") +
  theme_bw(base_family = "HiraKakuProN-W3")

summary(df10a$x1[df10a$t1==0])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  66.00   70.00   74.00   74.67   78.00   88.00 

summary(df10a$x1[df10a$t1==1])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   73.0    77.5    83.0    83.0    89.0    92.0 

• 処置群 (t = 1) と統制群 (t = 0) における x1 の分布はかなり異なることがわかる

5.3 バランシングスコアを使った調整

• ここでバランシングスコア b(X) を使って、共変量 x1 の分布が処置群と統制群で同じになるように調整してみる
  
• このときのバランシングスコア b(X) は共変量 x1
  → 共変量 x1 の値が小さい順に並べ替えてバランスをとる

• 被験者 2 被験者 3 は共分散の値 x1 = 70 が同じ
  
• しかし、処置 t1 がどちらも 0 なので群間比較できない
• 被験者 4 と被験者 5 は共分散の値 x1 = 73 が同じ
  被験者 4 は t1 = 1（処置群）、被験者 5 は t1 = 0（統制群）なので群間比較できる
  
• 被験者 4 と被験者 5 の組を事実上「同一人物」と見なして分析する
  
• 被験者 11 と被験者 12 は共分散の値 x1 = 78 が同じ
  被験者 11 は t1 = 1（処置群）、被験者 12 は t1 = 0（統制群）なので群間比較できる
  
• 被験者 11 と被験者 12 の組を事実上「同一人物」と見なして分析する

処置群 (t1 = 1) における共変量 x1 の平均値

• 処置群 (t1 = 1) における比較可能な共変量 x1 は 73 と 78 なのでその平均値を求める

(73 + 78)/2

[1] 75.5

統制群 (t1 = 0) における共変量 x1 の平均値

• 処置群 (t1 = 0) における比較可能な共変量 x1 は 73 と 78 なのでその平均値を求める

(73 + 78)/2

[1] 75.5

• どちらの平均値も 75.5 で共変量のバランスがとれている
  →　処置群と統制群の間で共変量 x1 の分布が同じという意味

平均処置効果 (ATE) の推定値を計算してみる

処置群 (t1 = 1) における結果 y の平均値

(77 + 81)/2

[1] 79

統制群 (t1 = 0) における結果 y の平均値

(73 + 65)/2

[1] 69

ATE は両者の差だから

79-69

[1] 10

人工的に設定したATT と ATE の真値 (data10a.csv) ・ATE の真値 = 9.95
・ATT の真値 = 9.1
・ナイーブな推定値 = 17.4
・ATE の推定値 = 10

• ナイーブな推定値 (17.4) と比較するとかなり ATE の真値 (9.95) や ATT の真値 (9.1) にかなり近い値 (10) を推定できている

5.4 傾向スコアの使い方

• ここで使うのは次の 5 つの変数

y0t	潜在的結果変数: \(Y_i(0)\)	観察不可能（理論上の変数）
y1t	潜在的結果変数:\(Y_i(1)\)	観察不可能（理論上の変数）
t1	「t1 = 1: 処置あり」, 「t1 = 0: 処置なし」	観察可能
y	結果変数	観察可能
x1	共変量	観察可能
x2	共変量	観察可能
x3	共変量	観察可能
ps1	傾向スコア	観察可能

• x1 に x2 と x3 という 2 つの新たな共変量が追加され、共変量が 3 つになった
  

• 共変量が x1 だけであれば、バランシングが達成できた
  

• 共変量が 3 つになると単純なバランシングが難しい
  

• どうするか？
  

• 解決策　→　傾向スコア ps1 を使う
  

• 次に、傾向スコアを作成する

• ここでは、高橋将宣著『統計的因果推論の理論と実装』(2022) のサポートページで公開されている「本書で使用するデータ（zip）」をダウンロードし data10b.csv を使う

• data10b.csv を読み取る

df_10b <- read_csv("data/takahashi/data10b.csv")

DT::datatable(df_10b)

• 処置群 (t = 1) と統制群 (t = 0) における \(x1, x2, x3\) の分布を可視化してみる

• 処置群と統制群のばらつきに系統的な差がある
  →　このままの状態では結果変数 y を処置群と統制群で比較できない

解決策 →　傾向スコアを使う 　　

• 傾向スコア \(e_i(X_i)\) とは
  ・バランシングスコア (blancing score) と呼ばれる関数の一種
  ・バランシングスコアとは、観測される共変量 \(X_i\) の関数 \(b(X_i)\) が与えられた時の \(X_i\) の条件付き分布が、処置群と統制群において同じになる関数のこと
  
• 「観測される共変量 \(X\) で条件付けた処置 \(D\) が 1 になる）確率 」
  = 予測された共変量 \(X_i = (X_{1i}, X_{2i},..., X_{kj})\) で条件付けた、処置される（\(D=1\)）確率
  
• 処置群に割り当てられやすい個体と統制群に割り当てられやすい個体を区別してスコア化したもの

\[e_i(X_i) = Pr(D_i = 1|X_i)\] ただし \[(0≦e_i(X_i)≦1)\]

• 傾向スコアは確率　→　 0 から 1 の間の値をとる

• 処置を受けやすい個体同士、あるいは処置を受けにくい個体同士を比較する
  → 傾向スコアが同じ(＝近い)もの同士を処置群と統制群として比較
  →平均処置効果 (ATE) が推定できる
  

• ここでは傾向スコアは ps1（計算方法は後で解説）
  

• ps1 の値が小さい順に並べ替えている

傾向スコア値が近い被験者をペアにして比較する

• 被験者 13 被験者 11 は傾向スコア値 ps1 = 0.26 が同じ
  被験者 13 は t1 = 1（処置群）、被験者 11 は t1 = 0（統制群）なので群間比較できる
  
• 被験者 7 と被験者 5 は傾向スコア値が近い：
  
• 被験者 7 の ps1 = 0.26、被験者 5 の ps1 = 0.35
  
• しかし、処置 t1 がどちらも 0 なので群間比較できない
  
• ここで群間比較可能な 4 つのペア： 13 & 11, 6 & 4, 10 & 9, 15 & 17

処置群 (t1 = 1) における結果 y の平均値

(73 + 77 + 77 + 85)/4

[1] 78

統制群 (t1 = 0) における結果 y の平均値

(65 + 66 + 72 + 75)/4

[1] 69.5

ATE は両群の平均値の差だから

78-69.5

[1] 8.5

人工的に設定したATT と ATE の真値 (data10b.csv) ・ATE の真値 = 9.95
・ATT の真値 = 9.1
・ナイーブな推定値 = 17.4
・ATE の推定値 = 8.5

• ナイーブな推定値 (17.4) と比較すると ATE の真値 (9.95) や ATT の真値 (9.1) にかなり近い値 (8.5) を推定できている
  
• 共変量の分布を確認する
  
• 処置群 (t = 1) と統制群 (t = 0) における x1, x2, x3 の分布を可視化してみる

• 傾向スコアを使ったおかげで、処置群 (t = 1) と統制群 (t = 0) における x1, x2, x3 の分布はほぼ同じになった
  →　バランシングしている

5.5 傾向スコア定理

• Rosenbaum & Rubin (1983) は傾向スコアに関して次の2 つが成り立つことを示した

定理1 （バランシング）

• 処置の割り付け \(T\) と観測された共変量 \(X\) は、傾向スコア \(e(X)\) が与えられた時、条件付き独立である
  →　傾向スコア \(e(X)\) が同じであれば、処置群と統制群の間で共変量 \(X\) の分布は同じである

\[X ⊥ T|e(X)\]

定理2 （条件付き独立）

• 傾向スコア \(e(X)\) が与えられれば、潜在的結果変数 \({Y(1), Y(0)}\) と割り付け変数 \(T\) は条件付き独立である
  →　傾向スコア \(e(X)\) が同じであれば、処置への割り付けは無作為だとみなせる

\[{Y(1), Y(0)} ⊥ T|e(X) \]

定理1と定理2が成り立つために必要な 2 つの前提条件

「強い意味での無視可能な割り付けのための条件」

1. 条件付き独立性（無交絡性）

\[{Y(1), Y(0)} ⊥ T|e(X)\]

• 処置の割り付けに影響を与えるのは観測された共変量だけということ

傾向スコアの利点 ・観察研究では、この条件を満たすためにできるだけ多くの変数をモデルに取り込む必要がある
・その際、モデルに含めるべき「交絡因子」や含めてはいけない「中間変数」に注意する必要あり
・傾向スコア（単変量）の値が同じであれば、\(X\)（多変量）の値自体は異なっていてもかまわない

2. 条件付き正値性

\[0 < Pr(T_i=1|X) < 1\]

• \(Pr(T_i=1|X)\) は「傾向スコア」そのもの
  
• 傾向スコアは 0 と 1 の間の値をとる
  →　傾向スコアは 0 と 1 の値をとらない
  →　処置群と統制群の両方に観測値が必要
  →　処置群と統制群に「共通部分」が必要

傾向スコアのモデリングと実験研究の違い ・傾向スコアのモデリングでは処置の割り付けに影響を与えるのは観測された共変量 \(X\) のみと規定
・観測されない共変量 \(u\) のバランシングは保証しない
・実験研究では観測されない共変量 \(u\) のバランシングも保証する

• 「1. 条件付き独立性（無交絡性）」と「2. 条件付き正値性」
  ・両方を満たす場合　→　傾向スコアを使った因果推論
  ・満たさない場合　→　RDD や IV を使った因果推論

5.6 傾向スコアのモデル化

傾向スコアは、ロジステック回帰モデル (logistic regression model) でモデル化できる

• ここでは架空データを使って考える
  
• 20人の生徒がおり、前期試験を受けるとする
  
• 前期試験の点数 x1 の結果を見て、生徒自身が補講を受けるかどうか決める
  
• ここでは「前期試験の点数が低い生徒ほど補講を受ける」と想定
  
• 処置の割合を次のように設定する
  ・70点 → 4/5 = 0.8
  ・75点 → 3/5 = 0.6
  ・85点 → 2/5 = 0.4
  ・90点 → 1/5 = 0.2
  ・つまり、前期試験点数が悪い生徒ほど補講を受ける確率が高い
  →　この確率が「傾向スコア」の真値（実際の観察研究ではこの割合が不明）
  
• 前期試験の点数次第で処置（補講を受けるかどうか）が決まる
  
• 生徒 1 から生徒 5 は前期試験結果 x1 は同じ 70 点
  
• この 5 人が割り付けられる確率（＝補講を自主的に受ける確率）は 4/5 = 0.8
  →　無作為に 5 人中 4 人が補講授業を受けている
  
• 生徒 6 から生徒 10 は前期試験結果 x1 は同じ 75 点
  
• この 5 人が割り付けられる確率（＝補講を自主的に受ける確率）は 3/5 = 0.6
  →　無作為に 5 人中 3 人が補講授業を受けている

• 生徒 11 から生徒 15 は前期試験結果 x1 は同じ 85 点
  
• この 5 人が割り付けられる確率（＝補講を自主的に受ける確率）は 2/5 = 0.4
  →　無作為に 5 人中 2 人が補講授業を受けている
  
• 生徒 16 から生徒 20 は前期試験結果 x1 は同じ 90 点
  
• この 5 人が割り付けられる確率（＝補講を自主的に受ける確率）は 1/5 = 0.2
  →　無作為に 5 人中 1 人が補講授業を受けている

無視可能な割付け ・共変量の値 \(x1\) が同じ個体同士では、処置の割付 \(t1\) 確率が同じであることを「無視可能な割付け」という
・共変量の値 \(x1\) に条件付けたとき、局所的に無作為割付が実現されている
→　前期試験点数 \(x1\) の点数が同じ個体なら、その点数内で割付けられる確率（＝傾向スコア）は等しい
→　処置群と統制群で共変量 \(x1\) の値が同じ個体を 2 群間で比較
→　共変量 \(x1\) が特定の値をとる場合（= 共変量 \(x1\) に条件づけられた場合）の処置効果 (ATE) が推定できる
例）\(x1\) = 90 の時の処置効果 (ATE) = 91 - (81 + 82 + 82 + 82)/4 = 9.25

• 実際の観察研究では傾向スコアの真値がわからない
  →　ロジステック回帰モデルを使ってデータからモデル化して傾向スコアを推定する
  
• 推定に使うデータは次の通り

x1	割当変数（前期試験の点数）	観察可能
y	結果変数（期末試験の点数）	観察可能
t1	処置変数（補講に参加 t1 = 1、参加せず t1 = 0）	観察可能

• ここでは、高橋将宣著『統計的因果推論の理論と実装』(2022) のサポートページで公開されている「本書で使用するデータ（zip）」をダウンロードし data03.csv を使う

• データを読み取る

df3 <- read_csv("data/takahashi/data03.csv")
df3 <- df3 |> 
  select(x1, y3, t1) # 分析に使う変数だけに絞る  

DT::datatable(df3)

• 変数の記述統計を示す

stargazer(as.data.frame(df3), type = "html")

Statistic	N	Mean	St. Dev.	Min	Max
x1	20	80.000	8.111	70	90
y3	20	77.250	7.025	63	91
t1	20	0.500	0.513	0	1

• ロジステック回帰モデルを使って、共変量 x1 を条件とした時の処置 t1 の割り付け確率を求める

psmodel <- glm(t1 ~ x1, 
               family = binomial(link = "logit"),
               data = df3)

• 共変量 x1 を条件とした時の処置 t1 の割り付け確率の推定値を可視化してみる

pred_ps <- data_frame(x1 = seq(70, max(df3$x1), length.out = 20))
pred_ps$fit <- predict(psmodel, type = "response", newdata = pred_ps)

ggplot(df3, aes(x = x1)) +
  geom_hline(yintercept = c(0, 1), color = "gray") +
  geom_jitter(aes(y = t1), pch = 16, size = 1, width = 0.05, height = 0.05) +
  geom_line(data = pred_ps, aes(y = fit), color = "red") +
  labs(x = "x1", y = "t1", title = "共変量 x1 を条件とした時の処置 t1 の割り付け確率")

• 前期試験の点数 x1 が高い人ほど処置 t1 を受けない（＝補講を受けない）ことがわかる
  
• 予測された確率を取り出し ps3 と名前を付ける

ps3 <- round(psmodel$fitted.values, 3) # 表示桁数を指定

• 傾向スコアの真値を示す変数 ps_true を作る

ps_true <- c(rep(0.8, 5), rep(0.6, 5), rep(0.4, 5), rep(0.2, 5))

• データフレーム df4 に ps3 と ps_true を加えて中身を表示させる

df4 <- df3 |> 
  mutate(ps3 = ps3, 
         ps_true = ps_true)
df4

# A tibble: 20 × 5
      x1    y3    t1   ps3 ps_true
   <dbl> <dbl> <dbl> <dbl>   <dbl>
 1    70    74     1 0.775     0.8
 2    70    63     0 0.775     0.8
 3    70    73     1 0.775     0.8
 4    70    71     1 0.775     0.8
 5    70    74     1 0.775     0.8
 6    75    67     0 0.65      0.6
 7    75    77     1 0.65      0.6
 8    75    68     0 0.65      0.6
 9    75    77     1 0.65      0.6
10    75    78     1 0.65      0.6
11    85    88     1 0.35      0.4
12    85    77     0 0.35      0.4
13    85    76     0 0.35      0.4
14    85    86     1 0.35      0.4
15    85    78     0 0.35      0.4
16    90    81     0 0.225     0.2
17    90    91     1 0.225     0.2
18    90    82     0 0.225     0.2
19    90    82     0 0.225     0.2
20    90    82     0 0.225     0.2

• ps3 ・・・モデルから推定した傾向スコア
  
• ps_true ・・・傾向スコアの真値
  
• これでロジステック回帰モデルを使ってデータからモデル化して傾向スコアを推定できた

5.7 傾向スコアの推定 (MatchIt パッケージ)

library(MatchIt)

• 共変量 x1 を条件とした時の処置 t1 の割り付け確率を求める

matchit <- matchit(t1 ~ x1, data = df3)

• 傾向スコアの推定値 (fitted.values) を抜き出し ps_mi と名前を付ける

ps_mi <- matchit$model$fitted.values

• 傾向スコアの真値 ps_true と 2 つの推定値 ps3, ps_mi を一つのデータフレームにまとめて表示させる

df5 <- df4 |> 
  mutate(ps_mi = ps_mi) 

df5

# A tibble: 20 × 6
      x1    y3    t1   ps3 ps_true ps_mi
   <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
 1    70    74     1 0.775     0.8 0.775
 2    70    63     0 0.775     0.8 0.775
 3    70    73     1 0.775     0.8 0.775
 4    70    71     1 0.775     0.8 0.775
 5    70    74     1 0.775     0.8 0.775
 6    75    67     0 0.65      0.6 0.650
 7    75    77     1 0.65      0.6 0.650
 8    75    68     0 0.65      0.6 0.650
 9    75    77     1 0.65      0.6 0.650
10    75    78     1 0.65      0.6 0.650
11    85    88     1 0.35      0.4 0.350
12    85    77     0 0.35      0.4 0.350
13    85    76     0 0.35      0.4 0.350
14    85    86     1 0.35      0.4 0.350
15    85    78     0 0.35      0.4 0.350
16    90    81     0 0.225     0.2 0.225
17    90    91     1 0.225     0.2 0.225
18    90    82     0 0.225     0.2 0.225
19    90    82     0 0.225     0.2 0.225
20    90    82     0 0.225     0.2 0.225

• \(ps_3\) と ps_mi が傾向スコアの真値 ps_true を正しく推定していることがわかる
  
• MatchIt パッケージを使って推定した傾向スコア ps_mi と 手計算による傾向スコアの推定値 ps3 が全く同一だとわかる

6. 傾向スコアマッチング

• 「傾向スコアマッチング」に関しては、高橋将宣著『統計的因果推論の理論と実装』(2022) の Chapter 11 の内容が最新で解説が大変わかりやすい
  
• ここでは、基本的に同書 Chapter 11 で使われているデータを使い、同書に沿って「傾向スコア」について解説する
  
• 「傾向スコアマッチング」をマスターしたい人は同書を熟読することを強くおすすめします

6.1 統計的因果推論におけるマッチング

• 同じ個体に関して「処置を行い」かつ「処置を行わない」ことはできない
  

• 実験研究では、複数の個体を用意し、それらの個体を無作為に「処置群」と「統制群」に割り付けて比較可能な集団を作り、平均処置効果 (ATE) を推定できる
  

• しかし観察研究では、処置の割り付けが無作為ではない
  

• 「処置群」と「統制群」をそのまま比較しても因果効果を適切に推定できない
  

• そこで、処置を受ける確率（= 傾向スコア）が同じ個体同士をペアで比較すればいいのではないかと考える

• 傾向スコアで条件付けしてマッチングし、観測される共変量が同じ個体をペアにして、異なる個体同士を同じ個体として取り扱う
  →　セレクションバイアスを除外して因果効果を推定できる

• 下の表は 2021年総選挙に出馬した候補者のうち10人の名簿である
  

• 処置変数は候補者が世襲かどうか（0 なら非世襲、1 なら世襲）
  

• 10人のうち 3 人が世襲、7 人が非世襲候補者
  

• 結果は候補者がそれぞれの小選挙区で獲得した得票率 (%)
  

• 「年齢」と「性別」が共変量 →　共変量 \(X\) は 2 次元
  

• 2 つのペアが可能

• 40 歳で男性の「小泉進次郎」と「小倉将信」
  

• 64 歳で男性の「岸田文雄」と「石破茂」

• 年齢と性別でマッチングした結果は次のとおり

氏名	処置	結果	性別	年齢
小泉進次郎	1	79.17	male	40
小倉将信	0	51.25	male	40

• この 2 人をマッチングした理由は、この 2 人の「結果（=得票率）」の違い (79.17-51.25 = 27.92) は、年齢と性別が原因ではないと思われるから
  
• ここで未観測の交絡因子がないという仮定が正しければ、この 2 人の得票率の違い (27.92) は、「処置（世襲かどうか）」の違いが原因だといえる

氏名	処置	結果	性別	年齢
岸田文雄	1	80.67	male	64
石破茂	0	84.07	male	64

• この 2 人をマッチングした理由は、この 2 人の「結果（=得票率）」の違い (84.07 - 80.67 = 3.4) は、年齢と性別が原因ではないと思われるから
  
• ここで年齢と性別以外の未観測の交絡因子がないという仮定が正しければ、この 2 人の得票率の違い (3.4) は、「処置（世襲かどうか）」の違いが原因だといえる

因果推論の妥当性

• ここでの因果推論が妥当かどうかは、年齢と性別以外の未観測の交絡因子がないという仮定が正しければという条件に依存する
  
• つまり「条件付き独立性（無交絡性）」がみたされること

\[{Y(1), Y(0)} ⊥ T|e(X)\]

• 「処置の割り付けに影響を与えるのは、観測された共変量だけ」ということ
  
• そのためには共変量 \(X\) は多変量である必要がある
  
• 得票率に影響を与える多変量として「年齢」と「性別」だけでは不十分
  
• 交絡因子として「候補者の経験」「所属政党」「学歴」「費やした選挙費用」「美顔度」など様々な要因が考えられる
  
• 上の 2 つのマッチングにおける得票率の差が大きく異なるのは (27.9 vs 3.4)、このような多くの交絡因子が存在していることの証左
  
• 多変量 \(X\) を条件としたときに、処置に割り当てられる確率として単変量の傾向スコアに情報を集約して因果効果を推定できる

推定対象

• 統計的因果推論における推定対象は次の 2 つ：

平均処置効果	\(τ_{ATE} = E[Y_i(1)-Y_i(0)] = E[Y_i(1)] - E[Y_i(0)]\)
処置群の平均処置効果	\(τ_{ATT} = E[Y_i(1)-Y_i(0)|T_i=1] = E[Y_i(1)|T_i=1] - E[Y_i(0)|T_i=1]\)

• 傾向スコアマッチングにおける推定対象・・・処置群の平均処置効果 (ATT)

その理由：

• 処置群における個体にマッチする候補者を対照群から選ぶ
  →　マッチング後のデータは処置群の個体を中心として構成されるため
  →　処置群の平均処置効果 (ATT) を推定する

6.2 人工的に生成したデータと因果効果の真値

• ここでは人工的に生成した次の 10 の変数を使う

y0t	潜在的結果変数: \(Y_i(0)\)	通常、観察不可能だが、ここでは観察可と仮定する
y1t	潜在的結果変数:\(Y_i(1)\)	通常、観察不可能だが、ここでは観察可と仮定する
\(y3\)	結果変数	観察可能
t1	「t1 = 1: 処置あり」, 「t1 = 0: 処置なし」	観察可能
x1	共変量	観察可能
x2	共変量	観察可能
x3	共変量	観察可能
x4	共変量	観察可能
x5	共変量	観察可能
x6	共変量	観察可能

• x1,....,x6 は平均1、標準偏差1の多変量対数正規分布に従う乱数

• 誤差項 \(ε_i\) は標準正規乱数である

• ここでは、高橋将宣著『統計的因果推論の理論と実装』(2022) のサポートページで公開されている「本書で使用するデータ（zip）」をダウンロードし data11.csv を使う

• 生成したデータを読み取る

df11 <- read_csv("data/takahashi/data11.csv")

DT::datatable(df11)

• \(Y_i(0)\) に関しては次のように設定してデータを生成してある

\[Y_i(0) = 1 + 1X_{1i} - 3X_{2i} + 5X_{3i} + 7X_{4i} + 9X_{5i} + 11X_{6i} + ε_i\]

• 試しに y0t を応答変数として重回帰分析をしてみる　　　

m0 <- lm(y0t ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11)

tidy(m0)

# A tibble: 7 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    0.993    0.0610      16.3 6.53e- 53
2 x1             0.982    0.0376      26.1 8.27e-115
3 x2            -2.99     0.0396     -75.6 0        
4 x3             5.04     0.0396     127.  0        
5 x4             6.96     0.0418     167.  0        
6 x5             8.99     0.0400     225.  0        
7 x6            11.0      0.0380     291.  0        

• \(Y_i(1)\) に関しては次のように設定してデータを生成してある
  \[Y_i(1) = 5 + 2X_{1i} + 4X_{2i} - 6X_{3i} + 8X_{4i} + 10X_{5i} + 12X_{6i} + ε_i\]

• 試しに y1t を応答変数として重回帰分析をしてみる　　　

m1 <- lm(y1t ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11)

tidy(m1)

# A tibble: 7 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)     4.99    0.0610      81.8 0        
2 x1              1.98    0.0376      52.7 6.11e-290
3 x2              4.01    0.0396     101.  0        
4 x3             -5.96    0.0396    -150.  0        
5 x4              7.96    0.0418     191.  0        
6 x5              9.99    0.0400     250.  0        
7 x6             12.0     0.0380     317.  0        

処置 t1 の割り付けルール

• \(med(Y_i(0))\) は \(Y_i(0)\) の中央値
  
• \(u_{1i}\) と \(u_{2i}\) は区間 \([0, 1]\) の一様乱数

\[T_i = 1 \\if Y_i(0) > med(Y_i(0)) \& u_{1i} ≦ 0.50 \\or Y_i(0) ≦ med(Y_i(0)) \& u_{2i} > 0.75\]

• \(Y_i(0)\) が中央値より大きいとき 50% の確率で \(T_i = 1\)
  
• \(Y_i(0)\) が中央値以下のとき 25% の確率で \(T_i = 1\)

\[T_i = 0 \\if Y_i(0) > med(Y_i(0)) \& u_{1i} > 0.50 \\or Y_i(0) ≦ med(Y_i(0)) \& u_{2i} ≦ 0.75\]

• \(Y_i(0)\) が中央値より大きいとき 50% の確率で \(T_i = 0\)
  
• \(Y_i(0)\) が中央値以下のとき 75% の確率で \(T_i = 0\)

処置群の平均処置効果 (ATT) の真値

mean(df11$y1t[df11$t1==1]) - mean(df11$y0t[df11$t1==1])

[1] 2.888651

平均処置効果 (ATE) の真値

mean(df11$y1t) - mean(df11$y0t)

[1] 3.755947

人工的に設定したATT と ATE の真値 (data11.csv) ・ATT の真値 = 2.889
・ATE の真値 = 3.756

6.3 観察可能なデータを使った場合の比較

ナイーブな比較

obs_simple <- lm(y3 ~ t1, data = df11)
tidy(obs_simple)

# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)     26.3     0.906      29.0 1.84e-134
2 t1              15.8     1.45       10.9 4.02e- 26

共分散分析

obs_multi <- lm(y3 ~ t1 + x1 + x2 + x3 + x4 + x5 + x6, data = df11)
tidy(obs_multi)

# A tibble: 8 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    1.05      0.309      3.41 6.78e-  4
2 t1             3.47      0.327     10.6  4.24e- 25
3 x1             1.39      0.184      7.57 8.77e- 14
4 x2            -0.292     0.194     -1.50 1.33e-  1
5 x3             0.454     0.194      2.34 1.97e-  2
6 x4             7.47      0.205     36.5  8.05e-186
7 x5             9.66      0.197     49.1  8.91e-268
8 x6            11.2       0.187     59.5  0        

分析結果まとめ

推定方法	推定値
処置群の平均処置効果 (ATT) の真値	2.89
平均処置効果 (ATE) の真値	3.76
ナイーブな比較 (t1=1 vs t1=0)	15.8
共分散分析 (t1 の係数)	3.47

• 共分散分析では ATT を分析できない
  →　共分散分析の推定値 = 3.47（ATT の真値 2.89 とかけ離れた推定）
  
• 共分散分析が分析できるのは ATE
  
• しかし、共分散分析の推定値 = 3.47（ATE の真値 3.76 ともかけ離れた推定）
  →　その理由：モデルに全ての交差項を含めず、正しくモデリングされていないから

6.4 マッチング（復元・非復元） 　　

• 傾向スコアマッチングによる処置群の平均処置効果 (ATT) を推定するために必要な「復元と非復元」に関する考え方を R の sample() 関数を使って解説する

sample() 関数を使った復元と非復元 ・R 上で人工的に 5 人（小泉、岸田、安倍、森下、池内）の変数 list を作る

list <- c("小泉", "岸田", "安倍", "森下", "池内")
list

[1] "小泉" "岸田" "安倍" "森下" "池内"

• sample()関数を使って、5 人の中から無作為に 4 人を選ぶとする
  
• 選び方は 2 つある

1. 復元（選んだ人を list に戻して、2 人目以降を選ぶ方法）: replace = TRUE
   
2. 非復元（選んだ人を list に戻さずに、2 人目以降を選ぶ方法）: replace = FALSE

復元（replace = TRUEと指定する）

• 復元の方法で 4 人選んでみる

sample(list, 4, replace = TRUE)

[1] "池内" "森下" "安倍" "小泉"

sample(list, 4, replace = TRUE)

[1] "小泉" "森下" "池内" "小泉"

• 同一人物が選ばれることがある

非復元（replace = FALSEがデフォルト）

• 非復元の方法で 4 人選んでみる

sample(list, 4)

[1] "池内" "森下" "安倍" "小泉"

• もう一度、選んでみる

sample(list, 4)

[1] "森下" "小泉" "安倍" "池内"

• 何度やっても、同一人物が選ばれることはない

• 個体をペアにする場合、マッチング相手の選び方を「復元」もしくは「非復元」から選ぶ必要がある
  
• MatchItパッケージにおける matchit関数の引数として設定できる
  ・非復元の場合・・・replace = FALSEと指定
  ・復元の場合・・・replace = TRUEと指定
  
• 2021年総選挙に出馬した下記の 5 人の候補者を例として男女のマッチングについて考える

• この 5 人を年齢だけを基準にしてマッチングする
  
• 年齢の近い男女をペアにする
  
• 年齢が20歳以上離れている場合はマッチングペアとして不適切だとする

復元によるマッチング

• 森下は 40歳の女性なので、ペアの候補は岸田、安倍、小泉の 3 人の男性
  
• 3 人の男性の中で、森下と年齢が近いのは小泉（40歳）
  →　「森下・小泉ペア」が可能

• 復元によるマッチングでは、一度選ばれても再度選ばれる
  →　池内（39歳）と小泉（40歳）の年齢が近い
  →　「池内・小泉ペア」が可能

• 2 組のペアが可能

非復元によるマッチング

• 森下は40歳の女性なので、ペアの候補は岸田、安倍、小泉の 3 人の男性
  
• 3 人の男性の中で、森下と年齢が近いのは小泉（40歳）
  →　「森下・小泉ペア」が可能

• 非復元によるマッチングでは、一度選ばれた人は再度選ばれない
  → 池内（39歳）と小泉（40歳）はペアになれない
  
• 森下（40歳）や池内（39歳）と岸田（64歳）や安倍（64歳）は 20 歳以上年齢が離れているので、ペアにはなれない
  
• 1 組のペアだけが可能

復元によるマッチング方法の長所と短所

• 復元によるマッチング方法では多くのペアができるので、データに含まれる情報が無駄になりにくく、分析結果の偏りが小さい・・・長所
  
• しかし、同一人物が複数のペアに登場するため、マッチングされたデータは独立ではない・・・短所
  →　復元によるマッチング方法を使う場合、マッチングの重みを考慮する必要があり、推定方法が複雑になる・・・・短所

復元と非復元によるマッチング、どちらが良いのか？ ・復元と非復元はどちらも一長一短あるので、長所と短所に配慮して慎重に決める

6.4.1 個体間の距離 (distance)

• マッチングの対象となる個体が似ていればペアを作る
  
• その基準が「個体間の距離」
  
• MatchItパッケージにおける matchit関数の引数 distance として設定できる

distance の引数一覧

glm	一般化線形モデル	generalize linear model
gam	一般化加法モデル	generalize additive model
rpart	分類木	classification tree
randomforest	ランダムフォレスト	random forest
nnet	ニューラルネットワーク	neural network
cbps	共変量バランシング傾向スコア	covariate blancing propensity score
bart	ベイズ的加法回帰木	BART: Bayesian additive regression trees
mahalanobis	マハラノビスの距離	Mahalanlobis distance:傾向スコアは計算されず

• glm がデフォルト
  
• randomforest の場合・・・distance = randomforestと指定
  
• それぞれの手法からモデル化した傾向スコアを「距離」として用いることができる
  
• 傾向スコアの値が値が近いもの同士をペアにする
  
• 個体 \(i\) と個体 \(j\) の距離・・・傾向スコアの差の絶対値 \(|e_i - e_j|\) によって測定

6.4.2 マッチングの方法

• MatchItパッケージにおける matchit関数の引数 method として設定できる

method の引数一覧

nearest	最近隣法マッチング	nearest neighbor matching
optimal	最適ペアマッチング	optimal pair matching
genetic	遺伝的マッチング	genetic matching
exact	厳密マッチング	exact matching
cem	単純化厳密マッチング	coarsened exact matching
subclass	階層	subclassification
full	最適フルマッチング	optimal full matching

最近隣法マッチング

• 処置群の個体に対して傾向スコアの値が近い個体を 1 つずつ順番にマッチングする
  
• 必ずしもデータ全体での距離が最小化するとは限らない

最適ペアマッチング

• マッチングされた全てのペアに関して、平均絶対距離 (average absolute distance) が最小化するようにマッチングする

遺伝的マッチング

• 遺伝的探索アルゴリズムに基づいて、それぞれの共変量に対してマッチング後のバランスが最適になるように重み付けする方法

厳密マッチング

• 共変量の値が同じ値のペアを作成する方法
  
• 実際のデータ解析ではあまり使わない

単純化厳密マッチング

• 厳密マッチングの代替法として提案された方法
  
• 連続変数を使った厳密マッチングでは、全く同じ値の個体同士を見つけることができないことがある
  
• そういう場合、連続変数を複数のカテゴリに分割して厳密マッチングを行う
  
• 例えば、「人口」という連続変数を「少ない」「中程度」「多い」などの分類する方法

傾向スコアをモデリングする際の注意 ・処置群と統制群における各々の共変量の標準化平均差の絶対値は 0.1 未満が推奨されている

傾向スコアマッチングによる ATT の推定

library(MatchIt)

mit_1 <- matchit(t1 ~ x1 + x2 + x3 + x4 + x5 + x6,
                 data = df11,
                 replace = TRUE,   # 復元のマッチングを指定
                 distance = "glm", # ロジステック回帰によって計算された傾向スコアを距離として使う  
                 method = "nearest") # 最近隣法マッチングを指定  

• マッチング後のデータフレームに df12 という名前を付ける

df12 <- match.data(mit_1)

• 小数点を統一してデータの中身を確認する

DT::datatable(df12)

• distance の値が傾向スコア：範囲は 0.11（最小値）から 0.81（最大値）

• マッチング後のデータフレーム df12 を使って ATT を推定する
  

• 復元 (replace = TRUE) のマッチングなので、weights == weights とマッチングの重みを設定する
  

• 理論上、傾向スコアマッチングによって共変量の影響による偏りを除去したので ATT を推定できる
  

• マッチング後のデータ df12 を用いて平均の差を求めれば良い

model1 <- lm(y3 ~ t1,
             data = df12,
             weights = weights)

• ATT の推定結果を表示する

tidy(model1)

# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    39.5       1.46     27.0  2.26e-106
2 t1              2.59      1.85      1.41 1.60e-  1

ATT の真値と推定値 (model1) ・ATT の真値 = 2.889
・ATT の推定値 = 2.59 (model1)

• ATT の真値 (2.899) と比較すると、過小推定 (2.59) されている（差は 0.309）

その理由：

• 傾向スコアマッチングの後でも共変量の影響が残っているため

解決策：

・傾向スコアマッチングの時だけでなく、解析モデルにも共変量 (\(x1, ..., x6\)) を取り込んだ共分散分析を行い ATT を推定する
・復元(replace = TRUE) のマッチングなので、weights == weights とマッチングの重みを設定する

model2 <- lm(y3 ~ t1 + x1 + x2 + x3 + x4 + x5 + x6,
             data = df12,
             weights = weights)

• ATT の推定結果を表示する

tidy(model2)

# A tibble: 8 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)     1.85     0.505      3.66 2.73e-  4
2 t1              2.92     0.403      7.25 1.29e- 12
3 x1              1.73     0.247      7.01 6.32e- 12
4 x2              1.10     0.248      4.42 1.18e-  5
5 x3             -2.12     0.244     -8.68 3.52e- 17
6 x4              7.85     0.255     30.7  5.20e-126
7 x5              9.99     0.248     40.3  3.96e-174
8 x6             11.3      0.241     46.7  8.91e-204

ATT の真値と推定値 (model2) ・ATT の真値 = 2.889
・ATT の推定値 = 2.92 (model2)

• ATT の真値 (2.899) がかなり正確に推定 (2.92) されている（差は 0.031）
  
• model1 よりも model2 の方が ATT の真値に近いATT 値を推定できることがわかる

ポイント ・共分散分析と比較すると、傾向スコアモデリングの方が、間違ったモデルを設定しても柔軟に対応できる
・傾向スコアマッチングの時だけでなく、解析モデルにも共変量 (\(x1, ..., x6\)) を取り込んだ共分散分析を行い ATT を推定すべき

クラスターに頑強な標準誤差を使う場合
・ 傾向スコアがモデルから計算されたことに起因する不確実性を標準誤差に反映させる方法
・ 標準誤差を推定する際、データがマッチングされたという性質を反映させるべきと主張する研究者もいる
・ どの個体がペアになっているのかという情報をクラスターとして、クラスターに頑強な標準誤差を使う
・ ここでは上で推定した model2 を事例としてクラスターに頑強な標準誤差を計算してみる
・ 2 つのRパッケージをロードする　　

library(lmtest)
library(sandwich) # 不均一分散に頑強な標準誤差を推定するためのパッケージ

・ cluster の設定方法：
・ 非復元 (replace = FALSE) なら cluster = ~ subclass
・ 復元 (replace = TRUE) なら cluster = ~ weights

coeftest(model2, 
         vcov. = vcovCL,
         cluster = ~ weights) 

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  1.84855    1.08109  1.7099    0.0878 .  
t1           2.91741    0.67454  4.3251 1.781e-05 ***
x1           1.72827    0.33283  5.1927 2.829e-07 ***
x2           1.09558    2.06205  0.5313    0.5954    
x3          -2.11607    3.03974 -0.6961    0.4866    
x4           7.84693    0.33587 23.3631 < 2.2e-16 ***
x5           9.99163    0.25037 39.9077 < 2.2e-16 ***
x6          11.27663    0.42177 26.7365 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

分析結果

	推定値	標準誤差
ATT の真値	2.889
通常の標準誤差	2.92	0.403
頑強な標準誤差	2.92	0.675

・ 点推定値は「通常の標準誤差」を使った場合も「頑強な標準誤差」を使ったばあいもどちらも 2.92
・ 「通常の標準誤差」を使った場合の標準誤差 (0.403)と比較すると、「通常の標準誤差」を使った場合の標準誤差 (0.675) は若干大きくなっていることがわかる
・ t1 の p-value = 1.781e-05 なので、処置効果 (ATT) はあるという結論が得られる

傾向スコアによるバランシング
・ 傾向スコアを使ってマッチングする際、処置群と統制群が「十分に似ているかどうか」を確認する必要がある
→　共変量の分布が処置群と統制群で十分にバランシングしているかどうか
・ ここでは mit_1 の結果のサマリーを出力すればバランシングの評価指標を表示できる

mit_1 <- matchit(t1 ~ x1 + x2 + x3 + x4 + x5 + x6,
                 data = df11,
                 replace = TRUE,   # 復元のマッチングを指定
                 distance = "glm", # ロジステック回帰によって計算された傾向スコアを距離として使う  
                 method = "nearest") # 最近隣法マッチングを指定  

・ Std.Mean Diff の値が 0 に近いほど、バランスがとれている
・ Var.Ratio の値が 1 に近いほど、バランスがとれている

summary(mit_1)

Call:
matchit(formula = t1 ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11, 
    method = "nearest", distance = "glm", replace = TRUE)

Summary of Balance for All Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.4371        0.3584          0.6022     0.9897    0.1755
x1              0.9655        0.9919         -0.0297     0.7966    0.0210
x2              0.9596        0.9810         -0.0217     0.8841    0.0200
x3              1.1451        0.9162          0.2196     1.0303    0.0616
x4              1.2246        0.8986          0.3429     0.9352    0.1017
x5              1.2874        0.7913          0.4963     1.0065    0.1461
x6              1.2896        0.8382          0.4860     0.9166    0.1376
         eCDF Max
distance   0.3232
x1         0.0540
x2         0.0524
x3         0.1002
x4         0.1797
x5         0.2387
x6         0.2212


Summary of Balance for Matched Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.4371        0.4369          0.0012     1.0028    0.0014
x1              0.9655        0.9759         -0.0117     0.7990    0.0241
x2              0.9596        0.9644         -0.0049     0.8047    0.0273
x3              1.1451        1.1281          0.0162     0.9851    0.0100
x4              1.2246        1.2465         -0.0231     0.8135    0.0241
x5              1.2874        1.2614          0.0260     1.0968    0.0121
x6              1.2896        1.3208         -0.0336     0.9959    0.0120
         eCDF Max Std. Pair Dist.
distance   0.0129          0.0076
x1         0.0668          1.2125
x2         0.0771          1.2002
x3         0.0386          1.0269
x4         0.0566          0.9933
x5         0.0540          0.6018
x6         0.0488          0.6997

Percent Balance Improvement:
         Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance            99.8       73.1      99.2     96.0
x1                  60.7        1.3     -14.6    -23.7
x2                  77.6      -76.3     -36.3    -47.2
x3                  92.6       49.7      83.8     61.5
x4                  93.3     -208.3      76.3     68.5
x5                  94.8    -1328.8      91.7     77.4
x6                  93.1       95.2      91.3     77.9

Sample Sizes:
              Control Treated
All            611.       389
Matched (ESS)  167.21     389
Matched        229.       389
Unmatched      382.         0
Discarded        0.         0

・ これだとわかりにくいので love plot を描いて、バランシングの改善を可視化してみる

plot(summary(mit_1), position=NULL)

・ ○ が傾向スコアマッチング前
・ ● が傾向スコアマッチング後
・ 横軸・・・標準化平均差 (Standardized Mean Differences)
・ 標準化平均差 = 処置群と統制群の分布の違い
・ この値が小さいほど望ましい
・ この値が 0 であることが理想
・ この値は 0.1 以下であることが望ましい
・ マッチング前（○）と比べると、マッチング後（●）の標準化平均差が小さくなっている
・ 傾向スコアマッチング後に、処置群と統制群の分布の違いが小さくなっている
→　傾向スコアマッチング後の（●）は全て 0.1 以下

→　ここでの傾向スコアマッチングでは、共変量の偏りが十分是正された

標準化平均差の絶対値が 0.1 以上になった時の対策 ・両群がバランスするように、傾向スコアのモデルに高次の項、非線形の変換、交差項などを含める

• 傾向スコアマッチングを使えば、ATT を推定できる
  
• 「処置の割り付けに影響を与えるのは、観測された共変量だけ」という仮定が満たされるなら、共分散分析を使って ATE を推定できる
  
• しかし、現実的にこの仮定が満たされることは難しい
  →　傾向スコアマッチングを使っても、ATE を推定できない

解決策：

• 傾向スコアを使った「層化解析法」と「重み付け法」

7. 層化解析法

7.1 標本調査における層化抽出法

無作為抽出法 (simple ramdom sampling)

• 全ての個体が同じ確率で母集団から標本に選ばれる抽出法
  「2.1 調査・観察データの場合」で紹介したケース
  
• 100人の有限母集団（男性50人、女性50人）から無作為に 50 人抽出することを考える 　　
• 最悪の場合・・・偶発的に男性50人、女性50人ずつ割り振られてしまうことがあり得る

• 「ランダムな割り付けは、自然には実現されない」（ローゼンバウム 2021）
  →　男女比にばらつきが出る　→　性別による差が生じてしまう
  →　この状況は交換可能ではない（＝グループを入れ替えられない）
  →　平均因果効果 (ATE) を推定できない

• 「自然に任せていたのでは公平な比較は望むべくもない」（ローゼンバウム 2021）

解決策 →　層化無作為抽出法

• ランダム化比較試験 (RCT) を実現するため、事前にセレクションバイアスが生じないように母集団をブロッキング（層化）する

層化無作為抽出法 (stratified ramdom sampling)

• 単純無作為抽出法 (simple random sampling) では、通常、母集団はひとつ
  
• 層化無作為抽出法では、母集団を複数の集団に分ける
  → 例：性別でブロッキングする（= 層化する） →　それぞれのブロック内でランダムに半数を処置群と統制群に割り付ける
  →　ブロックした要素（この場合は性別）に関して、処置群と統制群の違いがなくなる

• 層化無作為抽出法の事例:
  
• 例えば日本国民全体から 28,200人のサンプルを得るとする
  
• 無作為抽出法であれば、下図の上側のグラフのように、都道府県ごとにサンプル数 (N) に偏りがでる可能性がある

• しかし、ブロッキング（層化）によって 47 都道府県からそれぞれ N = 600人 無作為に選べば、都道府県ごとにサンプル数 (N)が揃うので都道府県ごとの偏りを防止ででき、推定の精度を良くできる

7.2 層化解析法

• 調査・観察研究における割り付けは無作為ではない
  →　処置群と統制群を単純比較しても因果効果は推定できない

解決策　→　交絡因子を使って層に分ける

→　層化に使った交絡因子の影響を除去できる
- 2021年総選挙に出馬した10人の候補者のリストを「性別」と「年齢」という二つの交絡因子で層化してみる

• 左側の名簿の情報を、50歳以上と50歳未満で層化

• 4 つの層に分類できる
  

• 同じ層内では「性別」は同じで、年齢もほぼ同じ
  →　層内での結果（得票率）が異なるのは、性別と年齢が理由ではない
  →　他の理由（例えば、処置＝世襲かどうか）のはず
  

• それぞれの層ごとに推定値を計算してみる
  

• ここでの知りたい「世襲候補者であることが得票率に与える因果効果」= 平均処置効果 (ATE) は次の式で表すことができる

\[\hat{τ}_{ATE}=\sum_{k=1}^K \frac{n_k}{N}[\bar{Y}_k(1) - \bar{Y}_k(0)]\]

K	層の数
N	標本サイズ
\(n_k\)	k 番目の層における標本サイズ
\(\bar{Y}_k(1)\)	k番目の層における処置群（世襲）の平均値
\(\bar{Y}_k(0)\)	k番目の層における統制群（非世襲）の平均値
\(\bar{Y}_k(1) - \bar{Y}_k(0)\)	k番目の層における平均処置効果 (ATE)

• 上図（「性別」と「年齢」で層化）を使って、平均処置効果 (ATE) を計算してみる

層 1 における ATE:

• 処置軍 (処置 = 1) における結果変数の平均値 = 79.17
  
• 統制軍 (処置 = 0) における結果変数の平均値 = 51.25
• ATE = 79.17 - 51.25 = 27.92

層 2 における ATE:

• 処置軍 (処置 = 1) における結果変数の平均値 = 58.06
  
• 統制軍 (処置 = 0) における結果変数の平均値 = 43.11
  
• ATE = 58.06 - 43.11 = 14.95

層 3 における ATE:

• 処置軍 (処置 = 1) における結果変数の平均値 = 41.72
  
• 統制軍 (処置 = 0) における結果変数の平均値 = (33.03+52.35)/2 = 42.69
  
• ATE = 41.72 - 42.69 = - 0.97

層 4 における ATE:

• 処置軍 (処置 = 1) における結果変数の平均値 = (80.67 + 69.72)/2 = 75.2
  
• 統制軍 (処置 = 0) における結果変数の平均値 = 63.5
  
• ATE = 75.2 - 63.5 = 11.7

全体の ATE:

\[\hat{τ}_{ATE}=\sum_{k=1}^K \frac{n_k}{N}[\bar{Y}_k(1) - \bar{Y}_k(0)] \\ = \frac{2}{10} ・ 27.92 + \frac{2}{10} ・ 14.95 + \frac{3}{10}・ (-0.97)+ \frac{3}{10}・ 11.7 = 11.79 \]

因果推論の妥当性

• ここでの因果推論が妥当かどうかは、年齢と性別以外の未観測の交絡因子がないという仮定が正しければという条件に依存する
  
• つまり「条件付き独立性（無交絡性）」がみたされること

\[{Y(1), Y(0)} ⊥ T|e(X)\]

• 「処置の割り付けに影響を与えるのは、観測された共変量だけ」ということ
  
• そのためには共変量 \(X\) は多変量である必要がある
  
• 得票率に影響を与える多変量として「年齢」と「性別」だけでは不十分
  
• 交絡因子として「候補者の経験」「所属政党」「学歴」「費やした選挙費用」「美顔度」など様々な要因が考えられる
  
• 上の 2 つのマッチングにおける得票率の差が大きく異なるのは (27.9 vs 3.4)、このような多くの交絡因子が存在していることの証左
  
• 多変量 \(X\) を条件としたときに、処置に割り当てられる確率として単変量の傾向スコアに情報を集約して因果効果を推定できる

層の数

• 層化解析を行うためには、データをいくつかの層に分ける必要がある
  
• 一般的には 5 〜 10 が妥当
  
• 標本サイズが大きければ、10 〜 20 が妥当
  
• R パッケージ MatchIt のデフォルトは 6
  
• 各層に含まれる処置群の数がほぼ同じになるよう設定されている
  →　あまり神経質になる必要はない

7.3 R による傾向スコア層化解析法：ATE の推定

• ここでは人工的に生成した次の 10 の変数を使う

y0t	潜在的結果変数: \(Y_i(0)\)	通常、観察不可能だが、ここでは観察可と仮定する
y1t	潜在的結果変数:\(Y_i(1)\)	通常、観察不可能だが、ここでは観察可と仮定する
\(y3\)	結果変数	観察可能
t1	処置変数（1=処置あり, 0=処置なし）	観察可能
x1	共変量	観察可能
x2	共変量	観察可能
x3	共変量	観察可能
x4	共変量	観察可能
x5	共変量	観察可能
x6	共変量	観察可能

• x1,....,x6 は平均1、標準偏差1の多変量対数正規分布に従う乱数

• 誤差項 \(ε_i\) は標準正規乱数である

• ここでは、高橋将宣著『統計的因果推論の理論と実装』(2022) のサポートページで公開されている「本書で使用するデータ（zip）」をダウンロードし data11.csv を使う

• 生成したデータを読み取る

df11 <- read_csv("data/takahashi/data11.csv")

DT::datatable(df11)

• \(Y_i(0)\) に関しては次のように設定してデータを生成してある

\[Y_i(0) = 1 + 1X_{1i} - 3X_{2i} + 5X_{3i} + 7X_{4i} + 9X_{5i} + 11X_{6i} + ε_i\]

• 試しに y0t を応答変数として重回帰分析をしてみる　　　

m0 <- lm(y0t ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11)

tidy(m0)

# A tibble: 7 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    0.993    0.0610      16.3 6.53e- 53
2 x1             0.982    0.0376      26.1 8.27e-115
3 x2            -2.99     0.0396     -75.6 0        
4 x3             5.04     0.0396     127.  0        
5 x4             6.96     0.0418     167.  0        
6 x5             8.99     0.0400     225.  0        
7 x6            11.0      0.0380     291.  0        

• \(Y_i(1)\) に関しては次のように設定してデータを生成してある
  \[Y_i(1) = 5 + 2X_{1i} + 4X_{2i} - 6X_{3i} + 8X_{4i} + 10X_{5i} + 12X_{6i} + ε_i\]

• 試しに y1t を応答変数として重回帰分析をしてみる　　　

m1 <- lm(y1t ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11)

tidy(m1)

# A tibble: 7 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)     4.99    0.0610      81.8 0        
2 x1              1.98    0.0376      52.7 6.11e-290
3 x2              4.01    0.0396     101.  0        
4 x3             -5.96    0.0396    -150.  0        
5 x4              7.96    0.0418     191.  0        
6 x5              9.99    0.0400     250.  0        
7 x6             12.0     0.0380     317.  0        

処置 t1 の割り付けルール

• \(med(Y_i(0))\) は \(Y_i(0)\) の中央値
  
• \(u_{1i}\) と \(u_{2i}\) は区間 \([0, 1]\) の一様乱数

\[T_i = 1 \\if Y_i(0) > med(Y_i(0)) \& u_{1i} ≦ 0.50 \\or Y_i(0) ≦ med(Y_i(0)) \& u_{2i} > 0.75\]

• \(Y_i(0)\) が中央値より大きいとき 50% の確率で \(T_i = 1\)
  
• \(Y_i(0)\) が中央値以下のとき 25% の確率で \(T_i = 1\)

\[T_i = 0 \\if Y_i(0) > med(Y_i(0) \& u_{1i} > 0.50 \\or Y_i(0) ≦ med(Y_i(0)) \& u_{2i} ≦ 0.75\]

• \(Y_i(0)\) が中央値より大きいとき 50% の確率で \(T_i = 0\)
  
• \(Y_i(0)\) が中央値以下のとき 75% の確率で \(T_i = 0\)

平均処置効果 (ATE) の真値

mean(df11$y1t) - mean(df11$y0t)

[1] 3.755947

人工的に設定した ATE の真値 (data11.csv) ・ATE の真値 = 3.756

実際に私たちが観測できるのは次の変数だけ

\(y3\)	結果変数	観察可能
t1	処置変数（1: 処置あり、0: 処置なし）	観察可能
x1	共変量	観察可能
x2	共変量	観察可能
x3	共変量	観察可能
x4	共変量	観察可能
x5	共変量	観察可能
x6	共変量	観察可能

• 傾向スコアをモデル化する

library(MatchIt)

sub <- 5 # 層の数を指定  
m.out2 <- matchit(t1 ~ x1 + x2 + x3 + x4 + x5 + x6,
                 data = df11,
                 method = "subclass",
                 subclass = sub,
                 estimand = "ATE",
                 min.n = 2)
m.out2

A matchit object
 - method: Subclassification (5 subclasses)
 - distance: Propensity score
             - estimated with logistic regression
 - number of obs.: 1000 (original), 1000 (matched)
 - target estimand: ATE
 - covariates: x1, x2, x3, x4, x5, x6

m.data2 <- match.data(m.out2)

names(df11)

 [1] "y0t" "y1t" "y3"  "t1"  "x1"  "x2"  "x3"  "x4"  "x5"  "x6" 

names(m.data2)

 [1] "y0t"      "y1t"      "y3"       "t1"       "x1"       "x2"      
 [7] "x3"       "x4"       "x5"       "x6"       "distance" "weights" 
[13] "subclass"

• もともとのデータ (df11) と比較すると、m.data2 では新たに 3 つの変数 (distance, weights, subclass) がつくられていることがわかる
  
• m.data2 の中を確認する　　

DT::datatable(m.data2)

• データを 5 つの層に分類できた

7.3.1 単純な重回帰分析による推定

model4 <- lm(y3 ~ t1 + x1 + x2 + x3 + x4 + x5 + x6, 　
               data = df11)　
tidy(model4)

# A tibble: 8 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    1.05      0.309      3.41 6.78e-  4
2 t1             3.47      0.327     10.6  4.24e- 25
3 x1             1.39      0.184      7.57 8.77e- 14
4 x2            -0.292     0.194     -1.50 1.33e-  1
5 x3             0.454     0.194      2.34 1.97e-  2
6 x4             7.47      0.205     36.5  8.05e-186
7 x5             9.66      0.197     49.1  8.91e-268
8 x6            11.2       0.187     59.5  0        

• seshu_dummy の係数・・・ 3.47
  
• 全ての共変量の値を平均値に固定した時、t = 1 であれば y3 が平均して 3.47%ポイント上がる
  
• p-value = 4.24e-25 →　統計的に有意

結論： →　t1 の推定値は 3.47

7.3.2 ATE の推定（層内で共分散分析を行わない場合）

• ここでやるべきことは次の作業
  

1. 5 つの subclass ごとに t1 = 1と t1 = 0 それぞれの y3 の平均値を求める
   
2. 平均値の差を求める（これが ATE）

subclass = 1 の ATE

df01 <- m.data2 |> 
  filter(subclass == 1)

sub11_ave <- mean(df01$y3[df01$t1==1]) # 処置群 (t1=1)での y3 の平均値
sub10_ave <- mean(df01$y3[df01$t1==0]) # 統制群 (t1=0)での y3 の平均値

sub11_ave - sub10_ave  # ATE  

[1] 8.86462

subclass = 2 の ATE

df02 <- m.data2 |> 
  filter(subclass == 2)

sub21_ave <- mean(df02$y3[df02$t1==1])
sub20_ave <- mean(df02$y3[df02$t1==0])

sub21_ave - sub20_ave

[1] 3.321628

subclass = 3 の ATE

df03 <- m.data2 |> 
  filter(subclass == 3)

sub31_ave <- mean(df03$y3[df03$t1==1])
sub30_ave <- mean(df03$y3[df03$t1==0])

sub31_ave - sub30_ave

[1] 7.223671

subclass = 4 の ATE

df04 <- m.data2 |> 
  filter(subclass == 4)

sub41_ave <- mean(df04$y3[df04$t1==1])
sub40_ave <- mean(df04$y3[df04$t1==0])

sub41_ave - sub40_ave

[1] 2.714279

subclass = 5 の ATE

df05 <- m.data2 |> 
  filter(subclass == 5)

sub51_ave <- mean(df05$y3[df05$t1==1])
sub50_ave <- mean(df05$y3[df05$t1==0])

sub51_ave - sub50_ave

[1] -1.364288

• 全体の ATE を計算する

K	層の数 = 5
N	標本サイズ = 1000
\(n_k\)	k 番目の層における標本サイズ = 200
\(\bar{Y}_k(1)\)	k番目の層における処置群（世襲）の平均値
\(\bar{Y}_k(0)\)	k番目の層における統制群（非世襲）の平均値
\(\bar{Y}_k(1) - \bar{Y}_k(0)\)	k番目の層における平均処置効果 (ATE)

• 層ごとの標本サイズ \(n_k\) を求める

table(m.data2$subclass)

  1   2   3   4   5 
200 200 200 200 200 

\[\hat{τ}_{ATE}=\sum_{k=1}^K \frac{n_k}{N}[\bar{Y}_k(1) - \bar{Y}_k(0)] \\ = \frac{200}{1000} ・ 8.86462 + \frac{200}{1000} ・ 3.321628 + \frac{200}{1000}・7.223671 \\+ \frac{200}{1000}・2.714279 + \frac{200}{1000}・(-1.364288) \\ = 4.15 \]

200*(8.86462 + 3.321628 + 7.223671 + 2.714279 - 1.364288)/1000

[1] 4.151982

7.3.2 ATE の推定（層内で共分散分析を行う場合）

• 傾向スコア層化解析による ATE の推定
  
• 分析に必要な 2 つの R パッケージをロードする

library(lmtest)
library(sandwich)

• 5 つの層ごとに分析する →　その結果を統合する
  
• 結果を統合する際に使う 4 つの空の入れ物 (= 変数) を作る: psp, psvar, nps, robustvar

n <- nrow(df11)    # データに含まれる標本サイズ (n = 1000) を設定
psp <- NULL　　　　# t1 の偏回帰係数を入れる
ps_var <- NULL      # t1 の標準誤差を二乗した値（＝分散）を入れる   
nps <- NULL　　　　# 層ごとの標本サイズ 
robust_var <- NULL  # 各層内でのクラスターに頑強な分散を入れる  

• 同じ作業を 5 回繰り返すので、for ループを使って計算を自動化する
  
• 傾向スコアをモデル化した（= 層化した）データフレーム m.data2 の 1 から 5 までの層化番号を j と指定する
  
• 観測された結果変数 y3 を 1 つの処置変数 t1 と 6 つの共変量 x1, x2, x3, x4, x5, x6 で回帰する
  
• j を 1 〜 sub まで (1 から 5 まで) 変化させる

for(j in 1:sub){
  data_ps <- m.data2[m.data2$subclass == j,]    # j 番目の層のデータを取り出して data_ps と名前を付ける  
  model4 <- lm(y3 ~ t1 + x1 + x2 + x3 + x4 + x5 + x6, 　
               data = data_ps)　　　　          # data_ps を使って j 番目の層内で共分散分析を実行
  psp[j] <- summary(model4)$coefficients[2, 1] # t1 の偏回帰係数を取り出して psp に入れる  
  ps_var[j] <- summary(model4)$coefficients[2, 2]^2 # t1 の標準誤差を二乗して (=分散) 取り出して psvar に入れる 
  nps[j] <- nrow(data_ps) # data_ps の標本サイズ (N) を入れる  
  robust_var[j] <- coeftest(model4, # 各層内でのクラスターに頑強な分散  
                           vcov. = vcovCL,
                           cluster = ~weights)[2,2]
}

• 各層内での平均処置効果 psp は次のとおり

psp

[1] 8.2114114 4.4809301 3.8157096 2.1218324 0.1144802

• 層ごとの標本サイズ nps は

nps

[1] 200 200 200 200 200

• これらの情報から 5 つの層で得られた平均処置効果 psp を統合する

\[\hat{τ}_{ATE}=\sum_{k=1}^K \frac{n_k}{N}[\bar{Y}_k(1) - \bar{Y}_k(0)] \\ = \frac{200}{1000} ・ 8.2114114 + \frac{200}{1000} ・ 4.4809301 + \frac{200}{1000}・3.8157096 \\ + \frac{200}{1000}・2.1218324 + \frac{200}{1000}・(0.1144802) \\ = 3.748873 \]

200*(8.2114114 + 4.4809301 + 3.8157096 + 2.1218324 + 0.1144802)/1000

[1] 3.748873

• 次のように入力しても同様の結果が得られる

tau_hat <- sum((nps/n) * psp)
tau_hat

[1] 3.748873

統計的有意性をチェック

• 統計的有意性を確認するために必要なのは次の 2 つの情報
  1. ATE の推定値: tau_hat
  2. 統合したクラスターに頑強な標準誤差

• ATE の推定値: tau_hat は上で得られたので、統合したクラスターに頑強な標準誤差を計算する

• 各層内での標準誤差の 2 乗（＝分散）ps_var は次のとおり

ps_var

[1] 0.4013066 0.5674602 0.4879605 0.4931660 0.4106398

• ps_var を使って統合した標準誤差を計算する

var_tau <- sum((nps/n)^2 * ps_var)
sqrt(var_tau)

[1] 0.3072805

• 各層内でのクラスターに頑強な分散 robust_var は次のとおり

robust_var

[1] 0.2943865 0.4497133 0.4868898 0.6962178 0.8831411

• robust_var を使って統合したクラスターに頑強な標準誤差を計算する

robust_var_tau <- sum((nps/n)^2 * robust_var)
robust_se_tau <- sqrt(robust_var_tau)
robust_se_tau

[1] 0.3352819

• tau_hat と robust_se_tau を使って 95% 信頼区間を計算する

ci_95 <- qt(0.975, n - 8)
tau_hat + ci_95 * robust_se_tau

[1] 4.406816

tau_hat - ci_95 * robust_se_tau

[1] 3.09093

• 95%信頼区間は 3.0 〜 4.4 で 0 をまたがない
  
• 母集団でも帰無仮説は棄却され、5% で統計的に有意

ATE の真値と推定値 (data11.csv) ・ATE の真値 = 3.756
・単純な重回帰分析の推定値 = 3.47（p-value = 4.24e-25）
・共分散分析なしのATE の推定値 = 4.15
・共分散分析ありのATE の推定値 = 3.75（統計的に有意）
→　層内で共分散分析を行った時のATE の推定値の方が、より ATE の真値 (3.765) に近い推定値 (3.75) を得られる

7.4 傾向スコアによるバランシング評価

・ 傾向スコアを使ってマッチングする際、処置群と統制群が「十分に似ているかどうか」を確認する必要がある
→　共変量の分布が処置群と統制群で十分にバランシングしているかどうか

・ ここでは傾向スコアをモデル化した m.out2 の結果のサマリーを出力してバランシングの評価指標を表示してみる

sub <- 5 # 層の数を指定  
m.out2 <- matchit(t1 ~ x1 + x2 + x3 + x4 + x5 + x6,
                 data = df11,
                 method = "subclass",
                 subclass = sub,
                 estimand = "ATE",
                 min.n = 2)
m.out2

A matchit object
 - method: Subclassification (5 subclasses)
 - distance: Propensity score
             - estimated with logistic regression
 - number of obs.: 1000 (original), 1000 (matched)
 - target estimand: ATE
 - covariates: x1, x2, x3, x4, x5, x6

・ Std.Mean Diff の値が 0 に近いほど、バランスがとれている
・ Var.Ratio の値が 1 に近いほど、バランスがとれている

summary(m.out2)

Call:
matchit(formula = t1 ~ x1 + x2 + x3 + x4 + x5 + x6, data = df11, 
    method = "subclass", estimand = "ATE", subclass = sub, min.n = 2)

Summary of Balance for All Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.4371        0.3584          0.6007     0.9897    0.1755
x1              0.9655        0.9919         -0.0280     0.7966    0.0210
x2              0.9596        0.9810         -0.0211     0.8841    0.0200
x3              1.1451        0.9162          0.2213     1.0303    0.0616
x4              1.2246        0.8986          0.3371     0.9352    0.1017
x5              1.2874        0.7913          0.4971     1.0065    0.1461
x6              1.2896        0.8382          0.4753     0.9166    0.1376
         eCDF Max
distance   0.3232
x1         0.0540
x2         0.0524
x3         0.1002
x4         0.1797
x5         0.2387
x6         0.2212

Summary of Balance Across Subclasses
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.3902        0.3882          0.0148     1.0303    0.0095
x1              0.9917        0.9727          0.0202     0.7583    0.0190
x2              0.9844        0.9713          0.0129     0.9283    0.0213
x3              1.0157        1.0030          0.0123     1.0600    0.0107
x4              1.0413        1.0228          0.0191     0.9487    0.0205
x5              0.9994        0.9760          0.0235     1.0893    0.0100
x6              1.0112        1.0121         -0.0009     0.9771    0.0083
         eCDF Max
distance   0.0339
x1         0.0536
x2         0.0474
x3         0.0408
x4         0.0553
x5         0.0324
x6         0.0328

Percent Balance Improvement:
         Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance            97.5       -4.1      94.6     89.5
x1                  27.9        4.8       9.5      0.8
x2                  38.6       -5.0      -6.2      9.6
x3                  94.4       -2.9      82.6     59.3
x4                  94.3       -1.4      79.8     69.2
x5                  95.3       -8.2      93.2     86.4
x6                  99.8       -6.6      94.0     85.1

Sample Sizes:
              Control Treated
All            611.       389
Matched (ESS)  570.38     325
Matched        611.       389
Unmatched        0.         0
Discarded        0.         0

・ これだとわかりにくいので love plot を描いて、バランシングの改善を可視化してみる

plot(summary(m.out2), position=NULL, xlim = c(-0.1, 0.6))

・ ○ が傾向スコアマッチング前
・ ● が傾向スコアマッチング後
・ 横軸・・・標準化平均差 (Standardized Mean Differences)
・ 標準化平均差 = 処置群と統制群の分布の違い
・ この値が小さいほど望ましい
・ この値が 0 であることが理想
・ この値は 0.1 以下であることが望ましい
・ マッチング前（○）と比べると、マッチング後（●）の標準化平均差が小さくなっている
・ 傾向スコアマッチング後に、処置群と統制群の分布の違いが小さくなっている
→　傾向スコアマッチング後の（●）は全て 0.1 以下

→　ここでの傾向スコアマッチングでは、共変量の偏りが十分是正された

8. 世襲候補者は選挙で有利なのか？

• ここでは、処置（世襲候補者か否か）が得票率に与える処置効果を傾向スコアマッチングを使って分析してみる
  

• 衆議院議員総選挙の得票データ hr96-21.csv をダウンロード
  

• 1996年から2021年までの総選挙データ(hr96-21.csv) を読み込む

df <- read_csv("data/hr96-21.csv", na = ".")

• 変数を確認

names(df)

 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"

• 二値変数を作成してデータフレームに追加する

df <- df %>% 
  mutate(male = ifelse(gender == "male", 1, 0), # 男性なら 1、女性なら 0 の変数 male
         inc = ifelse(status == 1, 1, 0),　　 　 # 現職なら 1、それ以外なら 0 の変数 inc
         ldp = ifelse(seito == "自民", 1, 0))# 選挙費用が平均以上なら 1、平均以下なら 0 の変数 more_exp

• 必要な変数だけを抜き出す
  
• 2009年総選挙データを使う

df1 <- df %>% 
  dplyr::filter(year == 2009) %>%  # 2009年総選挙データだけを抜き出す  
  dplyr::select(voteshare, seshu_dummy, male, inc, nocand, previous, age, ldp) # 必要な変数を指定

DT::datatable(df1)

• ここで使う変数は次のとおり

変数名	変数の役割	説明
voteshare	結果変数	候補者が得た得票率 (%)
seshu_dummy	処置効果	候補者が世襲なら 1、それ以外なら 0
male	共変量	候補者が男性なら 1、女性なら 0
inc	共変量	候補者が現職（小選挙区当選者）なら 1、そうでなければ 0
nocand	共変量	小選挙区内における立候補者の数
previous	共変量	当選回数
age	共変量	候補者の年齢
ldp	共変量	候補者が自民党所属なら 1、それ以外なら 0

• 傾向スコアをモデル化する

library(MatchIt)

• 欠損値を取り除く

df1 <- na.omit(df1)

sub <- 5 # 層の数を指定  
m.out3 <- matchit(seshu_dummy ~ male + inc + nocand + previous + age + ldp,
                 data = df1,
                 method = "subclass",
                 subclass = sub,
                 estimand = "ATE",
                 min.n = 2)
m.out3

A matchit object
 - method: Subclassification (5 subclasses)
 - distance: Propensity score
             - estimated with logistic regression
 - number of obs.: 1135 (original), 1135 (matched)
 - target estimand: ATE
 - covariates: male, inc, nocand, previous, age, ldp

m.data3 <- match.data(m.out3)
names(df1)

[1] "voteshare"   "seshu_dummy" "male"        "inc"         "nocand"     
[6] "previous"    "age"         "ldp"        

names(m.data3)

 [1] "voteshare"   "seshu_dummy" "male"        "inc"         "nocand"     
 [6] "previous"    "age"         "ldp"         "distance"    "weights"    
[11] "subclass"   

• もともとのデータ (df1) と比較すると、m.data3 では新たに 3 つの変数 (distance, weights, subclass) がつくられていることがわかる
  
• m.data3 の中を確認する　

DT::datatable(m.data3)

• データを 5 つの層に分類できた

8.1 単純な重回帰分析による推定

model5 <- lm(voteshare ~ seshu_dummy + male + inc + nocand + previous + age + ldp, 　
               data = df1)
tidy(model5)

# A tibble: 8 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   32.1      3.41        9.40 2.84e-20
2 seshu_dummy    2.79     1.80        1.55 1.21e- 1
3 male           3.90     1.37        2.84 4.56e- 3
4 inc           20.9      1.67       12.5  7.79e-34
5 nocand        -3.02     0.518      -5.82 7.46e- 9
6 previous       2.60     0.318       8.16 8.58e-16
7 age           -0.146    0.0509     -2.87 4.19e- 3
8 ldp           -2.74     1.56       -1.75 7.98e- 2

• seshu_dummy の係数・・・ 2.79
  
• 全ての共変量の値を平均値に固定した時、候補者が世襲であれば候補者の得票率が平均して 2.79%ポイント上がる
  
• p-value = 0.121 →　統計的に有意ではない

結論： →　候補者が世襲かどうかということと得票率は無関係

8.2 ATE の推定（層内で共分散分析を行う場合）

• 傾向スコア層化解析による ATE の推定
  
• 分析に必要な 2 つの R パッケージをロードする

library(lmtest)
library(sandwich)

• 5 つの層ごとに分析する →　その結果を統合する
  
• 結果を統合する際に使う 4 つの空の入れ物 (= 変数) を作る: psp, psvar, nps, robustvar

n <- nrow(df1)    # データに含まれる標本サイズ (n = 1294) を設定
psp <- NULL　　　　# t1 の偏回帰係数
ps_var <- NULL      # t1 の標準誤差を二乗した値（＝分散）   
nps <- NULL　　　　# 層ごとの標本サイズ 
robust_var <- NULL  # 各層内でのクラスターに頑強な分散  

• 同じ作業を 5 回繰り返すので、for ループを使って計算を自動化する
  
• 傾向スコアをモデル化した（= 層化した）データフレーム m.data3 の 1 から 5 までの層化番号を j と指定する
  
• 観測された結果変数 voteshare を 1 つの処置変数 female と 8 つの共変量 seshu_dummy, inc, nocand, previous, age, exp, male, ldp で回帰する
  
• j を 1 〜 sub まで (1 から 5 まで) 変化させる

for(j in 1:sub){
  data_ps <- m.data3[m.data3$subclass == j,]    # j 番目の層のデータを取り出して data_ps と名前を付ける  
  model5 <- lm(voteshare ~ seshu_dummy + male  + inc + nocand + previous + age + ldp, 　
               data = df1)　　　　          # data_ps を使って j 番目の層内で共分散分析を実行
  psp[j] <- summary(model5)$coefficients[2, 1] # female  の偏回帰係数を取り出して psp に入れる  
  ps_var[j] <- summary(model5)$coefficients[2, 2]^2 # female の標準誤差を 2 乗して (=分散) 取り出して ps_var に入れる 
  nps[j] <- nrow(data_ps) # data_ps の標本サイズ (N) を入れる  
  robust_var[j] <- coeftest(model5, # 各層内でのクラスターに頑強な分散  
                           vcov. = vcovCL)[2,2]
}

• 各層内での平均処置効果 psp は次のとおり

psp

[1] 2.793121 2.793121 2.793121 2.793121 2.793121

• 層ごとの標本サイズ nps は

nps

[1] 222 233 226 227 227

• これらの情報から 5 つの層で得られた平均処置効果 psp を統合する

\[\hat{τ}_{ATE}=\sum_{k=1}^K \frac{n_k}{N}[\bar{Y}_k(1) - \bar{Y}_k(0)] \\ = \frac{222}{1294} ・ 2.793121 + \frac{233}{1294} ・ 2.793121 + \frac{226}{1294}・2.793121 \\ + \frac{227}{1294}・2.793121 + \frac{227}{1294}・2.793121 \\ = 2.793121 \]

共分散分析ありの ATE の推定値

tau_hat <- sum((nps/n) * psp)
tau_hat

[1] 2.793121

統計的有意性をチェック

• 統計的有意性を確認するために必要なのは次の 2 つの情報
  1. ATE の推定値: tau_hat
  2. 統合したクラスターに頑強な標準誤差

• ATE の推定値: tau_hat は上で得られたので、統合したクラスターに頑強な標準誤差を計算する

• 各層内での標準誤差の 2 乗（＝分散）ps_var は次のとおり

ps_var

[1] 3.245972 3.245972 3.245972 3.245972 3.245972

• ps_var を使って統合した標準誤差を計算する

var_tau <- sum((nps/n)^2 * ps_var)
sqrt(var_tau)

[1] 0.805823

• 各層内でのクラスターに頑強な分散 robust_var は次のとおり

robust_var

[1] 1.589279 1.589279 1.589279 1.589279 1.589279

• robust_var を使って統合したクラスターに頑強な標準誤差を計算する

robust_var_tau <- sum((nps/n)^2 * robust_var)
robust_se_tau <- sqrt(robust_var_tau)
robust_se_tau

[1] 0.5638548

• tau_hat と robust_se_tau を使って 95% 信頼区間を計算する

ci_95 <- qt(0.975, n - 8) # パラメータの数は 8  
tau_hat + ci_95 * robust_se_tau

[1] 3.899444

tau_hat - ci_95 * robust_se_tau

[1] 1.686797

• 95%信頼区間は 1.686797 〜 3.899444 で 0 をまたがない
  
• 母集団では帰無仮説は棄却され、5% で統計的に有意

ATE の真値と推定値 (2009年総選挙：世襲議員の処置効果) ・ATE の真値 = 不明
・単純な重回帰分析の推定値 = 2.79（統計的に有意ではない）
・共分散分析ありのATE の推定値 = 2.79（統計的に有意）

8.3 傾向スコアによるバランシング評価

・ 傾向スコアを使ってマッチングする際、処置群と統制群が「十分に似ているかどうか」を確認する必要がある
→　共変量の分布が処置群と統制群で十分にバランシングしているかどうか

・ ここでは傾向スコアをモデル化した m.out2 の結果のサマリーを出力してバランシングの評価指標を表示してみる

sub <- 5 # 層の数を指定  
m.out3 <- matchit(seshu_dummy ~ male + inc + nocand + previous + age + ldp,
                 data = df1,
                 method = "subclass",
                 subclass = sub,
                 estimand = "ATE",
                 min.n = 2)
m.out3

A matchit object
 - method: Subclassification (5 subclasses)
 - distance: Propensity score
             - estimated with logistic regression
 - number of obs.: 1135 (original), 1135 (matched)
 - target estimand: ATE
 - covariates: male, inc, nocand, previous, age, ldp

・ Std.Mean Diff の値が 0 に近いほど、バランスがとれている
・ Var.Ratio の値が 1 に近いほど、バランスがとれている

summary(m.out3)

Call:
matchit(formula = seshu_dummy ~ male + inc + nocand + previous + 
    age + ldp, data = df1, method = "subclass", estimand = "ATE", 
    subclass = sub, min.n = 2)

Summary of Balance for All Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.3881        0.0805          1.6100     3.0786    0.4162
male            0.9470        0.8285          0.3820          .    0.1185
inc             0.8864        0.2772          1.5701          .    0.6092
nocand          3.7197        4.0379         -0.3477     0.7583    0.0436
previous        4.1818        0.9521          1.2347     2.8052    0.1899
age            54.8712       49.4666          0.4949     0.9193    0.0958
ldp             0.7576        0.1884          1.3874          .    0.5691
         eCDF Max
distance   0.6626
male       0.1185
inc        0.6092
nocand     0.1660
previous   0.5743
age        0.2058
ldp        0.5691

Summary of Balance Across Subclasses
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.1345        0.1067          0.1458     1.6670    0.0607
male            0.9347        0.8414          0.3008          .    0.0933
inc             0.3565        0.3468          0.0250          .    0.0097
nocand          4.1982        4.0057          0.2104     1.2364    0.0404
previous        1.9766        1.2076          0.2940     1.4411    0.0388
age            49.5520       49.9827         -0.0394     0.7247    0.0383
ldp             0.2554        0.2520          0.0083          .    0.0034
         eCDF Max
distance   0.2720
male       0.0933
inc        0.0097
nocand     0.1759
previous   0.2328
age        0.1309
ldp        0.0034

Percent Balance Improvement:
         Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance            90.9       45.9      85.4     58.9
male                21.3          .      21.3     21.3
inc                 98.4          .      98.4     98.4
nocand              39.5      -63.1       7.3     -5.9
previous            76.2       48.6      79.6     59.5
age                 92.0       21.2      60.0     36.4
ldp                 99.4          .      99.4     99.4

Sample Sizes:
              Control Treated
All           1003.    132.  
Matched (ESS)  956.82   20.43
Matched       1003.    132.  
Unmatched        0.      0.  
Discarded        0.      0.  

・ これだとわかりにくいので love plot を描いて、バランシングの改善を可視化してみる

plot(summary(m.out3), position=NULL, xlim = c(-0.5, 1.7))

・ ○ が傾向スコアマッチング前
・ ● が傾向スコアマッチング後
・ 横軸・・・標準化平均差 (Standardized Mean Differences)
・ 標準化平均差 = 処置群と統制群の分布の違い
・ この値が小さいほど望ましい
・ この値が 0 であることが理想
・ この値は 0.1 以下であることが望ましい
・ マッチング前（○）と比べると、マッチング後（●）の標準化平均差が小さくなっている
・ 傾向スコアマッチング後に、処置群と統制群の分布の違いが小さくなっている
→　傾向スコアマッチング後のバランス（●）は大幅に改善

→　ここでの傾向スコアマッチングでは、共変量の偏りが十分是正された

9. Exercise

Question:

「8. 世襲候補者は選挙で有利なのか？」を参考に、小泉純一郎総理が主導して大勝利をもたらした 2005 の「郵政選挙」において、世襲候補者は得票において有利だといえたのかどうか、傾向スコアマッチングを使って調べたい

Q1: 2005年の「郵政選挙」と民主党による政権交代が実現した 2009 年総選挙と比較すると、どちらの総選挙において世襲候補者がより有意だといえるか？

Q2: 2005年の「郵政選挙」と2009年の「政権交代選挙」において世襲候補者の優位性に差があるとすれば、それはなぜであろうか？分析結果から推定できる範囲であなたの考えを書きなさい

Reference

• 高橋将宣『統計的因果推論の理論と実装』共立出版 2022年
• Paul Rosenbaum, Observation and Experiment: An Introduction to Causal Inference, 2019
• ローゼンバウム『統計的因果推論入門ー観察研究とランダム化実験』共立出版、2021年
• Rosenbaum and Rubin, The central role of the propensity score in observational studies for causal effects, Biometrika, 70, 1, pp.41-55, 1983
• 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
• 浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
• 浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年
• Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
• Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017