对于有多个特征的数据，我们一般的处理方式是构建特征函数，计算每个特征向量的系数，从而将其影响纳入到研究量中，但对于简单的问题，也这样做的话未免有点小题大做。这时我们可以考虑用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)
resultlibrary(ggplot2)
ggplot(df, aes(x = Drug, fill = Effect)) +geom_bar(position = "fill") +facet_wrap(~ Age) +labs(y = "Proportion") +theme_minimal()

输出：

, ,  = MiddleImproved Not ImprovedA       15           84B       51           59, ,  = OldImproved Not ImprovedA        7           63B       60           18, ,  = YoungImproved Not ImprovedA       52           24B        5           62 Mantel-Haenszel chi-squared test with continuity correctiondata:  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倍。