rm(list = ls(all.names = TRUE))
install.packages("geepack")
install.packages("broom")
install.packages("devtools")
devtools::install_github("BS1125/CMAverse")
#install.packages("remotes") #if the above code doesn't  work to install CMAverse
#remotes::install_github("LindaValeri/CMAverse", force = TRUE)

library("geepack")
library("broom")
library("CMAverse")

#データを作成する
set.seed(777)
N <- 100000
ID<- (1:N)
C1 <-rbinom(N,1,0.3)  #E-M confounder
C2 <-rbinom(N,1,0.3)  #M-Y confounder
C3 <-rbinom(N,1,0.3)  #E-Y confounder
E <-rbinom(N,1,1/(1+exp(-(log(0.30/0.70) + log(2.0)*C1+ log(2.0)*C3))))
M <-rbinom(N,1,1/(1+exp(-(log(0.20/0.80) + log(3.0)*E + log(2.0)*C1 + log(2.0)*C2))))
Y <-rbinom(N,1,1/(1+exp(-(log(0.03/0.97) + log(1.5)*E + log(1.5)*M 
                          + log(1.2)*M*E  + log(3.0)*C2+ log(3.0)*C3)))) 
d<-as.data.frame(cbind(C1,C2,C3,E,M,Y,ID))
summary(d)

## 回帰モデル
# Total effect
te_out<-glm(data=d,Y ~ E + C1+C3,family = binomial)
summary(te_out)
exp(cbind(OR=coef(te_out),confint(te_out)))

# Controlled direct effect
cde_out<-glm(data=d,Y ~ E + C1 + C2 +C3  + M + E*M,family = binomial)
summary(cde_out)
exp(cbind(OR=coef(cde_out),confint(cde_out)))

#Total effect of E: OR = 1.83872059  
#Controlled direct effect of E: OR =1.47904848

## CMAverseパッケージ
# https://bs1125.github.io/CMAverse/articles/overview.html#the-marginal-structural-model-1
# https://bs1125.github.io/CMAverse/articles/single_mediator.html#the-marginal-structural-model-1
set.seed(123)
fit_msms <- cmest(data = d, model = "msm", outcome = "Y", exposure = "E", #Marginal Structure approach
                  mediator = "M", basec = c("C1", "C2", "C3"), EMint = TRUE,
                  ereg =  glm(E ~ C1+C3, family = binomial(link="logit"), data = d), 
                  yreg = glm(Y ~ E*M, family = binomial(link="logit"), data = d), 
                  mreg = list(glm(M ~ E, family = binomial(link="logit"), data = d)),
                  
                  wmnomreg = list(glm(M ~ E, family = binomial(link="logit"), data = d)), 
                  
                  wmdenomreg = list(glm(M ~ E+C1+C2, family = binomial(link="logit"), data = d)),
                  astar = 0, a = 1, mval = list(0), 
                  estimation = "imputation", inference = "bootstrap", nboot = 2)
summary(fit_msms)
#Total effect of E: RR = 1.8233943 
#Controlled direct effect of E: RR = 1.4421112
#See the summary output for other causal effects of interest

## IPWアプローチ
#IPW for exposure
fitipwe<-glm(data=d,E ~ C1+C3,family = binomial)
d$p_e <- predict(fitipwe, d, type="response")
fitipwen<-glm(data=d,E ~ 1,family = binomial)
d$p_en <- predict(fitipwen, d, type="response")
d$sw_e<-ifelse(E==1,yes= d$p_en / d$p_e,
               no= ifelse(E==0, yes= (1-d$p_en) / (1-d$p_e), no=NA))

#IPW for mediator
fitipwm<-glm(data=d,M ~ E+C1+C2+C3,family = binomial)
d$p_m <- predict(fitipwm, d, type="response")
fitipwmn<-glm(data=d,M ~ 1,family = binomial)
d$p_mn <- predict(fitipwmn, d, type="response")
d$sw_m<-ifelse(M==1,yes=d$p_mn / d$p_m,
               no= ifelse(M==0, yes= (1-d$p_mn) / (1-d$p_m),no=NA))
# Combine IPWs
d$sw_e_m<- d$sw_e*d$sw_m   

# Total effect of E
msm.te <- geeglm(
  Y ~ E,
  data = d,
  weights = sw_e,
  id = ID,
  family=binomial(link=logit), 
  corstr = "independence"
)
summary(msm.te)
msm.te.out<-broom::tidy(x = msm.te, conf.int = TRUE)
exp(cbind(RR=msm.te.out[2,2],msm.te.out[2,6:7]))

# Controlled direct effect
msm.cde <- geeglm(
  Y ~ E*M,
  data = d,
  weights = sw_e_m,
  id = ID,
  family=binomial(link=logit),
  corstr = "independence"
)
summary(msm.cde)
msm.cde.out<-broom::tidy(x = msm.cde, conf.int = TRUE)
exp(cbind(RR=msm.cde.out[2,2],msm.cde.out[2,6:7]))
#Total effect of E: RR = 1.81 
#Controlled direct effect of E: RR = 1.44