对于有多个特征的数据,我们一般的处理方式是构建特征函数,计算每个特征向量的系数,从而将其影响纳入到研究量中,但对于简单的问题,也这样做的话未免有点小题大做。这时我们可以考虑用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倍。