R help - Jul 2007 - generating correlated Bernoulli random variables

Hi all,
I was wondering how to generate samples for two RVs X1 and X2.

X1 ~ Bernoulli (p1)
X2 ~ Bernoulli (p2)

Also, X1 and X2 are correlated with correlation \rho.

Regards,
Vineet

	[[alternative HTML version deleted]]

On 3/07/2007, at 1:28 PM, Vineet Kumar wrote:
> Hi all,
> I was wondering how to generate samples for two RVs X1 and X2.
>
> X1 ~ Bernoulli (p1)
> X2 ~ Bernoulli (p2)
>
> Also, X1 and X2 are correlated with correlation \rho.
Back-of-envelope calculation, not checked, not tested:

	Let q1 = p2 + (rho/p1) * sqrt((1-p1)*(1-p2))

                q2 = (p2 - q1*p1)/(1-p1)

	Generate X1 with success probability p1.
	If X1 = 1, generate X2 with probability q1, else generate X2 with  
probability q2.

	Clearly there must be restrictions on the value of rho (in terms of  
p1 and p2) in
	order that the values of q1 and q2 lie between 0 and 1.

	I.e. for some values of rho your goal will be impossible to achieve.

				cheers,

					Rolf Turner

######################################################################
Attention:\ This e-mail message is privileged and confidenti...{{dropped}}

Typo in my previous email:

	sqrt((1-p1)*(1-p2)) should have read sqrt(p1*(1-p1)*p2*(1-p2))

I think!

			cheers,

				Rolf

######################################################################
Attention:\ This e-mail message is privileged and confidenti...{{dropped}}

>
> > Hi all,
> > I was wondering how to generate samples for two RVs X1 and X2.
> >
> > X1 ~ Bernoulli (p1)
> > X2 ~ Bernoulli (p2)
> >
> > Also, X1 and X2 are correlated with correlation \rho.
>   
You can use the rmvbin() function in the bindata package, e.g.,

require(bindata)
n=10
p1=0.5
p2=0.3
rho=0.2
rmvbin(n, c(p1,p2), bincorr=(1-rho)*diag(2)+rho)
?rmvbin

However, as pointed out before, rho is bounded below and above by some 
function of the marginal probabilities. (Try above code with rho=0.9)
You may want to use the odds ratio (which is unrestricted) to specify 
the association between the two binary variables and then convert this 
odds ratio, for given marginal probabilities p1 and p2, into a (valid) 
correlation rho to be used in rmvbin().

Here is some ad hoc code to do that:

bincorr <- function(OR, p1, p2) {    #from odds ratio to binary correlation
    if (OR==1) p11=p2-p2+p1*p2 else {
        p11_1=p2-(1/2/(1-OR)*(1-p1+OR*p1+p2-OR*p2-
              sqrt((-1+p1-OR*p1-p2+OR*p2)^2-4*(1-OR)*(p2-p1*p2))))
        p11_2=p2-(1/2/(1-OR)*(1-p1+OR*p1+p2-OR*p2-
              sqrt((-1+p1-OR*p1-p2+OR*p2)^2)-4*(1-OR)*(p2-p1*p2)))
        if (p11_1>0 && p11_1<=p1 && p11_1<p2) p11=p11_1
else p11=p11_2
    }
    bincorr=(p11-p1*p2)/sqrt(p1*(1-p1)*p2*(1-p2))
    return(bincorr)
}

For instance, try

sapply(c(0,0.5,1,1.5,3,10,100),function(x) bincorr(x,p1,p2))

to see the range of valid correlations for odds ratios between 0 and 
100, with p1 and p2 as above.

Bernhard Klingenberg
Dept. of Mathematics and Statistics
Williams College, MA
www.williams.edu/~bklingen