分析师入门常见错误 "幸存者偏差"，如何用匹配和加权法规避
在日常功能迭代分析中，一般会直接看使用该功能和未使用该功能的用户在成功指标上的表现，将两组数据求个差异值就得出功能的效果结论。但是有敏锐的分析师会发现，功能大部分情况下有筛选效应，即使用该功能的用户可能本身质量比较高，活跃比较频繁。用以上的方法估计会导致效果评估失真，那么如何规避混杂因素导致的幸存者偏差。优先考虑的做法是探究一些相关关系因素，用 A/B 测试验证，把因果推断作为备选或探索式分析的手段，但有些场景无法进行 A/B 测试。这里介绍因果推断中的两个方法——匹配和逆概率加权。并将其和直接回归方法的结论进行对比，看看相关和因果的结论到底会差异多少。
因果推断方法一般的通俗理解，A/B 测试不行，创造条件近似 A/B 测试。
案例背景
分析师评估某功能是否可以降低用户流失率。这里的数据不是实验性的，工程上谁也无法控制用户去使用新功能。数据集包含以下列：
- 流失率（Churn_rate）：用户流失的可能性。值越高表示流失的概率越高。
- 是否使用该功能（is_using）：二进制变量，用户是否使用了该功能。
- 日均使用时长（avg_used_time）：用户使用产品的日均时长。
- 用户信用评分（health）：用户自我报告的信用评分状况。以0–100的等级进行测量，数值越高表示信用状况越好。
- 最近一次使用时间（recency）：最近一次使用距离现在的天数。
- 活跃天数（active_days）：最近一个月使用产品的天数。
我们仅使用观察数据来估计该功能使用对流失风险的因果关系。因为这不是随机控制实验（RCT），宣称因果关系似乎有些不准确。但是，如果我们可以绘制正确的 DAG 并针对正确的节点进行调整，则可以得到使用功能→流失的因果关系并求出接近的因果效应。因为这是模拟数据，所以我们知道真正因果关系，即设置真正的平均处理效应（ATE）为 -0.1。平均而言，使用该功能可使流失风险降低 0.1 个单位。
加载数据
library(tidyverse)
library(ggdag) # Make DAGs
library(dagitty) # Do DAG logic with R
library(broom) # Convert models to data frames
library(MatchIt) # Match things
library(modelsummary) # Make side-by-side regression tables
set.seed(1234) # Make all random draws reproducible
data <- read_csv("churn_causal_relationship.csv")
我们建立了这个因果模型，以显示功能使用和流失风险之间关系背后的数据生成过程：
dag <- dagify(
Churn_rate ~ is_using + avg_used_time + health + recency + active_days,
is_using ~ avg_used_time + recency + active_days,
active_days ~ avg_used_time,
exposure = "is_using",
outcome = "Churn_rate",
coords = list(x = c(Churn_rate = 7, is_using = 3, avg_used_time = 4, active_days = 5,
recency = 6, health = 8.5),
y = c(Churn_rate = 2, is_using = 2, avg_used_time = 3, active_days = 1,
recency = 3, health = 2)),
labels = c(Churn_rate = "Risk of Churn", is_using = "Is Using", avg_used_time = "average used time",
health = "Health", recency = "Recency",
active_days = "Active Days")
)
ggdag_status(dag, use_labels = "label", text = FALSE) +
guides(fill = FALSE, color = FALSE) + # Disable the legend
theme_dag()
通过上面的代码，我们得到变量间的因果关系图 DAG：
遵循 do-calculus 规则，我们可以找到所有混淆了功能使用和流失风险间关系的节点变量，控制 混淆变量 可阻断 Is Using 到 Risk of Churn 的 backdoor path, 帮助正确分析 功能使用对流失风险的影响。图看起来很复杂，我们可以直接使用 R 中的方法 adjustmentSets 来找出影响功能使用和流失风险间关系的混淆变量，得到活跃天数 active_days、日均使用时长 avg_used_time 和最近一次使用时间 recency 为混淆变量，将这三个变量进行控制以确定功能使用和流失风险间的真正因果效应量。
adjustmentSets(dag)
# { active_days, avg_used_time, recency }
相关关系不是因果关系
为了方便对比，我们可以计算出那些使用/不使用该功能的用户平均流失风险的差异。注意这绝对不是实际的因果关系，因为该模型并未控制任何的混淆变量，这里可以直接分组计算（使用=1，不使用=0）该功能的用户平均流失风险的差异。
data %>%
group_by(is_using) %>%
summarize(number = n(),
avg = mean(Churn_rate))
# # A tibble: 2 x 3
# is_using number avg
# <dbl> <int> <dbl>
#1 0 1071 0.419
#2 1 681 0.256
或者我们可以通过回归来做到这一点：
model_wrong <- lm(Churn_rate ~ is_using, data = data)
tidy(model_wrong)
# A tibble: 2 x 5
# term estimate std.error statistic p.value
#1 (Intercept) 0.419 0.00405 104. 0.
#2 is_using -0.163 0.00649 -25.1 2.25e-119
根据这一平均估计，使用该功能可将流失风险降低 0.163 点。虽然这是我们日常分析工作中经常容易得出的结论，但是，这不是实际两者关系的因果效应，因为混淆变量在这里并未处理。其实对于混淆变量不多的情况，更简单的方法是 standardize，非参数 standardize 方法就是直接分层得出因果效应，然后按照比例最后加权得出全局的平均因果效应 ATE，但是对于多维的混淆变量，通过以下两个方法可以解决分层方法的维度爆炸问题，但最终得到的结果和standardize 方法相等[1]。
匹配 Matching
我们可以使用匹配方法将相似的样本配对，并提出无混淆的假设，即如果我们看到两个观测样本几乎相同，而一个样本使用了一个功能，而一个样本则没有使用，那么控制到是否使用该功能的选择是随机的。
我们从 DAG 得知活跃天数 active_days、日均使用时长 avg_used_time和最近一次使用时间 recency 会同时影响功能使用和流失风险（即混淆了这两者的关系），所以我们将尝试找到具有相同活跃天数，日均使用时长和最近一次使用时间的观测样本，部分样本使用了该功能，而剩下的没有使用该功能。
我们可以使用 MatchIt R 包中的 matchit() 函数根据马氏距离来进行样本匹配。还有许多其他选项可用，有关详细信息，请参见在线文档。使用 replace = TRUE 可以实现重复匹配（即一对多匹配）。
不可重复匹配使得每个控制组只能匹配一次，即使该控制组是多个处理组的最佳匹配，这就使得匹配质量降低和样本变小。相反，重复匹配则可以有效避免这些问题，但是在估计处理效应时，需进行加权和调整标准误，以反映匹配次数的影响。
预处理
所有 681 个使用该功能的用户都与其相似的未使用该功能的用户（其中 431 个）进行匹配。640 人不匹配，将被丢弃。现在根据样本的混淆变量特征数据已经匹配，排除了混淆变量的影响，可以用关键变量进行建模：
matched_data <- matchit(is_using ~ avg_used_time + active_days + recency,
data = data,
method = "nearest",
distance = "mahalanobis",
replace = TRUE)
summary(matched_data)
matched_data_for_real <- match.data(matched_data)
model_matched <- lm(Churn_rate ~ is_using,
data = matched_data_for_real)
tidy(model_matched)
# A tibble: 2 x 5
# term estimate std.error statistic p.value
#1 (Intercept) 0.376 0.00594 63.3 0.
#2 is_using -0.120 0.00759 -15.8 7.48e-51
这里的 -0.120 点要比未处理混淆变量的估计要准确，但这不是真正的因果效应。可能是因为匹配效果不佳，或丢弃了太多数据。实际上，不准确估计的最大原因是数据中存在一些不平衡，即在完成匹配后需要检验匹配结果是否真的实现了平衡两组的混淆变量水平。因为我们设置 replace = TRUE，我们并没有做到 1:1 匹配，未使用该功能的观察样本与一个及以上的使用该功能的观察样本配对。结果，被多次匹配的观测样本在模型中的重要性太大。matchit() 为我们提供了一个名为 weights 的列，该列使我们可以在运行模型时按比例缩小因过度匹配而引起不平衡的观察值。如果您使用 replace=FALSE 并实施 1:1匹配，则整个 weights 列将仅为 1。我们可以将这些权重合并到模型中，以获得更准确的估算值：
model_matched_wts <- lm(Churn_rate ~ is_using,
data = matched_data_for_real,
weights = weights)
tidy(model_matched_wts)
# A tibble: 2 x 5
# term estimate std.error statistic p.value
#1 (Intercept) 0.357 0.00585 61.0 0.
#2 is_using -0.101 0.00748 -13.5 1.24e-38
在平衡处理过度匹配问题后，我们发现因果效应为 -0.101 点。这比到目前为止我们尝试过的任何其他估计要好得多！
逆概率加权 Inverse probability weighting
匹配方法的一个潜在弊端是，通常必须丢弃大量的数据，即不匹配的数据都不会包含在最终估计数据集中。
逆概率加权方法是首先为每个观察样本分配接受处理（这里是使用该功能）的概率，然后按其相反的概率对每个观察值进行加权，即对于实际得到处理的观测样本，预测大概率将没有得到处理（预测大概率不会使用该功能但实际使用了），比预测很可能得到处理的样本会赋予更大的权重。
生成这些逆概率权重需要两步过程： （1）首先生成倾向得分或接受处理的概率； （2）使用公式将倾向得分转换为权重。一旦有了逆概率权重，就可以将它们合并到回归模型中。
步骤1：倾向得分
有多种方法可以生成倾向得分（例如逻辑回归，概率回归，甚至是机器学习技术，例如随机森林和神经网络），但是逻辑回归可能是最常见的方法。
逻辑回归模型中的结果变量必须是二进制的。logistic 回归中的 Y 是概率的对数比，这迫使模型的输出在0-1范围内，由于是否使用该功能变量是二进制结果，这里采用逻辑回归来计算倾向得分：
当我们在生成倾向得分的模型中包含变量时，就像在匹配中所做的那样，我们处理了混淆变量。但是与匹配不同，该方法不会丢弃任何数据！只是使一些观察样本变得更重要，而另一些则变得不那么重要。
首先，我们建立一个基于活跃天数 active_days、日均使用时长 avg_used_time 和最近一次使用时间 recency 预测是否使用功能的模型（因为这些变量是 DAG 中的混杂因素），然后，我们可以使用数据集中的活跃天数、日均使用时长和最近一次使用时间数据，并使用以下模型生成概率：
model_logit <- glm(is_using ~ avg_used_time + active_days + recency,
data = data,
family = binomial(link = "logit"))
using_probabilities <- augment_columns(model_logit,
data,
type.predict = "response") %>% rename(propensity = .fitted)
# Look at the first few rows of a few columns
using_probabilities %>%
select(id, is_using, avg_used_time, active_days, recency, propensity) %>%
head()
# A tibble: 6 x 6
# id is_using avg_used_time active_days recency propensity
# <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#1 1 1 781 21 2 0.309
#2 2 0 974 27 4 0.421
#3 3 0 502 26 3 0.161
#4 4 1 671 21 5 0.351
#5 5 0 728 19 5 0.413
#6 6 0 1050 25 1 0.380
倾向得分为 propensity 列。考虑到他们的活跃天数 active_days、日均使用时长 avg_used_time 和最近一次使用时间 recency，某些人（如第3个人）不太可能使用该功能（只有 16.1％ 的机会）。部分用户由于活跃天数、日均使用时长较高（如第2个人），他们使用该功能的可能性更高（42.1％）。
接下来，我们需要将倾向得分转换为逆概率权重，这使那些奇怪的观察样本在模型中变得更加重要（比如使用功能的可能性很高但没有使用的人，反之亦然），该公式将创建权重：
我们用 mutate() 创建逆概率权重列：
using_ipw <- using_probabilities %>%
mutate(ipw = (is_using / propensity) + ((1 - is_using) / (1 - propensity)))
# Look at the first few rows of a few columns
using_ipw %>%
select(id, is_using, avg_used_time, active_days, recency, propensity, ipw) %>% head()
# A tibble: 6 x 7
# id is_using avg_used_time active_days recency propensity ipw
#1 1 1 781 21 2 0.309 3.24
#2 2 0 974 27 4 0.421 1.73
#3 3 0 502 26 3 0.161 1.19
#4 4 1 671 21 5 0.351 2.85
#5 5 0 728 19 5 0.413 1.70
#6 6 0 1050 25 1 0.380 1.61
对于不使用功能（即is_using=0）的用户，但propensity 值高，即使用概率高，因此根据实际真实情况，对该类样本需要赋予低权重；而对于真实使用功能，但概率低的样本赋予高权重。比如第4个人！他们只有35.1％的机会使用功能，但是他们却真实使用了！因此，它们具有较高的逆概率权重（3.81）。
步骤2：估算
现在我们已经基于混淆因素生成了逆概率权重，我们可以运行一个模型来查找使用功能对流失风险的因果关系。同样，我们不需要在模型中加入混淆因素，因为我们在创建倾向得分和权重时已经使用了它们：
model_ipw <- lm(Churn_rate ~ is_using,
data = using_ipw,
weights = ipw)
tidy(model_ipw)
# A tibble: 2 x 5
# term estimate std.error statistic p.value
#1 (Intercept) 0.396 0.00464 85.3 0.
#2 is_using -0.103 0.00654 -15.8 1.33e-52
在使用逆概率权重之后，我们发现因果效应为 -0.103 。这与 0.1 的真实值相比还差一点。检查逆概率权重的值很重要。有时它们可能变得太大，例如某个用户的日均使用时长低较低，几乎不活跃，但仍然使用了该功能，那么会拥有非常高的 IPW 值，这种情况可能会推高估算值。要解决此问题，我们可以将权重截断到较低的水平。截断阈值没有通用的经验法则，一般使用 10 作为阈值：
using_ipw <- using_ipw %>%
# If the IPW is larger than 10, make it 10, otherwise use the current IPW
mutate(ipw_truncated = ifelse(ipw > 10, 10, ipw))
model_ipw_truncated <- lm(Churn_rate ~ is_using,
data = using_ipw,
weights = ipw_truncated)
tidy(model_ipw_truncated)
# A tibble: 2 x 5
# term estimate std.error statistic p.value
#1 (Intercept) 0.396 0.00460 86.0 0.
#2 is_using -0.104 0.00649 -16.1 1.76e-54
现在，因果效应为 -0.104，该因果关系系数稍低，准确性可能较低，因为真实情况没有任何极端数据会拉高 IPW 估算值。
所有模型的结果
全文我们只是使用观察数据来估计因果关系。没有随机控制实验（ A/B 实验）的因果关系！让我们比较一下刚刚计算出的所有模型的平均处理效应 ATE 值：
modelsummary(list("Naive" = model_wrong,
"Matched" = model_matched, "Matched + weights" = model_matched_wts,
"IPW" = model_ipw, "IPW truncated at 10" = model_ipw_truncated))
哪一个是对的？在这种情况下，由于这是伪造的模拟数据，我在其中建立了 0.1 的因果效应量，因此我们可以看到其中哪个模型最接近：在这里，Matched+weights 模型获胜（-0.101）。但在现实中，我们不会知道真正的值，匹配和 IPW 都可以很好地对混杂因素进行调整。因此可以尝试多种方式得到多个值评估。
后台回复“ 匹配 ”获取数据。
参考资料
[1] 因果推断数据分析实战: https://zhuanlan.zhihu.com/p/171703047/