从0开始学习R语言--Day37--CMH检验

对于有多个特征的数据，我们一般的处理方式是构建特征函数，计算每个特征向量的系数，从而将其影响纳入到研究量中，但对于简单的问题，也这样做的话未免有点小题大做。这时我们可以考虑用CMH来分析变量在每个特征下的影响，这个方法可以通过分层控制不同的无关特征和变量，凸显变量真实的关联关系。

以下是一个例子：

set.seed(123)
n <- 500  # 增大样本量
Age <- sample(c("Young", "Middle", "Old"), n, replace = TRUE, prob = c(0.3, 0.4, 0.3))
Drug <- sample(c("A", "B"), n, replace = TRUE)

# 改进：药物B在Old组更可能有效，但允许例外
Effect <- ifelse(
  (Drug == "B" & Age == "Old" & runif(n) > 0.2) |  # 80% 有效
    (Drug == "A" & Age == "Young" & runif(n) > 0.3) | # 70% 有效
    (Age == "Middle" & Drug == "B" & runif(n) > 0.6) | # Middle组B药40%有效
    (runif(n) > 0.9),  # 10% 的全局随机有效
  "Improved", 
  "Not Improved"
)
df <- data.frame(Age, Drug, Effect)
head(df)

# 三维列联表（Age × Drug × Effect）
table_array <- table(df$Drug, df$Effect, df$Age)
table_array

# 使用mantelhaen.test()
result <- mantelhaen.test(table_array)
result


library(ggplot2)
ggplot(df, aes(x = Drug, fill = Effect)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Age) +
  labs(y = "Proportion") +
  theme_minimal()

输出：

, ,  = Middle

   
    Improved Not Improved
  A       15           84
  B       51           59

, ,  = Old

   
    Improved Not Improved
  A        7           63
  B       60           18

, ,  = Young

   
    Improved Not Improved
  A       52           24
  B        5           62 
	Mantel-Haenszel chi-squared test with continuity correction

data:  table_array
Mantel-Haenszel X-squared = 12.072, df = 1, p-value = 0.000512
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.4232736 0.8208328
sample estimates:
common odds ratio 
        0.5894378

从输出可以看到，在Middle和Old组药物B更有效，Young组则是药物A，而检验的结果p为0.0005，说明在调查年龄分组后，药物与疗效的关系十分显著，而公共比值则意味着使用药物B的患者获得改善的几率是药物A的0.59倍。