使用mantelhaen.test()时遇NA/NaN/Inf外部函数错误及整数溢出问题

我来帮你拆解这个问题——我之前在处理大样本CMH检验时也碰到过几乎一模一样的报错，咱们从原因到解决方法一步步梳理：

触发原因分析 #

这个问题其实是两个因素共同作用的结果：

1. 整数溢出导致的NA值
   你看到的警告NAs produced by integer overflow是关键线索。R默认的整数类型是32位，最大只能存储约21亿的数值。当你的数据有10万行时，计算ntot * rowsums（总样本数乘以各行列的频数和）很容易超过这个上限，溢出后变成NA。这些NA值被传递到后续的qr.default函数中，就触发了NA/NaN/Inf in foreign function call的错误。

2. 极端稀疏的分层列联表
   即使你用na.omit处理了缺失值，某些性别分组（sex的某个水平）里，education和score.grouped的交叉表可能存在全0的行/列，或者某层的总频数极小。这种情况会导致计算过程中出现奇异矩阵，进而让qr.solve（CMH检验内部用来计算统计量的函数）调用Fortran代码时失败。

而且有意思的是，这两个问题只会在特定变量组合下出现——比如你提到的education和imc.cl，可能它们的交叉表在某些分层里刚好出现了极端稀疏+大数值相乘的组合，而其他变量没有。

解决方法 #

1. 先解决整数溢出问题 #

最直接的方式是把频数转换为64位数值型（而非默认的32位整数），避免溢出。你可以手动构建分层列联表并转换类型后再传入检验：

# 先生成分层列联表，同时自动移除缺失值
xtab <- xtabs(~ education + score.grouped + sex, data = db, na.action = na.omit)
# 转换为数值型（64位），避免整数溢出
xtab <- as.matrix(xtab)
# 执行CMH检验
mantelhaen.test(xtab)

2. 检查并修复稀疏分层 #

遍历每个性别分组，查看交叉表是否存在全0的行/列：

# 遍历所有性别分组
for (s_level in unique(db$sex)) {
  # 提取当前分组的子数据
  sub_data <- db[db$sex == s_level, ]
  # 生成交叉表
  sub_tab <- table(sub_data$education, sub_data$score.grouped)
  
  cat("=== 性别分组：", s_level, " ===\n")
  print(sub_tab)
  
  # 检查全0行/列
  row_zero <- any(rowSums(sub_tab) == 0)
  col_zero <- any(colSums(sub_tab) == 0)
  
  if (row_zero) cat("⚠️ 该分组存在全0的教育水平行\n")
  if (col_zero) cat("⚠️ 该分组存在全0的得分组列\n")
}

如果发现全0的行/列，可以：

• 合并频次极低的因子水平（比如把education中样本量<10的水平合并为“其他”）
• 如果某个性别分组的样本量极小（比如只有几行），可以考虑移除该分组（前提是不影响分析结论）

3. 换用更鲁棒的第三方函数 #

推荐使用epitools包中的cmh_test函数，它对大样本和稀疏数据的处理更友好，还会给出更详细的错误提示：

# 安装并加载包
install.packages("epitools")
library(epitools)

# 执行CMH检验
cmh_test(xtabs(~ education + score.grouped + sex, data = db, na.action = na.omit))

这个函数内部已经处理了整数溢出的问题，而且如果某个分层有问题，会明确告诉你具体是哪个分组导致的，方便定位修复。

内容的提问来源于stack exchange，提问作者Dan Chaltiel