mantelhaen.test() 的外部函数错误中的 NA/NaN/Inf

Question

我有一个 100k 行数据框，我想在其上计算 Cochran-Mantel-Haenszel 检验。

我的变量是教育水平和分位数的计算分数，我的分组变量是性别，代码行如下所示：

mantelhaen.test(db$education, db$score.grouped, db$sex)

此代码抛出此错误和警告：

qr.default(a, tol = tol) 中的错误：外部函数调用中的 NA/NaN/Inf (arg 1)
另外：警告信息：In ntot * rowsums : NAs generated by integer overflow

错误似乎是由我的第一个变量引起的，因为在测试的 7 个变量上，我只有其中 2 个出现问题，这似乎没有明显的共同点。

在抛出错误的变量和不抛出错误的变量之间，缺失值和因子水平似乎没有区别。我尝试了完整的案例（na.omit），但问题仍然存在。

是什么触发了这个错误？是什么意思？
我怎样才能摆脱它？

有趣的帖子：R: NA/NaN/Inf in foreign function call (arg 1)，What is integer overflow in R and how can it happen?

ADDENDUM：这是str的结果（失败是education和imc.cl）：

str(db[c("education","score.grouped","sex", ...)])
'data.frame':   104382 obs. of  7 variables:
 $ age.cl: Ord.factor w/ 5 levels "<30 ans"<"30-40 ans"<..: 5 2 1 1 3 4 2 3 4 4 ...
  ..- attr(*, "label")= chr "age"
 $ emploi2          : Factor w/ 8 levels "Agriculteurs exploitants",..: 3 5 6 8 8 8 8 3 3 3 ...
  ..- attr(*, "label")= chr "CSP"
 $ tabac            : Factor w/ 4 levels "ancien fumeur",..: 4 1 4 4 3 4 4 1 4 4 ...
  ..- attr(*, "label")= chr "tabac"
 $ situ_mari2       : Factor w/ 3 levels "Vit seul","Divorsé, séparé ou veuf",..: 3 2 1 1 1 3 1 3 2 3 ...
  ..- attr(*, "label")= chr "marriage"
 $ education        : Factor w/ 3 levels "Universitaire",..: 1 1 1 1 3 1 1 1 1 1 ...
 $ revenu.cl        : Factor w/ 4 levels "<1800 euros/uc",..: 3 4 2 NA 4 1 1 4 4 1 ...
 $ imc.cl           : Ord.factor w/ 6 levels "Maigre"<"Normal"<..: 2 2 1 2 3 1 3 2 2 3 ...
  ..- attr(*, "label")= chr "IMC"

EDIT：通过深入函数内部，错误和警告是由对qr.solve 的调用引起的。我对此一无所知，但我会尝试更深入地研究
EDIT2：在qr.solve 内部，错误是由Fortran 调用.F_dqrdc2 引发的。这远远超出了我的水平，我的鼻子开始流血了。
EDIT3：我尝试head我的数据以找出原因：

db2 = db %>% head(99787)   #fails at 99788
db2 = db %>% tail(99698)   #fails at 99699
mantelhaen.test(db2$education, db2$score.grouped, db2$sex)

这给我的信息不多，但也许它可以给你。

Answer 1

我能够通过扩大数据集来复制问题。

set.seed(101); n <- 500000
db <- data.frame(education=
                   factor(sample(1:3,replace=TRUE,size=n)),
                 score=
                   factor(sample(1:5,replace=TRUE,size=n)),
                 sex=
                   sample(c("M","F"),replace=TRUE,size=n))

在此之后，mantelhaen.test(db$education, db$score, db$sex) 给出报告的错误。

值得庆幸的是，真正的问题不在 QR 分解代码的内部：而是在 QR 分解之前设置矩阵时发生。有两个计算，ntot*colsums 和 ntot*rowsums，溢出了 R 的整数计算能力。通过创建函数的修改版本，有一种相对简单的方法可以解决此问题：

• 复制源码：dump("mantelhaen.test",file="my_mh.R")
• 编辑源代码
  • 升。 1：修改函数名称为 my_mantelhaen.test（避免混淆）
  • 第 199 和 200 行：将 ntot 更改为 as.numeric(ntot)，在溢出发生之前将整数转换为双精度
• source("my_mh.R") 读入新功能

现在

my_mantelhaen.test(db$education, db$score, db$sex)  

应该可以。 您应该绝对针对旧功能测试新功能，以确保您得到相同的答案。

现在发布到R bug list，我们看看会发生什么...

更新 2018 年 5 月 11 日：这是fixed in the development version of R（待定为 3.6）。