install.packages("sandwich")
install.packages("lmtest")
install.packages("tidyverse")
install.packages("MatchIt")
install.packages("Rcpp")
install.packages("geepack")
install.packages("cobalt")
library("sandwich")
library("lmtest")
library("tidyverse")
library("MatchIt")
library("Rcpp")
library("geepack")
library("cobalt")

#データを作成する
set.seed(777)
N <- 10000
id<- (1:10000)
age <- floor(rnorm(N, 60, 10))
sex <- rbinom(N, 1, 0.5)
dl<-rbinom(N, 1, 1/(1+exp(-(-7+0.1*age-0.5*sex))))
CHDhistory<-rbinom(N, 1, 1/(1+exp(-(-4.81+0.05*age-0.5*sex+1*dl))))
aspirin<-rbinom(N, 1, 1/(1+exp(-(-5.6+0.05*age-0.5*sex+0.5*dl+2.87*CHDhistory))))
CHDevent<-rbinom(N, 1, 1/(1+exp(-(-10-1.4675*aspirin+0.1*age-0.5*sex+1*dl+2*CHDhistory))))
data<-data.frame(id, age, sex, dl, CHDhistory, aspirin, CHDevent)
summary(data)

#傾向スコアを求める
model_ps <- glm(aspirin ~ CHDhistory+age+sex+dl, family = binomial(link = "logit"), data = data)
data$ps <- predict(model_ps, data)

#1.回帰モデル（修正ポアソンモデル）
model1<- glm(CHDevent ~ aspirin + CHDhistory+age+sex+dl, family = poisson(link = "log"), data = data)
summary(model1)
estimate<-coeftest(model1, vcov = sandwich)[2,1]%>% exp()%>% round(3)
ci.low<-(coeftest(model1, vcov = sandwich)[2,1] + qnorm(0.05/2)*coeftest(model1, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
ci.up<-(coeftest(model1, vcov = sandwich)[2,1]+ qnorm(1-0.05/2)*coeftest(model1, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
paste("アスピリンのリスク比：点推定値 (95%信頼区間下限, 95%信頼区間上限)=", estimate, " (", ci.low,",", ci.up, ")", sep="")

#2. PSを回帰モデルの共変量として用いる
model2 <- glm(CHDevent ~ aspirin + ps, family = poisson(link = "log"), data = data)
summary(model2)
estimate<-coeftest(model2, vcov = sandwich)[2,1]%>% exp()%>% round(3)
ci.low<-(coeftest(model2, vcov = sandwich)[2,1] + qnorm(0.05/2)*coeftest(model2, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
ci.up<-(coeftest(model2, vcov = sandwich)[2,1]+ qnorm(1-0.05/2)*coeftest(model2, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
paste("アスピリンのリスク比：点推定値 (95%信頼区間下限, 95%信頼区間上限)=", estimate, " (", ci.low,",", ci.up, ")", sep="")

#3. PSマッチング
fit.psm <- matchit(aspirin ~ age+sex+dl+CHDhistory, method = "nearest", caliper =0.005, 
                   distance = "linear.logit",  data=data)
data.psm <- match.data(fit.psm)

#3-0. バランスチェック
summary(fit.psm)
bal.tab(fit.psm, thresholds = c(m = 0.1), un = TRUE)
bal.plot(fit.psm, var.name = "age")
bal.plot(fit.psm, var.name = "age", mirror = TRUE, type = "histogram")
love.plot(fit.psm, stats = c("mean.diffs"), 
          thresholds = c(m = 0.1), abs = TRUE, binary = "std", var.order = "unadjusted")

#3-1. マッチングされた集団での因果リスク比（マッチング考慮なし）
model3 <- glm(CHDevent ~ aspirin, family = poisson(link = "log"), data = data.psm)
summary(model3)
estimate<-coeftest(model3, vcov = sandwich)[2,1]%>% exp()%>% round(3)
ci.low<-(coeftest(model3, vcov = sandwich)[2,1] + qnorm(0.05/2)*coeftest(model3, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
ci.up<-(coeftest(model3, vcov = sandwich)[2,1]+ qnorm(1-0.05/2)*coeftest(model3, vcov = sandwich)[2,2]) %>% exp() %>% round(3)
paste("アスピリンのリスク比：点推定値 (95%信頼区間下限, 95%信頼区間上限)=", estimate, " (", ci.low,",", ci.up, ")", sep="")

#3-2. マッチングされた集団での因果リスク比（※推奨：マッチング考慮あり）
model4_1_results <-coeftest(model3, vcov. = vcovCL, cluster = ~subclass)
estimate<-model4_1_results[2,1]%>% exp()%>% round(3)
ci.low<-(model4_1_results[2,1] + qnorm(0.05/2)*model4_1_results[2,2]) %>% exp() %>% round(3)
ci.up<-(model4_1_results[2,1]+ qnorm(1-0.05/2)*model4_1_results[2,2]) %>% exp() %>% round(3)
paste("アスピリンのリスク比：点推定値 (95%信頼区間下限, 95%信頼区間上限)=", estimate, " (", ci.low,",", ci.up, ")", sep="")