parametric bootstrap for low sample sizes

Question

I believe that this question is sufficiently different from previous related ones to warrant a new post. (I apologize if it has been answered already)

I need to decide between various resampling methods to "best" (highest power and correct type-I error to reject H0 for the right reasons) evaluate auto correlation in binary sequences of small to moderate lengths (20-50) with relatively low incidence rate. The time series are too short to justify any normal approximations, so to e.g. to decide whether or not the observed auto correlation is "real", I would like to generate my empirical Null distribution under the assumption of no auto correlation. With that goal in mind I can either

1. shuffle the sequence repeatedly (fixed number of 0s and 1s), or
2. estimate the probability $\hat{p}$ of an event and repeatedly draw from the respective Bernoulli process.

The latter is -I think- known as the parametric bootstrap. I struggle with the correct choice for the following reasons:

(i) The estimate for the true Bernoulli probability will be rather poor/noisy for low smple sizes. That bias will be part of the generated empirical NULL. (ii) The permutation procedure 1 will generate test statistics with considerably less variation than procedure 2. In fact, the distribution will tend to generate few discrete levels.

What line of reasoning will defend either choice ?

Here is an example:

x=c(0,1,0,0,0,1,1,0,0,1,0,0,0,0,1,0)
#observed acf:
aObs = acf(x,lag.max=1,plot=F)$acf[2]
    a1=a2=rep(NA,100)
    for (i in 1:100){
      a1[i]=acf(sample(x),lag.max=1,plot=F)$acf[2]
      a2[i]=acf(rbinom(length(x),1,mean(x)),lag.max=1,plot=F)$acf[2]
}
hist(a1);abline(v=aObs,lty=2,col=2)
hist(a2);abline(v=aObs,lty=2,col=2)

I am adding more code as a reply to my motivation. Let us begin with computing the exact p-value of 0.01011 from a Fisher test:

  Convictions <-
  matrix(c(2, 8, 10, 3),
         nrow = 2,
         dimnames =
           list(c("Dizygotic", "Monozygotic"),
                c("Convicted", "Not convicted")))
fisher.test(Convictions, alternative = "less")

Instead, could we not simply simulate (and hence NOT condition on the margins)

p=sum(Convictions[,"Convicted"])/sum(Convictions)
N = rowSums(Convictions)
ConvictionsSim = Convictions
OR0=prod(diag(Convictions)) / prod(as.vector(Convictions)[c(2:3)])

OR = rep(NA,1000)
for (i in 1:1000){
  ConvictionsSim[1,1] = rbinom(1,N[1],p=p)
  ConvictionsSim[1,2] = N[1]-ConvictionsSim[1,1]
  ConvictionsSim[2,1] = rbinom(1,N[2],p=p)
  ConvictionsSim[2,2] = N[2]-ConvictionsSim[2,1]
  OR[i] = prod(diag(ConvictionsSim)) / prod(as.vector(ConvictionsSim)[c(2:3)]) 
}
mean(OR<OR0)

which gives me a very different p-value of 0.004. Which one is "correct"?

Answer 1

Why don't you use 2), only via the beta-binomial scheme:

http://www.cs.cmu.edu/~10701/lecture/technote2_betabinomial.pdf

That is, the probability of success, $\hat{p}$, is drawn from the posterior distribution separately for each simulated sequence.

The advantages are:

1) Low/high incidence rate is not a problem, even if you have 0% or 100% of successes.

2) Uncertainty in estimating the probability of success, $\hat{p}$, is incorporated in the simulated series.

3) You don't need to rely on (Normal) approximations anywhere.

I then looked up some information here and there, and it looks like a permutation test is the way to go. Suppose you have two samples: one with autocorrelation, the other one without, but the probability of success is the same. Under $H_0$, all those observations are i.i.d. and therefore exchangeable, which is the main assumption for the permutation test.

To compare this to bootstrap, consider a situation when Sample1 is drawn from distribution $G$ and Sample2 is drawn from $F$. The null hypothesis is the equality of means:

$H_0: E_G = E_F$

Then, even when $H_0$ is true the observations are not exchangeable because their marginal distribution doesn't have to be the same. Therefore, one can't use permutation and has to bootstrap.