R中列联表的蒙特卡罗模拟

Question

我想从现有的列联表中生成1000个表，并通过蒙特卡洛估计精确的p值。为此，我知道列联表必须满足条件S={p(T(≤p(数组观察(}。然后，我必须计算生成的1000个表中唯一表的概率。

我已经开始了部分代码，但我在生成符合上述条件的表时遇到了困难。我正在寻找一个可以这样做的命令，而不使用命令fisher.test.

这是我的代码，其中infarctus.tables对应于从列联表生成的表，infarctus是我的数据。

infractus <- matrix(c(6, 4, 2, 0, 7, 4, 1, 0, 2,2, 3, 5, 2, 5, 3, 2), nrow = 4, byrow = T)
infarctus.tables ### how do I generate this?
##probabilities of unique tables among the 1000 generated
tables.prob= lapply( unique(infarctus.tables), 
FUN=function(x) { exp( a1 - sum(lgamma( x + 1) ) ) })
hist( unlist(tables.prob))
obs.prob= which(unlist(tables.prob) <= 0.0164141415 )

你知道怎么做吗？

Answer 1

Base R包含函数r2dtable。来自文档。

描述
使用Patefield算法生成具有给定边距的随机双向表。

infractus <- matrix(c(6, 4, 2, 0, 7, 4, 1, 0, 2,2, 3, 5, 2, 5, 3, 2), nrow = 4, byrow = T)
rsums <- rowSums(infractus)
csums <- colSums(infractus)
n <- 5
r2dtable(n, rsums, csums)
#> [[1]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    3    5    4    0
#> [2,]    6    3    1    2
#> [3,]    3    4    3    2
#> [4,]    5    3    1    3
#> 
#> [[2]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    8    3    1    0
#> [2,]    1    6    2    3
#> [3,]    4    3    4    1
#> [4,]    4    3    2    3
#> 
#> [[3]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    4    3    3    2
#> [2,]    5    4    2    1
#> [3,]    6    1    3    2
#> [4,]    2    7    1    2
#> 
#> [[4]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    4    5    2    1
#> [2,]    4    6    1    1
#> [3,]    7    1    3    1
#> [4,]    2    3    3    4
#> 
#> [[5]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    3    4    4    1
#> [2,]    5    4    2    1
#> [3,]    7    2    1    2
#> [4,]    2    5    2    3

^{创建于2022-10-02，reprex v2.0.2}

要生成1000个具有给定边距的矩阵，请运行

n <- 1000L
infarctus.tables <- r2dtable(n, rsums, csums)