##
##更新日：2021/12/27
##

##使用するパッケージのインストールとロード
install.packages("tableone")
library(tableone)
install.packages("tidyverse")
library("tidyverse")

##データの作成
df1 <- data.frame(E = 1, D=1, Z=1, ID = 1:60)
df2 <- data.frame(E = 1, D=0, Z=1, ID = 61:100)
df3 <- data.frame(E = 0, D=1, Z=1, ID = 101:700)
df4 <- data.frame(E = 0, D=0, Z=1, ID = 701:1100)
df5 <- data.frame(E = 1, D=1, Z=0, ID = 1101:1140)
df6 <- data.frame(E = 1, D=0, Z=0, ID = 1141:2000)
df7 <- data.frame(E = 0, D=1, Z=0, ID = 2001:6400)
df8 <- data.frame(E = 0, D=0, Z=0, ID = 6401:101000)
df <- rbind(df1,df2,df3,df4,df5,df6,df7,df8)

#作成したデータの確認
summary(df)

#本文中の表と値が一致することを確認
CreateTableOne(vars = c("Z","D"), strata = c("E"), data = df,test=FALSE)
count <-addmargins(table(df$D,df$E))
print(count)

#バイアス調整前のRR
RR<-(count[2,2]/count[3,2])/(count[2,1]/count[3,1]) %>% round(3)
print(paste("バイアス調整前RR=",RR))
#exp(glm(D~E, family=binomial(link = "log"),data=df)$coef) #glmを用いた分析

#穴埋めしたシナリオにおけるバイアス調整後のリスク比(RR)
AdjustedRR <-exp(glm(D~E+Z, family=binomial(link = "log"),data=df)$coef[2]) %>% round(3)# 共変量調整した結果
print(paste("バイアス調整後RR(共変量調整)=",AdjustedRR))
AdjustedRR1 <- exp(glm(D~E, family=binomial(link = "log"),data=df,subset=df$Z==1)$coef[2]) %>% round(3)　#Z=1に限定
print(paste("バイアス調整後RR(Z=1に限定)=",AdjustedRR1))
AdjustedRR0 <- exp(glm(D~E, family=binomial(link = "log"),data=df,subset=df$Z==0)$coef[2])%>% round(3)　#Z=0に限定
print(paste("バイアス調整後RR(Z=0に限定)=",AdjustedRR0))

#Bias formulaを用いてシナリオを与えたときのバイアス分析
rrdz <- 13.5
pz0 <- 0.01
pz1 <-0.1
bias <- ((rrdz*pz1 + 1 - pz1) / (rrdz*pz0 + 1 - pz0)) # Bias formula
Adjusted2RR <- RR/bias %>% round(3)
print(paste("バイアス調整後RR=",Adjusted2RR))

#様々なシナリオを与えた場合のバイアス分析
rrdz <- 1:20
pz0 <- 0.01
pz1 <-1:30*0.01
bias <- ((rrdz*pz1 + 1 - pz1) / (rrdz*pz0 + 1 - pz0)) # Bias formula
Adjusted3RR <- (RR/bias) %>% round(3)
print(paste("RRdz=",rrdz,"pz0=",pz0,"pz1=",pz1,":",Adjusted3RR))