5
$\begingroup$

I'm trying to generate a dataset with a pre-defined correlation between a normally distributed variable and a binary variable.

The method I had originally thought of was the following:

  1. Generate $X \sim Norm(0,1)$
  2. Generate $Y \sim Norm(0,1)$
  3. Generate $Q = \rho X + \sqrt{1-\rho^2}Y$, this will be the log-odds of success
  4. Generate $P = 1-\frac{1}{exp(Q) + 1}$, this is the probability of success
  5. Generate $U = Unif(0,1)$
  6. Generate $T = I(U < P) $

This method ensures that $Corr(X,P) = \rho$, however $Corr(X,T) \ne\rho$.

An alternative algorithm, replacing step 3 with:

  1. Generate $Q = \rho X$

provides similar results for $-1<\rho<1$, but in this second algorithm, we can take $\rho$ to be any value in $\mathbb{R}$, and still we have control of the correlation between $X$ and $T$.

I ran this algorithm through R using different values of $\rho$ and plotted $\rho$ against $Corr(X,T)$ (using the default methods in cor.test):

Plot of $\rho$ against $Corr(X,T)$

After altering $\rho$, it appears that the correlation between the continuous and binary variable is bounded, approximately $-0.8 \le Corr(X,T) \le 0.8$. While trying to figure out a relationship between $\rho$ and $Corr(X,T)$, I thought it looked like an arctangent and so the red line is the plot of $1.6*tan^{-1}(\rho)/\pi$, which doesn't quite match. when using the first algorithm (and limiting $-1<\rho<1$), the relationship between $\rho$ and $Corr(X,T)$ appears to be linear, with $Corr(X,T) = 0.43\rho$

My first question is whether there is any literature or sources on explicitly finding the relationship between $\rho$ and $Corr(X,T)$? That way I can predefine this correlation, rather than $\rho$. And my second is whether this is the best way to simulate this kind of data? Note that in the work that I'm doing, there is a causal relationship between X and T (X -> T)

$\endgroup$
5
$\begingroup$

To generate such a pair $(B,Y)$ with $B$ Bernoulli (with parameter $p$) and $Y$ normal, why not begin with a suitable binormal variable $(X,Y)$ and define $B$ to be the indicator that $X$ exceeds its $1-p$ quantile? By centering $(X,Y)$ at the origin and standardizing its marginals, the only question concerns what correlation $r$ should hold between $X$ and $Y$ so that the correlation between $B$ and $Y$ will be a given value $\rho$.

To this end, express $Y = r X + \sqrt{1-r^2}Z$ for independent standard Normal variables $X$ and $Z$. Set $x_0$ to be the $1-p$ quantile of $X$, so that $\Phi(x_0)=1-p$. (As is conventional, $\Phi$ is the standard Normal distribution and $\phi$ will be its density.)

Since the variance of $B$ is $p(1-p)$ and the variance of $Y$ is $1$, and $Y$ has zero mean, the correlation between $B$ and $Y$ is

$$\eqalign{ \rho&=\operatorname{Cor}(B,X) = \frac{E[BY] - E[B]E[Y]}{\sqrt{p(1-p)}\sqrt{1}}\\ &= \frac{E[B(rX+\sqrt{1-r^2}Z)]-0} {\sqrt{p(1-p)}} \\ &= \frac{rE[X\mid X \ge x_0]\Pr(X \ge x_0)}{\sqrt{p(1-p)}}. }$$

The conditional expectation is readily computed by integration, giving

$$E[X\mid X \ge x_0]\Pr(X \ge x_0) = \frac{1}{\sqrt{2\pi}}\int_{x_0}^\infty x e^{-x^2/2}dx = \frac{e^{-x_0^2/2}}{\sqrt{2\pi}} = \phi(x_0),$$

whence

$$\rho = \frac{r \phi(x_0)}{\sqrt{p(1-p)}}.$$

Solve this for $r$: by setting

$$r = \frac{\rho \sqrt{p(1-p)}}{\phi(x_0)},$$

$B$ and $Y$ will have correlation $\rho$.

Note that since it's necessary that $1-r^2\ge 0$, any values of $\rho$ that cause $|r|$ to exceed $1$ will not be achievable in this fashion. The figure plots feasible values of $r$ as a function of the desired correlation $\rho$ and Bernoulli parameter $p$: the contours range in increments of $1/10$ from $-1$ at the upper left through $+1$ at the upper right.

Figure

$\endgroup$
  • $\begingroup$ This is a great response and exactly what I'm looking for! Thanks! $\endgroup$ – Michael Barrowman Nov 29 '17 at 9:40
  • 1
    $\begingroup$ One thing I would ask for you to double check is whether the $p$ should be in the numerator in the final equation. I ran this through R and the correlation of the data was approximately $p * rho$ rather than $rho$ (e.g. when $p=0.5$, $rho$ = 0.5, cor(B,Y) = 0.25). When I removed $p$ from the numerator, it seems to be producing the correlation I want (although this may be a coincidence with errors in my code). $\endgroup$ – Michael Barrowman Nov 29 '17 at 9:50
  • 2
    $\begingroup$ @Michael Thank you for looking so closely and for your suggestion. You are correct that $p$ did not belong in the numerator. I had left out a factor of $p = \Pr(X\gt x_0)$ in the original derivation. That factor is now included and all subsequent formulas (and the figure) have been updated to reflect it. And--as I should originally have done--I checked the formula through some quick simulations, as you have been doing, and they indicate it's now correct. $\endgroup$ – whuber Nov 29 '17 at 15:41

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.