[R] 'simulate.p.value' for goodness of fit
Bob
thebobs at mistral.co.uk
Thu Jan 17 19:59:26 CET 2008
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
More information about the R-help
mailing list