# Not sure if standard error of p-values makes sense in Fisher Exact Test

I am working on implementing a Fisher Exact Test for some unemployment and wage data. The idea is to describe two populations (one receiving an assistance program (the "treatment") and one not receiving it) via a summary statistic.

The first statistic that I am using is the difference in average wage-level post-treatment between the two populations.

To calculate a p-value for such a set-up, the idea is to make a large number $N$ of random permutations of the assignment vector (the list of 0's and 1's that indicates whether a given observation belongs to the control population or not). I calculate the summary statistic for each of these random assignments, to get a distribution of summary statistics in counterfactual experiments where the control group had been different.

Since I also have the summary statistic for the observed assignment vector (the actual vector of 0's and 1's for the experiment that was truly observed), I can then count the number of simulated summary statistics which are more extreme than the observed summary statistic.

This proportion out of all of my simulated trials serves as the estimated p-value for my summary statistic.

My question is as follows: is there a standard way to get a standard error for such an estimated p-value? Obviously, I can calculate the Monte Carlo standard error for the summary statistics in my simulations, but since all of the simulations are used for just one single p-value calculation, it''s not clear how to get the standard error for that.

I had the following thought about what the standard error might be in this case.

In usual Monte Carlo, we have some function $f$ that we compute at each of the simulated draws $x_{i}$ (where here, the $x_{i}$ are understood to be the assignment vectors). If I define:

$$f(x_{i}) = \mathbb{1}_{ |stat(x_{i})| > |stat(x_{obs})| }$$

then it seems that the p-value I calculate is given by

$$\hat{p} = \frac{1}{N}\sum_{i=1}^{N}f(x_{i}) = \frac{\textrm{# more extreme}}{\textrm{total samples}}$$

And following the usual Monte Carlo formulas, would it then make sense to write the variance of the estimate as:

$$Var(\hat{p}) = \frac{1}{N}\sum_{i=1}^{N}[ f(x_{i}) - \hat{p} ]^{2}$$

and then take the square root to obtain the standard error?

The reason this confuses me is that for each $i$, $f(x_{i})$ will be binary, either the computed statistic was more extreme in that iteration or it wasn't. It seems like it would be error prone to sum up a bunch of binary things like that to estimate the variance in a p-value, but that could just be my unfamiliarity with this method.

Can anyone confirm that this is right? Also, if I exposed any other ignorance about what I am doing here, corrective comments are appreciated.

You have correctly computed the variance of $f$, not the variance of $\hat{p}$. To estimate the variance of $\hat{p}$ just go one step further:

$$\text{var}(\hat{p})= \text{var}\left(\frac{1}{N}\sum_{i=1}^{N}f(x_{i})\right) = \frac{1}{N} \text{var}(f).$$

Moreover, your formula simplifies greatly: it's an easy algebraic step to derive the estimator

$$\widehat{\text{var}(f)} = \hat{p}(1 - \hat{p}).$$

Whence, to find the standard error of $\hat{p}$, divide $\hat{p}(1 - \hat{p})$ by $N$ and take the square root.

After all, the distribution of $\hat{p}$ is binomial and these are the familiar formulas for sampling from a binomial distribution.

In general, if you are monte-carlo-ing a p-value, the only standard error that makes sense is the monte-carlo error in estimating it, and then only to show that you have done enough permutations. Under the null hypothesis, the p-value is just Uniform(0,1). Unless the p-value is really extreme, most people will be happy with you just reporting it and saying how many permutations you did - I definitely would not report a standard error in any results table.

• I totally agree. It threw me for a loop when my adviser asked for it. I thought the same thing... it doesn't give a confidence that the p-value is telling us something, just that we did enough iterations to believe our Monte Carlo approximation to FET is good. – ely Apr 7 '12 at 3:28
• @EMS, a Monte Carlo calculation is an approximation. It therefore behooves one to know how accurate that approximation is, which is why it is reasonable for your advisor to ask for some assessment of the quality of that approximation. More to the point, once you know how to quantify the connection between the number of MC iterations and precision of the p-value, you can determine (in advance) how many iterations may be necessary. – whuber Apr 7 '12 at 20:27
• Right, I get that. But in this case, it should have been pretty obvious that we were using far more iterations in the m.c. simulation than necessary. The s.e. I am computing is smaller than 10^(-9) for a p-value that realistically only needs a few digits of accuracy. If the simulation wasn't for a p-value, but for something truly needed many digits of accuracy, I would totally see the merit in needing to know a standard error, as has been the case for the many m.c. simulations I've created for other studies. – ely Apr 7 '12 at 21:03
• I should also add that the code is already very fast too. So learning that we could reduce the iterations doesn't allow us to do things faster. We'd be saving a fraction of a second on something we'll only run a few dozen times, if that. I think this s.e. was requested more out of a force of habit than anything else, which is totally fine and probably a good habit for me to develop too. – ely Apr 7 '12 at 21:06
• Yes, indeed it is blatantly obvious. Again, this is why I was puzzled. It's effectiveness of a particular unemployment program that is known to be extremely effective. The test statistic in question is extremely different for two populations. Just histogramming the outcomes shows you it's unreasonable under the Fisher Exact Test, which assumes a sharp null hypothesis that the treatment is ineffective for literally all participants. – ely Apr 8 '12 at 0:04