How to Sample from a Randomisation Distribution?

Question

I will use an example to explain my problem, but I think that it can be generalised to any randomization testing that does not allow one to create the null distribution by enumerating of all the possible orders of the data (that is the number of possible data permutation is too large to make the complete enumeration impossible), and a sampling of the randomised data should be made.

Example: Differential Correlation

Let be $\mathbf x_1$ and $\mathbf y_1$, both of length $n_1$, a pair of vectors that represents the age and weight of $n_1$ male individuals, and $\mathbf x_2$ and $\mathbf y_2$, both of length $n_2$, a pair of vectors that represents the age and weight of $n_2$ female individuals. The number of males may be different than the number of females.

Let $\rho_1$ be the Pearson's correlation coefficient between $\mathbf x_1$ and $\mathbf y_1$, and $\rho_2$ the Pearson's correlation coefficient between $\mathbf x_2$ and $\mathbf y_2$. To reject the null hypothesis that $\rho_1=\rho_2$ using a permutation test I should create random vectors where paired elements of $\mathbf x_1/\mathbf y_1$ are swapped with paired elements of $\mathbf x_2/\mathbf y_2$, in order to remove the relationship between the measurements and the labels (male/female).

The complete enumeration of all possible order of the data would be equal to all the possible order of the vector $\mathbf x (\mathbf y)$ built by concatenating the two vectors $\mathbf x_1 \mathbf x_2 (\mathbf y_1 \mathbf y_2)$, that is $(n_1+n_2)$! This makes my complete enumeration impossible, and I need to sample $N$ configurations from all the possible configuration of $\mathbf x$, of length $n_1+n_2$.

The problem: How many elements are swapped (in average) in each data configuration?

Each possible configuration will swap a random number $R$ of (paired) elements from $\mathbf x_1/\mathbf y_1$ to $\mathbf x_2/\mathbf y_2$, and I would like to know if $R$ follows any known distribution.

Let then $\mathbf v$ be a vector of length $n_1+n_2$, built by concatenating a vector $\mathbf v_1$, that includes $n_1$ elements set to $0$ representing $\mathbf x_1/\mathbf y_1$, and $\mathbf v_2$ that includes $n_2$ elements set to $1$ representing $\mathbf x_2/\mathbf y_2$. Let then shuffle the element of $\mathbf v$ for $N$ times, and count, at each step $i$ (with $i \in [i, N]$), how many elements $c_i$ of $n_1$ are now set to $1$, that is, how many elements were swapped at that permutation step. Let then build the distribution of my permuted elements $c$.

These steps are performed by the following function, while the histogram will show the distribution:

count.permuted.elements <- function(n1, n2, N=10000) {
    c <- numeric(N)
    for (i in 1:N) {   
        v <- c(rep(0, n1), rep(1, n2))
        v <- sample(v)
        c[i] <- sum(v[1:n1])
    }

    c
}

hist(count.permuted.elements(n1, n2))

I have run some simulations with different values of $n_1/n_2$, and what I obtained (see plots below) are normal distributions, as also confirmed by the function fitdist (in the R package fitdistrplus). Let me underline that the right tail is limited by the length of the shortest vector, that is, I can't permute more elements than those included in the shorter vector (see the plot on the bottom right for an example).

The problem (restated): How to estimate the parameters of the normal distributions ($\mu$ and $\sigma$) showed above, by using only $n_1$ and $n_2$?

After running several simulation, varying $n_1/n_2$, it seems to me that the mean $\mu$ is a function of the length of the smaller vector (let's say $n_1$) and the square root of the ratio between the two lengths, that is:

$$\mu = n_1 \cdot \left(-0.604 \cdot \sqrt{\frac{n_1}{n_2}} + 1.101\right)$$

However, I can't find a way to calculate $\sigma$ --and I am not sure that what I did to evaluate $\mu$ is very correct.

My question is: does anyone know how to evaluate, in a formal or empirical way, the parameters of the normal distributions ($\mu$ and $\sigma$) build by counting how many elements are swapped (in average) between the two arrays in each data configuration?

Thank you very much for your help!

Answer 1

You have $n_1$ $0$'s and $n_2$ $1$'s, so probability of sampling $1$'s from the whole sample is $p = \frac{n_2}{n_1+n_2}$. You randomly sample from it $n_1$ cases, where probability of sampling $1$'s is $p$, when $n_1 \ll n_2$ this resembles binomial distribution with parameters $n_1$ and $p$.

However your data does not follow binomial distribution since you are sampling without replacement from finite population. Your data follows a hypergeometric distribution. In this case mean and variance are

$$ \mu = n_1 \frac{n_2}{n_1+n_2} \\ \sigma^2 = n_1 \frac{n_2}{n_1+n_2} \frac{n_1}{n_1+n_2} \frac{n_2}{n_1+n_2-1} $$

Below you can see data created with your function with added pmf of binomial (red curve) and hypergeometric (blue curve) distribution.