从0开始学习R语言--Day14--贝叶斯统计与结构方程模型

贝叶斯统计

在很多时候，我们经常会看到在统计分析中出现很多反直觉的结论，比如假如有一种病，人群中的患病率为1%，患者真患病时，检测结果为阳性的概率是99%，如果没有，则检测结果为阳性的概率是5%，但是如果让我们计算一个人检测结果为阳性时真正患病的概率为多少，根据概率公式P(患病|检测为阳性) = [P(检测为阳性|患病)× P(患病)]/P(检测为阳性) 我们可以计算，而P(检测为阳性)可以通过全概率公式P(检测为阳性)=P(检测为阳性|患病)P(患病) + P(检测为阳性|没患病)P(没患病)得出，计算得到的结果为16.67%，可能有人会疑惑，为什么患真病的前提下，检测为阳性时99%，反过来就只有16.67%，其实会有这样的错句，一是因为没抓住基础的概率，即本身该病的患病率就是很低的1%，二是在人群中没病但检测为阳性的人（5%×99%≈4.95%）比有病且检测为阳性的人（1%×99%≈0.99%）多很多。

这种先有对数据样本的初始判断，观察数据在给定假设下的表现，通过后续检验后得出新判断的过程，就是贝叶斯统计，其中，初始的判断叫做先验概率（Prior），即人群中的患病率1%；在某个假设成立的情况下，观察到的表现叫做似然值（Likehood），即患者有病时检测为阳性的概率和没病时检测为阳性的概率，这里是否有病就是假设；在看到计算结果后得出的新判断叫做后验概率（Posterior）。

下面我们依然是给出一个例子来解释：

# 加载必要的库
library(tidyverse)
library(ggplot2)

# 定义参数
p_disease <- 0.01       # 患病先验概率(1%)
p_true_positive <- 0.99 # 真阳性率(99%)
p_false_positive <- 0.05 # 假阳性率(5%)

# 计算检测为阳性的总概率
p_positive <- (p_true_positive * p_disease) + (p_false_positive * (1 - p_disease))

# 计算后验概率(贝叶斯定理)
p_disease_given_positive <- (p_true_positive * p_disease) / p_positive

# 打印结果
cat(sprintf("检测为阳性时实际患病的概率: %.2f%%", p_disease_given_positive * 100))

# 可视化先验、似然和后验
results <- data.frame(
  Scenario = c("Prior", "Likelihood (Disease)", "Likelihood (No Disease)", "Posterior"),
  Probability = c(p_disease, p_true_positive, p_false_positive, p_disease_given_positive)
)

ggplot(results, aes(x = Scenario, y = Probability, fill = Scenario)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = scales::percent(Probability)), vjust = -0.5) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "贝叶斯分析: 疾病检测案例",
       subtitle = "先验、似然和后验概率比较",
       y = "概率") +
  theme_minimal()

# 更详细的模拟分析
set.seed(123)
population_size <- 10000

# 生成模拟人群
sim_data <- data.frame(
  id = 1:population_size,
  disease = rbinom(population_size, 1, p_disease),
  test_result = NA
)

# 模拟检测结果
sim_data <- sim_data %>%
  mutate(
    test_result = case_when(
      disease == 1 ~ rbinom(n(), 1, p_true_positive),
      disease == 0 ~ rbinom(n(), 1, p_false_positive)
    )
  )

# 计算实际阳性中患病的比例
empirical_posterior <- sim_data %>%
  filter(test_result == 1) %>%
  summarise(
    p_disease_given_positive = mean(disease)
  ) %>%
  pull()

cat(sprintf("\n模拟数据中阳性结果实际患病的比例: %.2f%%", empirical_posterior * 100))

# 创建列联表展示结果
table(sim_data$disease, sim_data$test_result) %>%