'simulate.p.value' for goodness of fit

R Help on 'chisq.test' states that
    "if 'simulate.p.value' is 'TRUE', the p-value is computed by Monte
     Carlo simulation with 'B' replicates.

     In the contingency table case this is done by random sampling  
from
     the set of all contingency tables with given marginals, and works
     only if the marginals are positive...

     In the goodness-of-fit case this is done by random sampling from
     the discrete distribution specified by 'p', each sample being of
     size 'n = sum(x)'."

The last paragraph suggests that in the goodness-of-fit case, if p  
gives the

expected probability for each cell, this random sampling is  
multinomial.

Unfortunately, as the following examples reveal, the sampling model is
neither

multinomial nor hypergeometric - at least when it is applied to a 4- 
fold
table.

This is rather sad as some people assume that R's chisq.test  
function can

perform a Monte Carlo test of X-squared, employing a multinomial  
model - in

other words, assuming that your data are a single random sample.


### Example 1.

x=matrix(c(1,2,3,4),nc=2)
# To begin with, let us apply the large-sample approximations
chisq.test(x,correct=TRUE)$p.value

[1] 0.6726038
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(x,  
correct = TRUE)

chisq.test(x,correct=FALSE)$p.value

[1] 0.7781597
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(x, correct =
FALSE)

# So let us apply a 2-tailed test of O.R.=1, using a  
hypergeometric model
fisher.test(x)$p.value

[1] 1

# This should also apply a hypergeometric model
chisq.test(x,simulate.p.value=TRUE,B=500000)$p.value

[1] 1

# Now we work out the expected probability for each cell
p=outer(c(sum(x[1,]),sum(x[2,])),c(sum(x[,1]),sum(x[,2])))/sum(x)^2
# But this applies a hypergeometric model, presumably because p is  
not
scalar
chisq.test(x,p=p,simulate.p.value=TRUE,B=500000)$p.value

[1] 1

# This seems to do something different,
# at any rate it is much slower, and needs more memory
chisq.test(x[1:4],p=p[1:4],simulate.p.value=TRUE,B=10000)$p.value

[1] 1

# Which would appear to be using the same model as above

# Now let us apply an X2 test using a multinomial model
# (The code for this x2.test function is in Appendix 1, below.)
x2.test(x,R=200000)

with cc P =  0.7316812
conventional-P =  0.8838786
mid-P =  0.8423058

# All of these P-values are higher than those given by the Chi- 
squared

approximation,

# but they certainly do not equal 1.
# But is this is an artefact of our very small sample?




### Example 2.
# Let us try a larger sample
x=matrix(c(56,35,23,42),nc=2)

# First we apply the asymptotic model
chisq.test(x,correct=TRUE)$p.value

[1] 0.002222425

chisq.test(x,correct=FALSE)$p.value

[1] 0.001276595

# Now for the hypergeometric (fixed margin totals model)
fisher.test(x)$p.value

[1] 0.001931078

chisq.test(x,simulate.p.value=TRUE,B=500000)$p.value

[1] 0.001913996

p=outer(c(sum(x[1,]),sum(x[2,])),c(sum(x[,1]),sum(x[,2])))/sum(x)^2
chisq.test(x,p=p,simulate.p.value=TRUE,B=500000)$p.value

[1] 0.001891996

Next comes what we had hoped to be a multinomial test
chisq.test(x[1:4],p=p[1:4],simulate.p.value=TRUE,B=10000)$p.value

[1] 0.01639836

# This is obviously not the same hypergeometric model as used for  
a > #

chi-squared test.

# The P-value is about 10x of the approximate tests (above)
#  or the exact tests (below).

x2.test(x,R=200000)

with cc P =  0.002059990
conventional-P =  0.001184994
mid-P =  0.001172494

# Whatever that chi-squared test model IS, it is certainly not
multinomial!
# Could it possibly be Poisson and, if so, why???

######## Appendix 1:

# We have used these  functions to do a 2x2 multinomial test of X2:

x2=function(y,cc=FALSE){
y=y*1.;n=sum(y);C=cc*n/2
a=y[1];b=y[2];c=y[3];d=y[4]
ab=a+b;cd=c+d;ac=a+c;bd=b+d
D=ab*cd*ac*bd
if(D==0)x2=NA else x2=n*(abs(a*d-b*c)-C)^2/D
x2}

x2.test=function(x,R=5000){
n=sum(x)
p=outer(c(sum(x[1,]),sum(x[2,])),c(sum(x[,1]),sum(x[,2])))/n/n
Q=sort(apply(rmultinom(R,n,p),2,x2))
q=x2(x,cc=TRUE)
pl=rank(c(q,Q),ties.method='max')[1]/(length(Q)+1)
pe=sum(c(q,Q)==q)/(length(Q)+1);pu=1-pl+pe
cat('with cc P = ',pu,'\n')
q=x2(x)
pl=rank(c(q,Q),ties.method='max')[1]/(length(Q)+1)
pe=sum(c(q,Q)==q)/(length(Q)+1);pu=1-pl+pe
cat('conventional-P = ',pu,'\n')
cat('mid-P = ',pu-pe/2,'\n')}

Bob