rm(list = ls(all.names = TRUE))
install.packages("lmtest")
install.packages("sandwich")
install.packages("tidyverse")
install.packages("boot")
library("lmtest")
library("sandwich")
library("tidyverse")
library(boot)


#データを作成する
set.seed(777)
N <- 100000
id<- (1:100000)
age <- floor(rnorm(N, 60, 10))
sex <- rbinom(N, 1, 0.5)
ldl_base<-rnorm(N, 120+0.1*age-1.5*sex, 15)
CHDhistory<-rbinom(N, 1, 1/(1+exp(-(-6+0.05*age-0.5*sex+0.01*ldl_base))))
statin_base<-rbinom(N, 1, 1/(1+exp(-(-5.6+0.05*age-0.5*sex+0.01*ldl_base+1*CHDhistory))))
ldl_follow<-rnorm(N, 120+0.1*age-30*sex+0.05*ldl_base+10*CHDhistory-20*statin_base, 15)
statin_follow<-rbinom(N, 1, 1/(1+exp(-(-5+0.05*age-0.5*sex+2*statin_base+0.01*(ldl_base-100)+0.03*(ldl_follow-100)+0.2*CHDhistory))))
CHDevent<-rbinom(N, 1, 1/(1+exp(-(-12-0.2*statin_base-0.2*statin_follow-0.05*statin_base*statin_follow+0.1*age-0.5*sex+0.01*ldl_base+0.03*ldl_follow+1*CHDhistory))))
data<-data.frame(id, age, sex, ldl_base, ldl_follow, CHDhistory, statin_base, statin_follow, CHDevent)
summary(data)


#IPTW
iptw<-function(d0, indices) {
  d <- d0[indices,] 
  #T1に対する重みを求める
  psmodel1<-glm(statin_base ~ age+sex+ldl_base+CHDhistory, 
                data=d, family=binomial(link="logit"))
  d$ps1<-predict(psmodel1, type = "response")
  psmodel2<-glm(statin_base ~ 1, data=d, family=binomial(link="logit"))
  d$ps2<-predict(psmodel2, type = "response")
  d$weight1 <- ifelse(d$statin_base==1, d$ps2/d$ps1, (1-d$ps2)/(1-d$ps1))
  #d$weight1[d$weight1>quantile(d$weight1, c(.99))]<- quantile(d$weight1, c(.99))  #必要に応じてtruncateする
  #d$weight1[d$weight1<quantile(d$weight1, c(.01))]<- quantile(d$weight1, c(.01))  #必要に応じてtruncateする
  
  #T2に対する重みを求める
  psmodel3<-glm(statin_follow ~ age+sex+ldl_base+CHDhistory+ldl_follow+statin_base, 
                data=d, family=binomial(link="logit"))
  d$ps3<-predict(psmodel3, type = "response")
  psmodel4<-glm(statin_follow ~ 1, data=d, family=binomial(link="logit"))
  d$ps4<-predict(psmodel4, type = "response")
  d$weight2 <- ifelse(d$statin_follow==1, d$ps4/d$ps3, (1-d$ps4)/(1-d$ps3))
  #d$weight2[d$weight2>quantile(d$weight2, c(.99))]<- quantile(d$weight2, c(.99))  #必要に応じてtruncateする
  #d$weight2[d$weight2<quantile(d$weight2, c(.01))]<- quantile(d$weight2, c(.01))  #必要に応じてtruncateする
  
  #T1, T2に対する重みを求める
  weight<-d$weight1*d$weight2
  #hist(weight)
  #summary(weight)
  #(重みをApplyしたときに曝露群・非曝露群で共変量のバランスが取れているかを確認する) 
    model<-glm(CHDevent ~ statin_base*statin_follow, data=d, family = poisson(link = "log"),
               weights=weight)
    
    #Return value
    estimate_bt2 = exp(sum(coef(model)[c('statin_follow')]))
    estimate_bt1 = exp(sum(coef(model)[c('statin_base')]))
    estimate_bt12 = exp(sum(coef(model)[c('statin_base','statin_follow', 'statin_base:statin_follow')]))
    return(c(bt2=estimate_bt2, bt1=estimate_bt1, bt12=estimate_bt12))
}

#ブートストラップ法で95% CIを計算する。
results_iptw <- boot(data=data, statistic=iptw, R=200)
mean_iptw <- results_iptw$t0
se_iptw <- c(sd(results_iptw$t[,1]), sd(results_iptw$t[,2]), sd(results_iptw$t[,3]))
ll_iptw <- mean_iptw - qnorm(0.975)*se_iptw
ul_iptw <- mean_iptw + qnorm(0.975)*se_iptw
bootstrap_iptw <- data.frame(cbind(c("Statin-/+", "Statin+/-", 
                                     "Statin+/+"), mean_iptw, ll_iptw, ul_iptw))
bootstrap_iptw


#G-computations
gcomp =function(d0, indices) {
  d <- d0[indices,] 
  ##Step 1. M (LDL follow) と Y (CHDevent)の予測モデルを構築する
  d$statin_temp1<-d$statin_base
  d$statin_temp2<-d$statin_follow
  reg_M<-glm(ldl_follow ~ age+sex+ldl_base+CHDhistory+statin_temp1, data=d, family=gaussian)
  d$ldl_temp<-d$ldl_follow
  reg_Y<-glm(CHDevent ~ age+sex+ldl_base+CHDhistory+ldl_temp+statin_temp1*statin_temp2, data=d, family=binomial(link="logit"))
  
  #Step2. データをコピーして、スタチン内服の状態を新しく割り当てる
  data1<-data2<-data3<-data4<-d 
  #全員スタチン内服（研究開始時）している/全員スタチン内服（１年後フォロー）している場合
  data1$statin_temp1=1 
  data1$statin_temp2=1 
  #全員スタチン内服（研究開始時）している/全員スタチン内服（１年後フォロー）していない場合
  data2$statin_temp1=1 
  data2$statin_temp2=0 
  #全員スタチン内服（研究開始時）していない/全員スタチン内服（１年後フォロー）している場合
  data3$statin_temp1=0 
  data3$statin_temp2=1 
  #全員スタチン内服（研究開始時）していない/全員スタチン内服（１年後フォロー）していない場合
  data4$statin_temp1=0 
  data4$statin_temp2=0 
  #４つのコピーデータを結合する。
  newdata<-rbind(data1, data2, data3, data4)
  
  #Step3: 結合したデータで、新しく割り当てられたスタチン内服の状態のもと、ＭとＹを改めて予測する
  N<-nrow(d)
  newdata$ldl_temp<-predict(reg_M,newdata=newdata)
  newdata$pY<-rbinom(4*N, size = 1, prob=1/(1+exp(-(predict(reg_Y,newdata=newdata)))))
  
  #Step4: ステップ3で算出されたアウトカムの予測値に対して新しく割り当てられたスタチン内服の状態を回帰モデルに当てはめて因果効果を推定する
  reg_po<-glm(pY ~ statin_temp1*statin_temp2, data=newdata, family = poisson(link = "log"))
  
  #Step5: Return value
  reg_po_results <-coeftest(reg_po, vcov. = vcovCL, cluster = ~id)
  estimate_bt2 = exp(reg_po_results['statin_temp2', 1])
  estimate_bt1 = exp(reg_po_results['statin_temp1', 1])
  estimate_bt12 = exp(reg_po_results['statin_temp1', 1]+reg_po_results['statin_temp2', 1]+ reg_po_results['statin_temp1:statin_temp2', 1])
  return(c(bt2=estimate_bt2, bt1=estimate_bt1, bt12=estimate_bt12))
}

#ブートストラップ法で95% CIを計算する。
results_gcomp <- boot(data=data, statistic=gcomp, R=200)
mean_gcomp <- results_gcomp$t0
se_gcomp <- c(sd(results_gcomp$t[,1]), sd(results_gcomp$t[,2]), sd(results_gcomp$t[,3]))
ll_gcomp <- mean_gcomp - qnorm(0.975)*se_gcomp
ul_gcomp <- mean_gcomp + qnorm(0.975)*se_gcomp
bootstrap_gcomp <- data.frame(cbind(c("Statin-/+", "Statin+/-", 
                                      "Statin+/+"), mean_gcomp, ll_gcomp, ul_gcomp))
bootstrap_gcomp