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:

#observed acf:
aObs = acf(x,lag.max=1,plot=F)$acf[2]
    for (i in 1:100){

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)

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)]) 

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


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


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.

  • $\begingroup$ Thanks a lot for this very useful answer. Indeed, my problem of poor estimation of incidence rates is handled more gracefully by Bayesian type methods. And maybe that is all I need as a practical solution. However, my question is more conceptual: under what circumstances/assumptions would choice 2 be more sensible than choice 1? In other words, what settings lead naturally to a Fisher test type (fixed marginals) and when ist it more appropriate to allow for extra variation due to fluctuating total counts? $\endgroup$ Sep 4 '14 at 8:44
  • $\begingroup$ What I don't understand about 1) is why you want to shuffle it as opposed to sampling with replacement. If you sample with replacement, then the "less variation" issue will be taken care of and it seems to me both methods will produce similar results. $\endgroup$
    – James
    Sep 4 '14 at 14:22
  • $\begingroup$ Very sorry about my turtle speed of replies. I added a section in my question above that maybe shows the motivation. In the meantime I realized what a long and interesting controversy a question very similar to mine (on conditioning on the margins) has had, allow me to point out the following reference: "Testing the Equality of Two Independent Binomial Proportions", Roderick J. A. Little, The American Statistician, Vol. 43, No. 4 (Nov., 1989), pp. 283-288 $\endgroup$ Sep 26 '14 at 14:07
  • $\begingroup$ I think your simulation should have used the Fisher statistic that you used for the original data. That is, you simulate the data under H0, compute the statistic for each simulated data set, and see how often it is more extreme than the statistic for the original data. $\endgroup$
    – James
    Sep 26 '14 at 14:43
  • $\begingroup$ Thanks for the suggestion, replacing the last line my loop with OR[i] = fisher.test(ConvictionsSim, alternative = "less")$estimate yields a somewhat higher p-value of 0.008 but the fundamental question of whether or not the margins should be fixed for the baseline Null remains. $\endgroup$ Sep 27 '14 at 20:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.