检索卡方检验的 monte carlo 模拟值

Question

我正在尝试绘制卡方中的零分布图 test.In R 可以使用以下代码进行 monte carlo 模拟以获得经验 p 值：

chisq.test(d,simulate.p.value=TRUE,B=10000)

但它没有 return 分布图。有没有办法让 R 达到 return 测试的模拟值？

Answer 1

如果您查看 chisq.test 的函数定义（大约 capture.output(chisq.test) 的第 56 行），您将进入模拟部分：

    if (simulate.p.value && all(sr > 0) && all(sc > 0)) {
        setMETH()
        tmp <- .Call(C_chisq_sim, sr, sc, B, E)
        STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE))
        PARAMETER <- NA
        PVAL <- (1 + sum(tmp >= almost.1 * STATISTIC))/(B + 
            1)
    }

这是在调用 C 函数。首先生成一些虚拟数据

## Some data
x <- as.table(rbind(c(762, 327, 468), c(484, 239, 477)))
dimnames(x) <- list(gender = c("F", "M"),
                party = c("Democrat","Independent", "Republican"))

然后抓住你需要的一点

sr <- rowSums(x)
sc <- colSums(x)
n <- sum(x)
E <- outer(sr, sc, "*")/n
v <- function(r, c, n) c * r * (n - r) * (n - c)/n^3
V <- outer(sr, sc, v, n)
dimnames(E) <- dimnames(x)
B = 2000
tmp <- .Call(stats:::C_chisq_sim, sr, sc, B, E)
STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE))
almost.1 <- 1 - 64 * .Machine$double.eps                                                  
PVAL <- (1 + sum(tmp >= almost.1 * STATISTIC))/(B + 1)

变量tmp包含您想要的输出。变量 PVAL 匹配

的输出

chisq.test(x, simulate.p.value = T, B=2000)$p.value

请注意，我使用了 :::，因为函数 C_chisq_sim 未从统计信息中导出。

检索卡方检验的 monte carlo 模拟值

Retrieving the monte carlo simulation values for chi square test

simulation

r

chi-squared