Help with a one-sample permutation z-test

Question

This is related to another question I asked recently. To recap:

[I had 30 people call a number and then roll a 5 sided die. If the call matches the subsequent face then the trial is a hit, else it is a miss. Each subject completes 25 trials (rolls) and thus, each participant has a score out of 25. Since the die is a virtual one, it cannot be biased. Before the experiment was conducted I was going to compare the subjects score with a one-sample t-test (compared to mu of 5). However I was pointed towards the more powerful z-test, which is appropriate because we know the population parameters for the null hypothesis: that everyone should score at chance. Since NPQ means the binomial approximates to the normal or Gaussian distribution, we can use a parametric test. So I could just forget about it all and go back to the t-test I planned to use, but it seems to me that although the z-test is not often used in real research it is appropriate here. That was the conclusion from my previous question. Now I am trying to understand how to use resampling methods (either permutation or bootstrap) to compliment my parametric analysis.]

Okay. I am trying to program a one-sample permutation z-test, using the DAAG package onet.permutation as inspiration. This is as far as I've got:

perm.z.test = function(x, mu, var, n, prob, nsim){
  nx <- length(x)
  mx <- mean(x)
  z <- array(, nsim)
  for (i in 1:nsim) {
    mn <- rbinom(nx*1, size=n, prob=prob)
    zeta = (mean(mn) - mu) / (sqrt(var/nx))
    z[i] <- zeta
  }
  pval <- (sum(z >= abs(mx)) + sum(z <= -abs(mx)))/nsim
  print(signif(pval, 3))
}

Where: x is the variable to test, n is the the number of trials (=25) and prob is the probability of getting it correct (=.2). The population value (mu) of the mean number correct is np. The population standard deviation, var, is square-root(np*[1-p]).

Now I guess this compares x to an array composed of randomly generated binomial sample. If I centre x at 0 (variable-mu) I get a p-value. Can somebody confirm that it is doing what I think it is doing?

My testing gives this:

> binom.samp1 <- as.data.frame(matrix(rbinom(30*1, size=25, prob=0.2), 
                                     ncol=1))
> z.test(binom.samp1$V1, mu=5, sigma.x=2)
data:  binom.samp1$V1 
z = 0.7303, p-value = 0.4652

> perm.z.test(binom.samp1$V1-5, 5, 2, 25, .2, 2000)
[1] 0.892

> binom.samp1 <- as.data.frame(matrix(rbinom(1000*1, size=25, prob=0.2),
                                     ncol=1))
> perm.z.test(binom.samp1$V1-5, 5, 2, 25, .2, 2000)
[1] 0.937

Does this look right?

UPDATE:

Since this obviously doesn't do what I want, I do have another angle. This website offers this advice:

There is no reason whatsoever why a permutation test has to use any particular test statistic. Any test statistic will do! ... For one-sample or paired two-sample tests, in particular, for Wilcoxon signed rank tests, the permutations are really subsets. The permutation distribution choses an arbitrary subset to mark + and the complementary subset is marked -. Either subset can be empty.

What about an arbitrary subset with a one-sample z-test?

Answer 1

Caveat: I'm not sure I fully understand your question. With this in mind, your solution does, IMHO, not provide a one-sample permutation z-Test, as it does not use the original data while performing some re-labeling of the experimental units consistent with the Null hypothesis in the given experimental design. Actually I do not see how any re-labeling can be performed at all in your situation that appears to be one-sample-from-one-population with a test for the distribution's location parameter.

The onet.permutation() function you cite is, IMHO, misleadingly named as it refers to the test of two dependent samples. For each unit it randomizes which of its two values belongs to sample 1, and which one to sample 2. This is equivalent to randomizing the sign of the unit-wise difference between the two samples, as done by the lines

mn <- sample(c(-1, 1), n, replace = TRUE)
xbardash <- mean(mn * abs(x))

(x is the difference between two dependent samples) Your function does something else: it simulates new data and creates a simulated distribution that does not stem from the empirical data.

Answer 2

I don't think it is going a one sample Z; to me it looks like it is a test against a certain set of priors.

I'm confused, why are you doing a one sample Z using binomial data as your source data? You could simply create a distribution of N successes and see what quantile your actual data was in. However the above method doesn't look like a permutation test per-se to me; as your code doesn't actually permute the observed values between 1 and 0.

That being said, let me comment on your code - to me it looks as though z is defined as

zeta = (mean(mn) - mu) / (sqrt(var/nx))
z[i] <- zeta

Thus, each score in Z is like a Z score of the randomly created binomial vector using the priors you've selected as the null hypothesis. But then you compare that Z to abs(mx); where mx is defined as the mean of your observed binomial vector. At the very least this looks like a problem to me. Z scores should be either compared to some other Z score or means should be compared to means.

As I alluded to above, it is odd that you'd put all of this under a structure of a Z-test. The Z score is nominally a linear transform of the differences between means, as such the result of a test like this should be the same whether you use a Z score or simply look at the differences between means.

Moreover, what you are doing seems a like an attempt to test the observed value against some priors rather than an actual one sample permutation test. What you want to test against is something like permbinom (code provided below) where it could be the case for each observed value that it either was a success or it was not a success. This is in-line with Tukey's classic example with the lady who claimed she could tell whether tea or milk was added first. Critically different from your test, the assumption of this permutation test is that the null hypothesis is fixed at p = .5.

permbinom <- function(x)
{
  newx <- x
    nx <- length(x) 
    change <- rbinom(n=nx,size=1,prob=.5) 
    #This code is readable but inefficient
    #Swap the values between 1 and 0 if change == 1
    for (i in 1:nx)
    {
        if ((change[i] == 1) & (x[i] == 1)) {newx[i] <- 0} 
        if ((change[i] == 1) & (x[i] == 0)) {newx[i] <- 1} 
    }
  return(newx)
}
permtest <- function(x,nsim=2000)
{
   permref <- rep(NA,nsim)
   obsn <- sum(x)
   for (i in 1:nsim)
     {
      permref[i] <- sum(permbinom(x))
     }
     pval <- min(c(
              (sum(permref > obsn)*2/nsim),
              (sum(permref < obsn)*2/nsim)
                ))
  return(pval)
}

I'm not 100% confident regarding how I'm calculating the p-value here; so if someone would kindly correct me if I'm doing it wrong I'll incorporate that as an edit.

For reference, here is a faster permutation function for one-sample tests of binomial data.

permbinomf <- function(x)
{
  newx <- x
    nx <- length(x) 
    change <- rbinom(n=nx,size=1,prob=.5) 
    #This code is readable but inefficient
    #Swap the values between 1 and 0 if change == 1
newx <- x + change
newx <- ifelse(newx==2,0,newx)
  return(newx)
}

Edit: The question is also put forth, "What about an arbitrary subset with a one-sample z-test?". That would also work, assuming you had a large enough sample to subset. However, it would not be a permutation test, it would be more akin to a bootstrap.

Edit 2: Perhaps the most important answer to your question, is this: You are doing something acceptable (if you fix the Z vs mean computational error noted above), but you aren't doing what you think you are doing. You are comparing your results to results where the null hypothesis is true. This is essentially a Monte-Carlo simulation and if you correct the math (and I suggest you also simplify it) it is an acceptable technique for testing your hypothesis. Also note, my answer above is for a two-tailed test. As noted in the other question, you are ignoring the nesting of binomial observations under participants but independence isn't an assumption in a permutation or monte-carlo test so you should be more or less fine. Though, as also noted there you ignore the possibility that some people are doing better than chance and others are performing at chance.