Thanks, all for your help!
vineet


On 7/3/07, Bernhard Klingenberg <Bernhard.Klingenberg@williams.edu> 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.
> >
>
> 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
>
>

	[[alternative HTML version deleted]]

