匹配（Matching）是一種解決由自選擇（self-selection）導致內生性問題的方法。關於自選擇問題的成因，我們在理論篇一個隱藏的魔鬼——內生性中已經給出了詳細的解釋，這裡再簡單回顧一下。
自選擇問題的原因是研究中樣本的選擇不是隨機的，而是包含了個體自主選擇的結果。這就導致了實驗組和控制組之間不平衡，最後的結果存在誤差。
比如，研究是否參加某培訓對未來工資的影響，雖然參加培訓會增強工作能力，但在未參加者可能是因為本來自身能力就強而選擇的不參加，而能力和培訓都是影響未來工資的因素，因此直接比較參加者和未參加者並不合理。
從整個系統的角度看，這個自主選擇的過程說明是否選擇參加培訓是內生變量，忽視這個過程導致內生性問題。
正式的論述，可以回看之前的證明：
用
對隨機實驗而言，實驗組和控制組來自同一個總體，我們感興趣的平均處理效應是
進而有
此時可以用實驗組和控制組的數據分別估計
對自選擇模型而言，實驗組和控制組實際上來自於不同的總體，我們要估計的是實驗組的平均處理效應：
此時我們無法使用觀測組來估計
用匹配去創造平衡顧名思義，匹配就是為每個參與者尋找有著相同可觀測特點的未參與者（亦可為未參與者匹配參與者），使得他們之間唯一的差別只有是否參與實驗，這就解決了實驗組和控制組之間不平衡的問題。
隨機實驗是通過隨機性保證在樣本容量充足時實驗組與對照組的平衡性，包括所有可見和不可見因素；而匹配則是要人為去創造這種平衡。所以在使用匹配法時，我們假設所有的影響因素都是可見的。
正式的，協變量
其中，
最理想的情況是精確匹配，就是為每一個參與者對應選出除了是否參與培訓外，其他方面完全一樣的未參與者，但是這個方法最大的問題是，通常影響結果的可觀測變量會非常多，所以配對非常難找，這鐘情況被稱為「維度詛咒（curse of dimensionality）」。
傾向得分匹配（propensity score matching）是一種常用的降維方法。
定義每個個體「參與實驗的概率（probability of participation）」為傾向得分（propensity score），即
Rosenbaum 等人證明了
所以有
因此，依據傾向得分進行匹配是一個簡單易行的選擇。
例子：培訓對工資的影響說到傾向得分匹配，有兩篇不得不提的經典：Evaluating the econometric evaluation of training programs with Experimental data 和 Causal effects in nonexperimental studies: reevaluating the evaluation of training programs. Lalonde 和 Dehejia，Wahba 先後研究對比了匹配和隨機試驗的效果。下面我們也來體驗一下匹配的魅力。
這裡只收集到了部分數據，包含通過隨機實驗得到的實驗組和控制組數據，還有一份觀測組數據。簡單整理後數據的基本情況如下：
import pandas as pd
import numpy as np
import math
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.api as sm_api
from matplotlib import pyplot as plt
import seaborn as sns
import scipy
from sklearn.neighbors import NearestNeighbors
# 有兩個數據集
Data1 = pd.read_stata(r'matching_earnings.dta')
Data2 = pd.read_stata(r'ldw_exper.dta')
# 統一列名
Data2.columns = ['TREAT', 'AGE', 'EDUC', 'BLACK', 'HISP', 'MARR', 'NODEGREE', 'RE74', 'RE75', 'RE78', 'U74', 'U75']
# 將 Data2 的 RExx 數量級恢復到 Data1 級別，即 * 1000
for i in ['RE74', 'RE75', 'RE78']:
Data2[i] *= 1000
# 復現 Dehejia，Wahba 文章表 1
# 按組分離數據：
Data_RE74t = Data1.loc[Data1['TREAT'] == 1] # 實驗數據，treated
Data_PSID1 = Data1.loc[Data1['TREAT'] == 0] # 觀測數據，PSID-1
# Data2 中：TREAT == 1 為 RE74 subset treated，檢驗過和 Data_RE74t 一致
Data_RE74c = Data2.loc[Data2['TREAT'] == 0] # 實驗數據，control
# 各組數據的統計性質
def Cmpt_mean_se(df, v):
mean = np.mean( df[v] )
stdErr = np.std( df[v], ddof=1 ) # 自由度 ddof=1 樣本標準差，除以 n-1，無偏；ddof=0 總體標準差，除以 n，有偏
return ["{:.3f}".format(mean), "{:.3f}".format(stdErr)]
idx = pd.MultiIndex.from_product([
['RE74 subset: Treated', 'RE74 subset: Control', 'Comparison group: PSID-1'],
['mean', 'std.err']
])
clmn = ['Age', 'Education', 'Black', 'Hispanic', 'No degree', 'Married', 'RE74 (US$)', 'RE75 (US$)']
table1 = pd.DataFrame(index = idx, columns = clmn + ['No. of observation'])
vLst = ['AGE', 'EDUC', 'BLACK', 'HISP', 'NODEGREE', 'MARR', 'RE74', 'RE75']
for j in range(8):
lst = Cmpt_mean_se(Data_RE74t, vLst[j]) + Cmpt_mean_se(Data_RE74c, vLst[j]) + Cmpt_mean_se(Data_PSID1, vLst[j])
table1[clmn[j]] = lst
cnt_RE74t = len(Data_RE74t)
cnt_RE74c = len(Data_RE74c)
cnt_PSID1 = len(Data_PSID1)
table1['No. of observation'] = [cnt_RE74t, '', cnt_RE74c, '', cnt_PSID1, '']
table1隨機試驗的結果是相對可靠的，這裡以之為標準衡量匹配的效果。我們先看看不匹配直接拿實驗組和觀測組做回歸（具體公式見代碼或原文），結果全部慘不忍睹：
# 復現表 2，panal B C
Data2['AGESQ'] = Data2.apply(lambda x: (x['AGE'] ** 2), axis=1)
idx = pd.MultiIndex.from_product([
['NSW (experimental)', 'PSID-1 (observational)'],
['mean', 'std.err']
])
table2 = pd.DataFrame(index=idx, columns=[0,1,2,3,4,5])
# regression
formula_un_OLS = 'RE78 ~ TREAT'
formula_a_OLS = 'RE78 ~ TREAT + AGE + AGESQ + EDUC + BLACK + HISP + NODEGREE' # panal B: not use RE74
formula_a2_OLS = 'RE78 ~ TREAT + AGE + AGESQ + EDUC + BLACK + HISP + NODEGREE + RE74' # panal C: use RE74
# experimental
result_un_ep_OLS = smf.ols(formula_un_OLS, Data2).fit() # panal B & C
result_a_ep_OLS = smf.ols(formula_a_OLS, Data2).fit() # panal B
result_a2_ep_OLS = smf.ols(formula_a2_OLS, Data2).fit() # panal C
# observational
result_un_ob_OLS = smf.ols(formula_un_OLS, Data1).fit() # panal B & C
result_a_ob_OLS = smf.ols(formula_a_OLS, Data1).fit() # panal B
result_a2_ob_OLS = smf.ols(formula_a2_OLS, Data1).fit() # panal C
table2[0] = [
"{:.0f}".format(result_un_ep_OLS.params['TREAT']),
"{:.0f}".format(result_un_ep_OLS.bse['TREAT']),
"{:.0f}".format(result_un_ob_OLS.params['TREAT']),
"{:.0f}".format(result_un_ob_OLS.bse['TREAT'])
]
table2[1] = [
"{:.0f}".format(result_a_ep_OLS.params['TREAT']),
"{:.0f}".format(result_a_ep_OLS.bse['TREAT']),
"{:.0f}".format(result_a_ob_OLS.params['TREAT']),
"{:.0f}".format(result_a_ob_OLS.bse['TREAT'])
]
table2[2] = [
"{:.0f}".format(result_a2_ep_OLS.params['TREAT']),
"{:.0f}".format(result_a2_ep_OLS.bse['TREAT']),
"{:.0f}".format(result_a2_ob_OLS.params['TREAT']),
"{:.0f}".format(result_a2_ob_OLS.bse['TREAT'])
]
# quasi-difference in earnings between 1978 and 1975
# controlling RE75
formula_un_q = 'RE78 ~ TREAT + RE75'
formula_a_q = 'RE78 ~ TREAT + RE75 + AGE + AGESQ + EDUC + BLACK + HISP + NODEGREE'
formula_a2_q = 'RE78 ~ TREAT + RE75 + AGE + AGESQ + EDUC + BLACK + HISP + NODEGREE + RE74'
result_un_ep_q = smf.ols(formula_un_q, Data2).fit()
result_a_ep_q = smf.ols(formula_a_q, Data2).fit()
result_a2_ep_q = smf.ols(formula_a2_q, Data2).fit()
result_un_ob_q = smf.ols(formula_un_q, Data1).fit()
result_a_ob_q = smf.ols(formula_a_q, Data1).fit()
result_a2_ob_q = smf.ols(formula_a2_q, Data1).fit()
table2[3] = [
"{:.0f}".format(result_un_ep_q.params['TREAT']),
"{:.0f}".format(result_un_ep_q.bse['TREAT']),
"{:.0f}".format(result_un_ob_q.params['TREAT']),
"{:.0f}".format(result_un_ob_q.bse['TREAT'])
]
table2[4] = [
"{:.0f}".format(result_a_ep_q.params['TREAT']),
"{:.0f}".format(result_a_ep_q.bse['TREAT']),
"{:.0f}".format(result_a_ob_q.params['TREAT']),
"{:.0f}".format(result_a_ob_q.bse['TREAT'])
]
table2[5] = [
"{:.0f}".format(result_a2_ep_q.params['TREAT']),
"{:.0f}".format(result_a2_ep_q.bse['TREAT']),
"{:.0f}".format(result_a2_ob_q.params['TREAT']),
"{:.0f}".format(result_a2_ob_q.bse['TREAT'])
]
table2.columns = pd.MultiIndex.from_product([
['Regression', 'Quasi-difference'],
['Unadjusted', 'Adjusted (not use RE74)', 'Adjusted (use RE74)']
])
table2下面我們為實驗組中的參與者匹配觀測組中的未參與者。
第一步，計算傾向得分。通常使用 Logit 或 Probit 模型：
# logit
formula_logit = 'TREAT ~ AGE + AGESQ + EDUC + EDUCSQ + BLACK + HISP + MARR + NODEGREE + RE74 + RE75 + RE74SQ + RE75SQ + U74BLACK' # table 3 remark e
result_logit = sm.Logit.from_formula(formula_logit, data = Data1).fit()
X = Data1[['AGE', 'AGESQ', 'EDUC', 'EDUCSQ', 'BLACK', 'HISP', 'MARR', 'NODEGREE', 'RE74', 'RE75', 'RE74SQ', 'RE75SQ', 'U74BLACK']]
Data1['ps_logit'] = result_logit.predict(X)
# probit
Y = Data1['TREAT']
result_probit = sm.Probit(Y, sm.add_constant(X)).fit()
Data1['ps_probit'] = result_probit.predict(sm.add_constant(X))
Data1[['ps_logit', 'ps_probit']]通過箱線圖，我們可以看到實驗組和觀測組之間是非常不平衡的：
# plot boxplot
plt.figure(figsize=(10, 6))
sns.boxplot(x='TREAT', y='ps_logit', data=Data1)
plt.ylabel('Propensity Score')
plt.xlabel('Observed and Treated')第二步，檢查共同支撐（common support）。
我們丟棄觀測組中傾向得分低於實驗組傾向得分最小值、以及高於實驗組傾向得分最大值的數據，保證只比較有相似的傾向得分的實驗組和觀測組對象，然後分析數據的分布。
# Check the assumptions: common support
Data_t = Data1.loc[Data1['TREAT'] == 1]
max_t_l = np.max(Data_t['ps_logit'])
min_t_l = np.min(Data_t['ps_logit'])
Data_c = Data1.loc[Data1['TREAT'] == 0]
Data_c_ld = Data_c.loc[ (Data_c['ps_logit'] > min_t_l) & (Data_c['ps_logit'] < max_t_l) ]
# 觀察數據
print('treated min: ' + str(min_t_l))
print('observational min: ' + str(np.min(Data_c['ps_logit'])))
print('treated max: ' + str(max_t_l))
print('observational max: ' + str(np.max(Data_c['ps_logit'])))
print('----')
print('Discard ' + str( len(Data_c) - len(Data_c_ld) ) + ' PSID-1 units')
# [0, 1]，gap 0.05 --> stratum
# count PSID units in each stratum
gap = .05
N = int( 1 / gap)
cnt_t = np.zeros(N) # no. of treated units in each stratum
cnt_PSID = np.zeros(N) # no. of PSID units in each stratum
for i in Data_t.index:
n = int( np.floor( Data_t.loc[i, 'ps_logit'] / gap ) )
cnt_t[n] += 1
for i in Data_c_ld.index:
n = int( np.floor( Data_c_ld.loc[i, 'ps_logit'] / gap ) )
cnt_PSID[n] += 1
df_cnt = pd.DataFrame(
index=['0-.05', '.05-.1', '.1-.15', '.15-.2', '.2-.25', '.25-.3', '.3-.35', '.35-.4', '.4-.45', '.45-.5', '.5-.55', '.55-.6', '.6-.65', '.65-.7', '.7-.75', '.75-.8', '.8-.85', '.85-.9', '.9-.95', '.95-1'],
columns=['No. of treated units', 'No. of PSID-1 units']
)
df_cnt['No. of treated units'] = cnt_t
df_cnt['No. of PSID-1 units'] = cnt_PSID
print(df_cnt.T)
x = [.025 + .05 * i for i in range(N)]
plt.figure(dpi=300)
plt.bar(x = x, height = cnt_t, color = 'yellow', label="treated", alpha=0.5, width=.045)
plt.bar(x = x, height = cnt_PSID, color = 'steelblue', label="PSID-1", alpha=0.5, width=.045)
plt.ylim(0, 150)
plt.xlabel('Propensity Score')
plt.ylabel('Frequency')
plt.legend(loc=2, bbox_to_anchor=(1.01,1.0),borderaxespad = 0.)
plt.show()這裡丟棄了超過一半的觀測組數據，但是數據的分布仍然很不均勻，這種情況下 PSM 的效果可能不會好。另外，如果是在丟棄實驗組數據的情況中，如果過多實驗組數據被丟棄，結果可能有偏。
第三步，匹配。
匹配的方式有很多，這裡介紹最常見的四種：
1. 區間匹配：設定一系列區間，傾向得分在同一個區間內的參與者和未參與者為一個配對。上面的傾向得分頻數圖就是按照從 0 到 1，以 0.05 為區間長度的區間匹配。
2. 近鄰匹配：確定要選取的近鄰數量 n，對每個參與者，選取傾向得分離他最近的 n 個未參與者為一個配對。（可以要求每個未參與者只能被匹配一次，比如最近鄰匹配中，一個未參與者同時是兩個參與者的最近鄰，我們將它與更靠近的參與者配對，另一個參與者只能選擇第二靠近的。這樣可以避免未參與者過於稀疏，大量參與者都與同一個未參與者配對，從而造成偏誤的情況。）
# nearest neighbors matching
# KNN -> top5 nearest neighbors
def find_matches(treated, observed, topn=5):
neigh = NearestNeighbors(n_neighbors=topn, metric='euclidean')
neigh.fit(observed) # be matched
distance, location = neigh.kneighbors(treated) # location is not the index in df # can add "return_distance=False"
return location[:, :topn]
treated = np.matrix(Data_t.ps_logit).reshape(-1,1)
observed = np.matrix(Data_c_ld.ps_logit).reshape(-1,1)
location = find_matches(treated, observed)部分結果如下：
3. 半徑匹配：以每個參與者的傾向得分為圓心，設定一個半徑，每個圓內的參與者和未參與者為一個配對。這種方法類似於區間匹配。
# Radius Matching
def find_matches_Radius(Radius):
df_Radius = pd.DataFrame(columns = Data_t.index, index = ['matches', 'no. of matches'])
for i in Data_t.index:
lst = []
n = 0
for j in Data_c_ld.index:
if np.abs(Data_t.loc[i,'ps_logit'] - Data_c_ld.loc[j,'ps_logit']) < Radius:
lst.append(j)
n += 1
df_Radius[i] = [lst, n]
Data_m = pd.concat([Data_t, df_Radius.T], axis=1)
return Data_m.loc[Data_m['no. of matches'] != 0]部分結果如下：
4. 核匹配：選定一個核函數（一般選取像高斯核這樣的鐘形曲線），對第 i 個參與者，按下式計算每個未參與者的權重：
# kernel matching
# weight_ij = K[(p_i - p_j)/h] / sum_j{ K[(p_i - p_j)/h] }
# for treated_i and observed_j
# h: bandwidth
# kernel：Gaussian and Epanechnikov, 其他：biweight, triweight，cosinus
def Gaussian(a):
return math.exp(- a**2 / 2) / (2*math.pi)**.5
def Epanechnikov(a):
if np.abs(a) <= 1:
return .75 * (1 - a**2)
else:
return 0
h_G = .07
h_E = .2
df_weights_G = pd.DataFrame(index = Data_c_ld.index, columns = Data_t.index)
df_weights_E = pd.DataFrame(index = Data_c_ld.index, columns = Data_t.index)
for j in Data_c_ld.index:
for i in Data_t.index:
df_weights_G.loc[j,i] = Gaussian( (Data_t.loc[i,'ps_logit'] - Data_c_ld.loc[j,'ps_logit']) / h_G )
df_weights_E.loc[j,i] = Epanechnikov( (Data_t.loc[i,'ps_logit'] - Data_c_ld.loc[j,'ps_logit']) / h_E )
for i in Data_t.index:
df_weights_G[i] = df_weights_G[i]/np.sum(df_weights_G[i])
df_weights_E[i] = df_weights_E[i]/np.sum(df_weights_E[i])權重矩陣結果：
第四步，平衡性檢驗（下面僅以最近鄰匹配為例）。
匹配的目的是讓數據平衡，所以每次匹配完都要先檢查數據是否平衡。定性的，匹配前後實驗組觀測組年齡的頻數直方圖分別是：
可以看到明顯的改善。
定量的，我們要進行平衡性檢驗，有兩種方法：
1. 檢驗樣本的均值是否沒有顯著差異，所以先進行異方差 Levene-test，然後對不同結果進行同/異方差 t-test ：
test_table = pd.DataFrame(index = ['before matching', '1:1 matching', '1:2 matching', '1:4 matching'])
alpha = .01
# before matching
diffb = np.abs(np.mean(Data_c.AGE) - np.mean(Data_t.AGE))
lpb = CmptLevene(Data_t, Data_c, 'AGE') # 選用 Levene-test，不要求正態，結果更穩健
if lpb < alpha: # p值小，拒絕原假設，即方差不齊
tpb = CmptT_diffVar(Data_t, Data_c, 'AGE')
instruction0 = "Welch's unequal variances t-test"
else:
tpb = CmptT_sameVar(Data_t, Data_c, 'AGE')
instruction0 = 't-test'
if tpb < alpha:
conclusion0 = 'different'
else:
conclusion0 = 'not different'
# -
# 1:1 matching
diff1 = np.abs(np.mean(Data_t.AGE) - np.mean(m1Lst))
lp1 = scipy.stats.levene(Data_t.AGE, m1Lst).pvalue
if lp1 < alpha:
tp1 = scipy.stats.ttest_ind(Data_t.AGE, m1Lst, equal_var=False).pvalue
instruction1 = "Welch's unequal variances t-test"
else:
tp1 = scipy.stats.ttest_ind(Data_t.AGE, m1Lst).pvalue
instruction1 = 't-test'
if tp1 < alpha:
conclusion1 = 'different'
else:
conclusion1 = 'not different'
# -
# 1:2 matching
diff2 = np.abs(np.mean(Data_t.AGE) - np.mean(m2Lst))
lp2 = scipy.stats.levene(Data_t.AGE, m2Lst).pvalue
if lp2 < alpha:
tp2 = scipy.stats.ttest_ind(Data_t.AGE, m2Lst, equal_var=False).pvalue
instruction2 = "Welch's unequal variances t-test"
else:
tp2 = scipy.stats.ttest_ind(Data_t.AGE, m2Lst).pvalue
instruction2 = 't-test'
if tp2 < alpha:
conclusion2 = 'different'
else:
conclusion2 = 'not different'
# -
# 1:4 matching
diff4 = np.abs(np.mean(Data_t.AGE) - np.mean(m4Lst))
lp4 = scipy.stats.levene(Data_t.AGE, m4Lst).pvalue
if lp4 < alpha:
tp4 = scipy.stats.ttest_ind(Data_t.AGE, m4Lst, equal_var=False).pvalue
instruction4 = "Welch's unequal variances t-test"
else:
tp4 = scipy.stats.ttest_ind(Data_t.AGE, m4Lst).pvalue
instruction4 = 't-test'
if tp4 < alpha:
conclusion4 = 'different'
else:
conclusion4 = 'not different'
# -
test_table['diff'] = [diffb, diff1, diff2, diff4]
test_table['Levene-test p-value'] = [lpb, lp1, lp2, lp4]
test_table['instruction'] = [instruction0, instruction1, instruction2, instruction4]
test_table['t-tset p-value'] = [tpb, tp1, tp2, tp4]
test_table['conclusion (alpha = 1%)'] = [conclusion0, conclusion1, conclusion2, conclusion4]
test_table2. 計算均值絕對標準化差異（absolute standardized differences in means）：
其中，分子為實驗組和配對後觀測組均值之差，分母為實驗組和配對前觀測組方差均值的開方。一般我們認為，當這個值小於 0.25 且實驗組和控制組方差的比值在 0.5 和 2 之間時，我們認為沒有明顯差異。（個人認為這是一種簡化版的 t-test，其中方差比值在 0.5 和 2 之間這個條件也是 t-test 中粗略判定同方差的條件。）
# 1:1 matching
df = Data_match['knn_1']
m1_list = []
for i in df:
m1_list.append(Data_c.loc[i,"AGE"])
mean_control_after = np.mean(m1_list)
mean_treat_after = np.mean(Data_t.AGE) # take AGE as an example
var_treat_before = np.var(Data_t.AGE)
var_control_before = np.var(Data_c_ld.AGE)
asmd_age = np.abs(mean_treat_after - mean_control_after) / ( (var_treat_before + var_control_before)/2 ) ** .5
print(asmd_age)
var_ratio = var_treat_before/var_control_before
print(var_ratio)第五步，計算實驗組平均處理效應和對應方差。
分別計算每個區間內的實驗組平均處理效應，最後按照每個區間內實驗組個體的個數加權平均，得到總的處理效應。這裡將各種方法的所有結果整理展示如下：
可以看到，在匹配後結論或多或少有所改善，但是都沒有通過顯著性檢驗，因為匹配丟棄數據會增大標準差，降低結論的可信度，比如這裡有些匹配方法只用了約百個觀測組數據，標準差非常大。有些結論，比如最近鄰匹配和半徑為 0.002 的匹配效果不好的原因，可能是因為在傾向得分大於 0.8 的部分觀測數據太少，大量參與者無法找到匹配或者匹配並不接近。雖然如此，但是我們還是可以利用各種匹配方法，得到一個平均因果效應的估計範圍，下一些定性的結論。
總結總結一下，PSM 的基本步驟是：
PSM 的優點是通過傾向得分，平衡了控制組和幹預組的特徵，提供了較為可信的反事實結果。但是，越複雜的模型假設越強，設定條件越多，也就越脆弱，psm 是有很強的局限性的。這裡引用一些他人的評價：
陳強教授在《高級計量經濟學及stata應用》中指出：
PSM 通常要求比較大的樣本容量以達到高質量的匹配。
PSM 要求除立足於控制組的傾向得分有較大的共同取值範圍；否則會丟失較多的觀測值，導致剩下的樣本不具有代表性（比如這裡 Lalonde 的例子）。
PSM 只控制了可側變量的影響，如果仍然存在依不可測變量選擇，仍然會帶來隱形偏差。
Gary King 等人的批評更加激烈，他們在 2015 年的工作論文：Why propensity scores should NOT be used for matching 中指出，雖然 PSM 的初衷是減少實驗組和控制組在效應之前的協變量不平衡性，進而削減模型依賴性及偏誤，但是實際上 PSM 經常適得其反。文章在一些原本就平衡，以及專門為 PSM 設計的數據集上應用此方法，結果匹配後協變量反而不平衡了。
參考文獻[1] 陳強(2014)：高級計量經濟學及stata應用（第二版）。高等教育出版社。
[2] Robert J Lalonde(1986): Evaluating the econometric evaluation of training programs with Experimental data. America Econamic Review, 76, 406-620.
[3] Rajeev H. Dehejia and Sadek Wahba(1999): Causal effects in nonexperimental studies: reevaluating the evaluation of training programs. Journal of the American Statistical Association, 94(448), 1053 - 1062.
[4] Gary King and Richard Nielsen(2019): Why propensity scores should not be used for matching. Political Analysis, 27(4), 435-454.