python虽然与R一样都可以做数据分析，但是在计量方面较为薄弱，python更像是干脏活，清洗数据用的。现在慢慢的python也有一些在计量的包，比如causalinference，这个包可以做因果推断分析。
安装
!pip3 install causalinference
Lookinginindexes:https://pypi.tuna.tsinghua.edu.cn/simple
Collectingcausalinference
Downloadinghttps://pypi.tuna.tsinghua.edu.cn/packages/dc/7f/4504b42ef5a1158075954f54d08b95b2d5b2186da0ef9fcbcd0cf31411f2/CausalInference-0.1.3-py3-none-any.whl (51 kB)
[K|████████████████████████████████|51kB81kB/s eta0:00:0101
[?25hInstallingcollected packages:causalinference
Successfullyinstalled causalinference-0.1.3
数据导入
importpandasaspd
df=pd.read_csv('data.csv')
df
数据描述
x1，x2，x3 协变量(控制变量)
y 因变量
istreatment 处置变量D，标注每条数据隶属于treatment或control组。1为treatment， 0为control。
fromcausalinferenceimportCausalModel
Y=df['y'].values
D=df['istreatment'].values
X=df[['x1','x2','x3']].values
#CausalModel参数依次为Y， D， X。其中Y为因变量
causal=CausalModel(Y,D,X)
causal
描述性统计分析
print(causal.summary_stats)
SummaryStatistics
Controls(N_c=2509)Treated(N_t=2491)
VariableMeanS.d.MeanS.d.Raw-diff
--------------------------------------------------------------------------------
Y-1.0121.7424.9783.0685.989
Controls(N_c=2509)Treated(N_t=2491)
VariableMeanS.d.MeanS.d.Nor-diff
--------------------------------------------------------------------------------
X0-0.3430.9400.3360.9610.714
X1-0.3470.9360.3450.9580.730
X2-0.3130.9400.3060.9630.650
causal.summary_stats含有的指标字段名
causal.summary_stats.keys()
dict_keys(['N','K','N_c','N_t','Y_c_mean','Y_t_mean','Y_c_sd','Y_t_sd','rdiff','X_c_mean','X_t_mean','X_c_sd','X_t_sd','ndiff'])
使用OLS估计处置效应
估计处置效应最简单的方法是使用OLS方法，
CausalModel.est_via_ols(adj)
该方法有一个参数adj
adj=0 模型未使用X(协变量)
adj=1 模型使用了D(是否为处置组)和X(协变量)。
adj=2 模型使用了D(是否为处置组)、X(协变量)、D与X的交互
adj默认为2
causal.est_via_ols(adj=2)
print(causal.estimates)
TreatmentEffectEstimates:OLS
Est.S.e.z P>|z|[95%Conf.int.]
--------------------------------------------------------------------------------
ATE3.0170.03488.7400.0002.9503.083
ATC2.0310.04051.1830.0001.9532.108
ATT4.0100.039103.9640.0003.9344.086
参数解读
ATE 平均处置效应(average treatment effect)
ATC 控制组的平均处置效应(average treatment effect for the controls)
ATT 处置组的平均处置效应(average treatment effect for the treated)
你们再试试adj设置为0和1分别运行出什么结果
倾向得分估计
我们估计处置效应时，很希望处置组和控制组很类似。比如研究受教育水平对个人收入的影响，其他变量如家庭背景、年龄、地区等协变量存在差异，我们希望控制组和处置组的之间的协变量平衡性尽可能的好，这样两个组就会很像，当对这两个组的受教育水平进行操作时，两个组的收入差异可以认为是受教育水平带来的。
让两个组很像，这里就用到倾向得分估计。
causal.est_propensity_s()
print(causal.propensity)
EstimatedParametersofPropensityScore
Coef.S.e.z P>|z|[95%Conf.int.]
--------------------------------------------------------------------------------
Intercept0.0050.0350.1450.885-0.0630.073
X10.9990.04124.4950.0000.9191.079
X01.0000.04124.5430.0000.9201.080
X20.9330.04023.1810.0000.8551.012
分层方法估计处置效应
倾向得分估计，让两个组尽量相似，但实际上这个相似值范围有点大。比如假设受教育水平对个人收入的影响，身高、体重等颜值信息(协变量)其实对收入也是有影响的，那么就应该对人群进行分层，不同颜值水平(分组)下受教育水平对个人收入的影响。
分层方法估计CausalModel.stratify_s() 自动选择协变量
causal.stratify_s()
print(causal.strata)
StratificationSummary
PropensityScoreSampleSizeAve.PropensityOutcome
StratumMin.Max.ControlsTreatedControlsTreatedRaw-diff
--------------------------------------------------------------------------------
10.0010.04315350.0240.029-0.049
20.0430.06914880.0560.0590.142
30.0700.118283290.0930.0920.953
40.1190.178268450.1470.1471.154
50.1780.240247650.2080.2101.728
60.2400.3614511740.2990.3002.093
70.3610.4271961170.3930.3952.406
80.4270.4991531590.4650.4642.868
90.4990.53282750.5150.5152.973
100.5320.56865910.5510.5533.259
110.5680.6301141980.6000.6013.456
120.6300.7581804450.6930.6963.918
130.7580.818772360.7870.7894.503
140.8180.876572550.8450.8494.937
150.8760.933232890.9040.9045.171
160.9330.998123000.9570.9636.822
更多详细信息可阅读代码中说明论文，可在项目中下载到的。
