从0开始学习R语言--Day26--因果推断

很多时候我们在探讨数据的相关性问题时，很容易会忽略到底是数据本身的特点还是真的是因为特征的区分导致的不同，从而误以为是特征起的效果比较大。

这就好比测试一款新药是否真的能治病，假如吃药的患者康复的更快，那到底是因为药物本身的效果好，还是因为患者本身更健康，平时有控制饮食合理作息与运动，从而在患病后更快地凭借自身免疫力战胜病毒。这需要我们意识到对照试验还需要人为地补足某些条件，也就是探讨是否真的是X导致了Y。

以下是一个例子：

# 加载必要包
library(tidyverse)
library(broom)# 生成模拟数据集（1000名患者）
set.seed(123)
data <- tibble(# 年龄影响病情严重程度和治疗选择age = rnorm(1000, mean = 50, sd = 10),# 病情严重程度（年龄越大病情越重）severity = 0.3 * age + rnorm(1000, sd = 5),# 治疗选择（病情越重越可能接受治疗）treatment = rbinom(1000, 1, plogis(-2 + 0.05 * severity)),# 康复时间（治疗有效，但病情越重康复越慢）recovery_time = 30 - 5 * treatment + 0.5 * severity + rnorm(1000, sd = 3)
)# 查看前几行数据
head(data)data1 <- data
data1 %>% group_by(treatment) %>% summarise(mean_recovery = mean(recovery_time))
# 会发现治疗组康复时间更长！（因为治疗组病情更重）model <- lm(recovery_time ~ treatment + severity + age, data = data)
tidy(model) %>% filter(term == "treatment")
# 现在能看到治疗真实效果（约减少5天）library(MatchIt)# 计算倾向得分（基于年龄和病情）
match_model <- matchit(treatment ~ age + severity,data = data,method = "nearest"
)# 匹配后的数据
matched_data <- match.data(match_model)# 分析匹配后的数据
lm(recovery_time ~ treatment, data = matched_data) %>% tidy()ggplot(data, aes(x = severity, y = recovery_time, color = factor(treatment))) +geom_point(alpha = 0.5) +geom_smooth(method = "lm") +labs(title = "康复时间 vs 病情严重程度",subtitle = "控制病情后，治疗组(红色)康复更快",color = "是否治疗")

输出：

  treatment mean_recovery<int>         <dbl>
1         0          37.5
2         1          33.0term      estimate std.error statistic  p.value<chr>        <dbl>     <dbl>     <dbl>    <dbl>
1 treatment    -5.43     0.238     -22.8 3.69e-93term        estimate std.error statistic  p.value<chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    38.3      0.295     130.  0       
2 treatment      -5.36     0.418     -12.8 7.49e-32
`geom_smooth()` using formula = 'y ~ x'

从结果可以看到，单纯的比较康复时间，会隐藏在背后的很多原因，比如患者本身的病情，患者年龄较大，免疫力低等。进一步看线性回归控制变量，控制了相同的病情和年龄去看康复时间，治疗的方案会减少5天多，用倾向得分模拟的随机试验的得到的结果跟控制变量得到的结果差不多，加强了结果的可信度。而从图像来看，治疗组始终低于对照组，也就是说不管怎样，治疗组总是康复得更快，这也暗含了一种因果关系，只是以图像的形式表现。