The R function chisq.test has an simulate.p.value argument that uses permutation for resampling. I know how to implement this kind of simulation. However, I was wondering if it's also possible to do the following:
Let's say we have a sample $X_1, X_2, ..., X_n$ that distributes across disjoint categories $A$, $B$ and $C$ (i.e., $X_i$ is either $A$, $B$, or $C$). Now I sample (with replacement) from my sample and compute the $\chi^2$-statistic. Both steps are repeated many times, and I obtain the simulated $p$-value by computing the proportion of simulated $\chi^2$-statistic exceeding the original $\chi^2$-statistic.
The $p$-values are very different compared to using simulate.p.value. Certainly, we are comparing different approaches, so there should be no surprise if there are differences in the simulated $p$-values. However, the magnitude of difference makes me wonder whether there is a flaw in my reasoning.
