I want to generate the 1000 tables from an existing contingency table and estimate the exact p-value by Monte-Carlo. To do so, I know that the contingency tables have to satisfy the condition S = {p(T) ≤ p (array-observe)}. Then, I have to compute the probabilities of unique tables among the 1000 generated.
I have started a part of the code, but I am stuck in generating the tables respecting the condition mentioned above. I looking for a command that could to so without using the command fisher.test.
Here is my code, where infarctus.tables corresponds to the generated tables from a contingency table and infarctus is my data.
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 )
Any idea on how to do so?
Base R includes function r2dtable. From the documentation.
Description
Generate random 2-way tables with given marginals using Patefield's algorithm.
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
Created on 2022-10-02 with reprex v2.0.2
To generate 1000 matrices with the given marginals, run
n <- 1000L
infarctus.tables <- r2dtable(n, rsums, csums)