统计_结构因果模型：模型驱动大于数据驱动

全文均基于 《为什么：关于因果关系的新科学》

一、决定员工工资的更重要因素探索

在我们的生活当中会看到一些学历不高，由于经验非常丰富于是在一家公司里其工资很高；也会看到一些学历很高经验不足于是收入会不及一些低学历同事的情况。所以我们想要了解下到底是工作经验重要，还是学历更加重要。

为了便于研究和声明，我们用EX表示工作经验，ED表示学历，S表示工资。同时，假设只有三种学历水平：0=高中学历，1=大学学历，2=研究生学历。

雇员	EX(u) /年	ED(u)	$S_{ED=0}(u)$	$S_{ED=1}(u)$	$S_{ED=2}(u)$
爱丽丝	6	0	81000	？	？
博特	9	1	？	92500	？
卡罗琳	9	2	？	？	97500
戴维	8	1	？	91000	？
厄内斯特	12	1	？	100000	？
弗拉西斯	13	0	97000	？	？

1.1 简单线性回归（数据驱动）

一般简单粗暴的思路是，直接用数据做线性回归模型

from sklearn.linear_model import LinearRegression
import pandas as pd
df = pd.DataFrame({
    'name':['爱丽丝', '博特', '卡罗琳', '戴维', '厄内斯特', '弗拉西斯'],
    'ex':[6, 9, 9, 8, 12, 13],
    'ed': [0,1,2,1,1,0],
    's':[81000, 92500, 97000, 91000, 100000, 97000]
})

lr = LinearRegression()
lr.fit(df[['ex', 'ed']].values, df['s'].values)

lr_expr = f' S = {lr.intercept_:.0f} + EX*{lr.coef_[0]:.0f} + ED*{lr.coef_[1]:.0f}'
print(lr_expr)
"""
S = 67331 + EX*2301 + ED*4671
"""

当我们要去预测一个新员工（暂时命名为X）的收入时，可以依据其工作经验(假设6)，其学历(假设1)
。$S_x=67331+6\ast 2301+1\ast 4671=85808$
显然，新员工X的工作经验和爱丽丝是一样的，然后学历比爱丽丝高。那么，假如爱丽丝的学历为1，那么其收入是否就是刚刚计算出的85808呢？我们需要注意的是，对于同一个个体而言，其学历升高，在同等情况下其工作经验必定是会减少的。

我们需要从绘制一张因果图。我们的因果图会对数据背后的因果叙述进行编码，根据因果图，可以看出，工作经验“听从于”学历，并且工资“听从于“工作经验和学历两者。

1.2 结构因果模型

（1）工资“听从于“工作经验和学历两者

$$S=67331+EX\ast 2301+ED\ast 4671+U_s......(1)$$
考虑到个体差异，需要加入一个变量$U_s$。
$U_s$代表影响工资的其他未观测到的变量

（2）工作经验“听从于”学历

需要构建工作经验与学历的模型。这个时候直接从数据建立模型会存在很大的偏差，会发现学历和工作经验的模型和现实情况存在很大偏差(学历上升1，其工作经验减少的时间应该是大于1的)

lr_ex = LinearRegression()
lr_ex.fit(df[['ed']].values, df['ex'].values)
lr_ex_expr = f' EX = {lr_ex.intercept_:.2f} + ED*{lr_ex.coef_[0]:.2f} + U_EX'
print(lr_ex_expr)
"""
EX = 9.65 + ED*-0.18 + U_EX
"""

我们可以基于数据和显示情况调整模型，从仅有的两条ED=0的数据我们可以看出，这部分的数据中高中学历的员工平均工作时间是10年左右(两个数据的均值为9.5；以上模型的截距为9.65)。从高中学历到本科需要4年，从大学到研究生有3年，2年，1年这三种情况，为了简化模型，假设2个阶段的年限均为4年。这个时候我们可以给出以下工作经验关于学历的模型：
$$EX=10-4\ast ED+U_{EX}......(2)$$
(1)和(2)两个方程便是：结构因果模型

（3）结构因果模型反事实预测——假设爱丽丝有大学学历其现在工资是多少

1. (外展)利用关于爱丽丝和其他员工的数据来估计爱丽丝的特质因子：$U_S(爱丽丝)$，$U_{EX}(爱丽丝)$
2. (干预)利用do算子改变模型，以反映我们提出的反事实假设，在案例中即，假如爱丽丝有大学学位：$ED(爱丽丝)=1$
3. (预测)利用修改后的模型及有关外生变量的更新信息$U_S(爱丽丝)$，$U_{EX}(爱丽丝)$ 和 $ED(爱丽丝)$来估算爱丽丝的工资水平。新的工资水平就等于$S_{ED=1}(爱丽丝)$

步骤一：

def compute_Us(original_ex, original_ed, original_s):
    """
    计算特质因子U_s
    """
    return original_s - 67331 - original_ex*2301 - original_ed*4671

def compute_Uex(original_ex, original_ed):
    """
    计算特质因子U_ex
    """
    return original_ex - 10 + 4*original_ed

def recompute_s(original_ex, original_ed, original_s, assume_ed, name=None):
    """
    基于特质U_s, U_ex，假设学历为assume_ed的情况下，重新评估员工的薪资情况
    """
    us = compute_Us(original_ex, original_ed, original_s)
    uex = compute_Uex(original_ex, original_ed)
    print(f'{name}的特质因子：U_s={us}， U_ex={uex}')
    ex_ = 10 - 4 * assume_ed + uex
    return 67331 + ex_*2301 + assume_ed*4671 + us

original_ex = 6
original_ed = 0
original_s = 81000

recompute_s(original_ex, original_ed, original_s, assume_ed=1, name='爱丽丝')
"""
爱丽丝的特质因子：U_s=-137， U_ex=-4
76467
"""

二、将表格中的问号填充完整&同条件比较

2.1 将表格中的问号填充完整

from collections import namedtuple
assume_result = namedtuple('AssumeResult','name assume_ed s')

res_list = []
ed_set = {0, 1, 2}
gp = df.groupby('name')
for name, tmp_df in gp:
    print(name)
    original_ex, original_ed, original_s = tmp_df.values[0][1:]
    assume_ed_set = ed_set - {original_ed}
    for assume_ed in assume_ed_set:
        s = recompute_s(original_ex, original_ed, original_s, assume_ed=assume_ed, name=name)
        res_list.append(assume_result(name, assume_ed, s))

print(res_list)   
print(pd.DataFrame(res_list))
"""
[AssumeResult(name='博特', assume_ed=0, s=97033),
 AssumeResult(name='博特', assume_ed=2, s=87967),
 AssumeResult(name='卡罗琳', assume_ed=0, s=106066), 
 AssumeResult(name='卡罗琳', assume_ed=1, s=101533), 
 AssumeResult(name='厄内斯特', assume_ed=0, s=104533), 
 AssumeResult(name='厄内斯特', assume_ed=2, s=95467), 
 AssumeResult(name='弗拉西斯', assume_ed=1, s=92467), 
 AssumeResult(name='弗拉西斯', assume_ed=2, s=87934), 
 AssumeResult(name='戴维', assume_ed=0, s=95533), 
 AssumeResult(name='戴维', assume_ed=2, s=86467), 
 AssumeResult(name='爱丽丝', assume_ed=1, s=76467), 
 AssumeResult(name='爱丽丝', assume_ed=2, s=71934)]


    name  assume_ed       s
0     博特          0   97033
1     博特          2   87967
2    卡罗琳          0  106066
3    卡罗琳          1  101533
4   厄内斯特          0  104533
5   厄内斯特          2   95467
6   弗拉西斯          1   92467
7   弗拉西斯          2   87934
8     戴维          0   95533
9     戴维          2   86467
10   爱丽丝          1   76467
11   爱丽丝          2   71934
"""

雇员	EX(u) /年	ED(u)	$S_{ED=0}(u)$	$S_{ED=1}(u)$	$S_{ED=2}(u)$
爱丽丝	6	0	81000	76467	71934
博特	9	1	97033	92500	87967
卡罗琳	9	2	106066	101533	97500
戴维	8	1	95533	91000	86467
厄内斯特	12	1	104533	100000	95467
弗拉西斯	13	0	97000	92467	87934

本数据仅仅是作为反事实预测示意参考，笔者并不认同读书无用之说。提前工作和继续读书在社会上存在一些异议，这里不做其他讨论。

2.2 同条件比较

如果我们需要比较 爱丽丝 和 博特 在同一工作经验和学历的情况下的工资水平，那么我们需要对do算子进行修改，增加一个工作经验的修正：在特质因子$U_{EX}(爱丽丝)$计算出的$E{X}_{new}$的基准上增加上工作年限$E{X}_{extra}$以达到假设工作经验

def recompute_s_assume_edex(original_ex, original_ed, original_s, assume_ed, assume_ex, name=None):
    """
    基于特质U_s, U_ex，
    假设学历为assume_ed, 假设工作经验为assume_ex的情况下，
    重新评估员工的薪资情况
    """
    us = compute_Us(original_ex, original_ed, original_s)
    uex = compute_Uex(original_ex, original_ed)
    print(f'{name}的特质因子：U_s={us}， U_ex={uex}')
    ex_new = 10 - 4 * assume_ed + uex
    ex_extra = assume_ex - ex_new
    return 67331 + ex_new*2301 + assume_ed*4671 + us + ex_extra*2301

recompute_s_assume_edex(original_ex, original_ed, original_s, assume_ed=1, assume_ex=9, name='爱丽丝')

"""
爱丽丝的特质因子：U_s=-137， U_ex=-4
92574
"""

可以看出$S_{ED=1,EX=9}(爱丽丝)=92574$， $S_{ED=1,EX=9}(博特)=92500$。如果爱丽丝有大学学历且工作9年，其收入会高于博特，即爱丽丝的工作能力比博特强。

要做反事实的比较，我们需要构建结构因果模型而不是简单的数据拟合。