gbd <-
function(n=1000,p=1,intercept=TRUE,betas=rep(0,p+intercept),mu=rep(0,p),Sigma=diag(1,p),seed=1){
require(MASS)
set.seed(seed)
X<-mvrnorm(n,mu,Sigma)
Xm<-X
if(intercept){Xm<-cbind(rep(1,n),X)}
lp<-Xm %*% betas
pv<-expit(lp)
Y<-rbinom(n,1,pv)
return(list(Y=Y,X=X))
}