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

Question

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.

Answer 1

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.

Answer 2

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.