rm(list=ls())
set.seed(252)
n=1000
#Target population
U<-runif(n,min=0,max=1)
errorX<-rnorm(n,mean=0,sd=0.3) 
errorY<-rnorm(n,mean=0,sd=0.3)

#With causal effect
X<-3 + 2 * U + errorX
S2<-as.numeric(U<0.75)
Y<-1 + 2 * X + 3 * U  + errorY
continuousS1 <- X - Y + 2 * X * Y
sortedcontinuousS1<-sort(continuousS1)
S1<-as.numeric(continuousS1<=sortedcontinuousS1[1000])

#Data frame
df<-as.data.frame(cbind(U,X,Y,S2,S1,errorX,errorY))

#Effect of X on Y:
samplepopulation<-df[sample(c(1:n),500),]
population<-samplepopulation[samplepopulation$S2==1,]

#Crude linear regression
model<-glm(Y~X,data=population)
round(cbind(coef=model$coefficients["X"][[1]],lower=confint(model)["X","2.5 %"],upper=confint(model)["X","97.5 %"],pvalue=summary(model)$coefficients["X","Pr(>|t|)"]),2)


#Assuming no true causal effect:
#Target population
U<-runif(n,min=0,max=1)
errorX<-rnorm(n,mean=0,sd=0.3) 
errorY<-rnorm(n,mean=0,sd=0.3)

#With causal effect
X<-3 + 2 * U + errorX
S2<-as.numeric(U<0.75)
Y<-1 + 3 * U  + errorY
continuousS1 <- X - Y + 2 * X * Y
sortedcontinuousS1<-sort(continuousS1)
S1<-as.numeric(continuousS1<=sortedcontinuousS1[1000])

#Data frame
df<-as.data.frame(cbind(U,X,Y,S2,S1,errorX,errorY))

#Effect of X on Y:

#output summary
nsim<-1000
vectorname<-c("coef","pvalue")
simulationoutput<-as.data.frame(matrix(data=NA,nrow=nsim,ncol=length(vectorname)))
names(simulationoutput)<-vectorname

for (i in 1:nsim){
  #sample population
  samplepopulation<-df[sample(c(1:n),500),]
  population<-samplepopulation[samplepopulation$S2==1,]
  
  #Crude linear regression
  model<-glm(Y~X,data=population)
  simulationoutput$coef[i]<-model$coefficients["X"][[1]]
  simulationoutput$pvalue[i]<-summary(model)$coefficients["X","Pr(>|t|)"]
  
  if(i%%100==0){
    print(i)
  }
}
a<-sort(simulationoutput$coef)
alpha<-100*sum(as.numeric(simulationoutput$pvalue<0.05))/length(simulationoutput$pvalue)
round(cbind(coef=mean(a),lower=a[25],upper=a[975],type1errorPercent=alpha),2)