rm(list=ls())
set.seed(252)
n=1000
#Target population
exposure<-rnorm(n,mean=0,sd=0.3) #normal distribution for the exposure
effectmodifier<-as.numeric(runif(n,min=0,max=1)<0.25) #25% prevalence of the effect modifier
error<-rnorm(n,mean=0,sd=0.3)

#With causal effect
outcome<-1 + 2 * exposure + 3 * effectmodifier + 4 * exposure * effectmodifier + error

#Data frame
df<-as.data.frame(cbind(exposure,effectmodifier,outcome))

#output summary
nsim<-1000
vectorname<-c("Crude","lowCrudeCI","highCrudeCI","pvalueCrude","EM1","lowEM1CI","highEM1CI","pvalueEM1","EM0","lowEM0CI","highEM0CI","pvalueEM0")
simulationoutput<-as.data.frame(matrix(data=NA,nrow=nsim,ncol=length(vectorname)))
names(simulationoutput)<-vectorname

for (i in 1:nsim){
#sample population
population<-df[sample(c(1:n),100),]

#Crude linear regression
model<-glm(outcome~exposure,data=population)
simulationoutput$Crude[i]<-model$coefficients["exposure"][[1]]
simulationoutput$lowCrudeCI[i]<-confint(model)["exposure","2.5 %"]
simulationoutput$highCrudeCI[i]<-confint(model)["exposure","97.5 %"]
simulationoutput$pvalueCrude[i]<-summary(model)$coefficients["exposure","Pr(>|t|)"]

#Startified linear regression
populationEM1<-population[population$effectmodifier==1,]
modelEM1<-glm(outcome~exposure,data=populationEM1)
simulationoutput$EM1[i]<-modelEM1$coefficients["exposure"][[1]]
simulationoutput$lowEM1CI[i]<-confint(modelEM1)["exposure","2.5 %"]
simulationoutput$highEM1CI[i]<-confint(modelEM1)["exposure","97.5 %"]
simulationoutput$pvalueEM1[i]<-summary(modelEM1)$coefficients["exposure","Pr(>|t|)"]

populationEM0<-population[which(population$effectmodifier==0),]
modelEM0<-glm(outcome~exposure,data=populationEM0)
simulationoutput$EM0[i]<-modelEM0$coefficients["exposure"][[1]]
simulationoutput$lowEM0CI[i]<-confint(modelEM0)["exposure","2.5 %"]
simulationoutput$highEM0CI[i]<-confint(modelEM0)["exposure","97.5 %"]
simulationoutput$pvalueEM0[i]<-summary(modelEM0)$coefficients["exposure","Pr(>|t|)"]

  if(i%%100==0){
  print(i)
  }
}

#With bo causal effect:
outcome0<-1  + 3 * effectmodifier + error

#Data frame
df0<-as.data.frame(cbind(exposure,effectmodifier,outcome0))

#output summary
nsim<-1000
vectorname<-c("Crude","lowCrudeCI","highCrudeCI","pvalueCrude","EM1","lowEM1CI","highEM1CI","pvalueEM1","EM0","lowEM0CI","highEM0CI","pvalueEM0")
simulationoutput0<-as.data.frame(matrix(data=NA,nrow=nsim,ncol=length(vectorname)))
names(simulationoutput0)<-vectorname

for (i in 1:nsim){
  #sample population
  population<-df0[sample(c(1:n),100),]
  
  #Crude linear regression
  model<-glm(outcome0~exposure,data=population)
  simulationoutput0$Crude[i]<-model$coefficients["exposure"][[1]]
  simulationoutput0$lowCrudeCI[i]<-confint(model)["exposure","2.5 %"]
  simulationoutput0$highCrudeCI[i]<-confint(model)["exposure","97.5 %"]
  simulationoutput0$pvalueCrude[i]<-summary(model)$coefficients["exposure","Pr(>|t|)"]
  
  #Startified linear regression
  populationEM1<-population[population$effectmodifier==1,]
  modelEM1<-glm(outcome0~exposure,data=populationEM1)
  simulationoutput0$EM1[i]<-modelEM1$coefficients["exposure"][[1]]
  simulationoutput0$lowEM1CI[i]<-confint(modelEM1)["exposure","2.5 %"]
  simulationoutput0$highEM1CI[i]<-confint(modelEM1)["exposure","97.5 %"]
  simulationoutput0$pvalueEM1[i]<-summary(modelEM1)$coefficients["exposure","Pr(>|t|)"]
  
  populationEM0<-population[which(population$effectmodifier==0),]
  modelEM0<-glm(outcome0~exposure,data=populationEM0)
  simulationoutput0$EM0[i]<-modelEM0$coefficients["exposure"][[1]]
  simulationoutput0$lowEM0CI[i]<-confint(modelEM0)["exposure","2.5 %"]
  simulationoutput0$highEM0CI[i]<-confint(modelEM0)["exposure","97.5 %"]
  simulationoutput0$pvalueEM0[i]<-summary(modelEM0)$coefficients["exposure","Pr(>|t|)"]
  
  if(i%%100==0){
    print(i)
  }
}

hist(simulationoutput$Crude, breaks=15)
hist(simulationoutput$pvalueCrude, breaks=100)

hist(simulationoutput$EM1, breaks=15)
hist(simulationoutput$pvalueEM1, breaks=100)

hist(simulationoutput$EM0, breaks=15)
hist(simulationoutput$pvalueEM0, breaks=100)

hist(simulationoutput0$Crude, breaks=15)
hist(simulationoutput0$pvalueCrude, breaks=100)

hist(simulationoutput0$EM1, breaks=15)
hist(simulationoutput0$pvalueEM1, breaks=100)

hist(simulationoutput0$EM0, breaks=15)
hist(simulationoutput0$pvalueEM0, breaks=100)