Simulate phi-coefficient (correlation between dichotomous vars)

My situation is a bit more complicated and I'm not sure it is easily
solved.  The problem is that I must assume one of the vectors is
constant and generate one or more vectors that covary with the constant
vector.

One way is to sample from the 2x2 table with the specified means and
pearson correlation (phi):

for a fourfold table, a b
                      c d
with marginal proportions p1 and p2
cov <- phi * sqrt(p1*(1-p1)*p2*(1-p2))
a <- p1*p2 + cov
b <- p1*(1-p2) - cov
c <- (1-p1)*p2 - cov
d <- (1-p1)*(1-p2) + cov

Calculate the conditional probabilities from the above

P(X2=1|X1=1)= a/(a+b) = p2 + cov/p1
P(X2=1|X1=0)= c/(c+d) = p2 - cov/(1-p1)


condsim <- function(X, phi, p2, p1=NULL) {
  if (!all(X %in% c(0,1))) stop("expecting 1's and 0's")
  if (is.null(p1)) p1 <- mean(X)
  covar <- phi  * sqrt(p1*(1-p1)*p2*(1-p2))
  if (covar>0 && covar>(min(p1,p2)-p1*p2)) {
    warning("Specified correlation too large for given marginal proportions")
    covar <- min(p1,p2)-p1*p2
  }else if (covar<0 && covar < -min(p1*p2,(1-p1)*(1-p2))) {
    warning("Specified correlation too small for given marginal proportions")
    covar <- -min(p1*p2,(1-p1)*(1-p2))
  }
  Y <- X
  i1 <- X==1
  i0 <- X==0
  Y[i1] <- rbinom(sum(i1),1, p2 + covar/p1)
  Y[i0] <- rbinom(sum(i0),1, p2 - covar/(1-p1))
  data.frame(X,Y)
}

David Duffy.