library(MASS)

set.seed(43)

#Expit transform
expit<-function(x){return(1/(1+exp(-x)))}

n<-1000
p<-50

x.index<-1:p

#Forty two spurious
#Two instruments
#Three confounders
#Three predictors of outcome only 
al.select<-c(rep(0,42),1,1,1,1,1,0,0,0)
be.select<-c(rep(0,42),0,0,1,1,1,1,1,1)


#Confounders
confound.list<-x.index[al.select == 1 & be.select==1]

#Instruments
intrument.list<-x.index[al.select ==1 & be.select==0]

#Predictors of outcome only
outcomeonly.list<-x.index[al.select ==0  & be.select==1]

#Spurious variables
spurious.list<-x.index[al.select ==0  & be.select==0]

#Predictors of outcome
outcome.list<-x.index[be.select ==1]


#Predictors of exposure
treatment.list<-x.index[al.select ==1]

################################################

#Generate alpha parameters (for treatment allocation)
al<-rnorm(p+1,0,3.0)

#Generate beta parameters (for outcome regression)
be<-rnorm(p+1,0,3.0)

theta<-2.0
sig.Y<-1

p0<-42
p1<-2
p2<-3
p3<-3
Sig0<-rWishart(1,p0+1,diag(p0)/(p0+1))[,,1]
Sig1<-rWishart(1,p1+1,diag(p1)/(p1+1))[,,1]
Sig2<-rWishart(1,p2+1,diag(p2)/(p2+1))[,,1]
Sig3<-rWishart(1,p3+1,diag(p3)/(p3+1))[,,1]

Sig<-matrix(0,p,p)
Sig[1:p0,1:p0]<-Sig0
Sig[(p0+1):(p0+p1),(p0+1):(p0+p1)]<-Sig1
Sig[(p0+p1+1):(p0+p1+p2),(p0+p1+1):(p0+p1+p2)]<-Sig2
Sig[(p0+p1+p2+1):(p0+p1+p2+p3),(p0+p1+p2+1):(p0+p1+p2+p3)]<-Sig3