# SQL book
# cross sectional ---------------------------------------------------------
# 必要なライブラリをロード
library(tidyverse)
library(tableone)


# データ読み込み
# setwd("csvのあるフォルダのpath ")
# setwd("/Users/oonosachiko/Library/Mobile Documents/com~apple~CloudDocs/SQL本/統合原稿/取り出しテーブル")
df <- read_csv("crosssec.csv",na = c("", "NULL"), col_names = TRUE) #カラム名がなければcol_names = FALSE
# カラム名を設定
colnames(df) <- c("id", "malesex", "age", "SBP", "DBP", "diabetes", "drug_a", "drug_b") #すでにカラム名があれば不要

# Table 1 を作成（バイナリ変数は 1 のみ表示）
vars <- c("malesex", "age", "SBP", "DBP", "diabetes", "drug_a", "drug_b")
factorVars <- c("malesex", "diabetes", "drug_a", "drug_b")
table1 <- CreateTableOne(vars = vars, factorVars = factorVars, data = df)
# Table 1 をデータフレームに変換
table1_df <- print(table1, printToggle = FALSE) %>% as.data.frame()

# 行名の変更
rownames(table1_df) <- case_when(
  rownames(table1_df) == "n" ~ "Number",
  rownames(table1_df) == "malesex = 1 (%)" ~ "Male sex (%)",
  rownames(table1_df) == "age (mean (SD))" ~ "Age (years, mean (SD))",
  rownames(table1_df) == "SBP (mean (SD))" ~ "Systolic blood pressure (mmHg, mean (SD))",
  rownames(table1_df) == "DBP (mean (SD))" ~ "Diastolic blood pressure (mmHg, mean (SD))",
  rownames(table1_df) == "diabetes = 1 (%)" ~ "Diabetes (%)",
  rownames(table1_df) == "drug_a = 1 (%)" ~ "Drug A prescription (%)",
  rownames(table1_df) == "drug_b = 1 (%)" ~ "Drug B prescription (%)",
  TRUE ~ rownames(table1_df)
)

# CSVとして保存（行名を含める）
write.csv(table1_df, "table1.csv", row.names = TRUE)


# cohort ------------------------------------------------------------------
# 必要なライブラリをロード
library(tidyverse)
library(survival)
library(broom)
library(gt)

# データ読み込み
# setwd("csvのあるフォルダのpath ")
df <- read_csv("cohort.csv",na = c("", "NULL"), col_names = TRUE) #カラム名があればcol_names = FALSE
# カラム名を設定
colnames(df) <- c("id","HbA1c_first","SBP_first","DBP_first","malesex","age","MI_history","MI_outcome"
                  ,"lb_antiht","lb_antidiabetes","fwdays") #すでにカラム名があれば不要

# 除外基準
df <- df %>%
  filter(MI_history != 1, lb_antiht != 1, lb_antidiabetes != 1)
# **生存時間解析の準備**
df <- df %>%
  mutate(
    time = fwdays / 365.25,  # 日数を年単位に変換
    status = ifelse(MI_outcome == 1, 1, 0)) # 事象発生フラグ

# 発症率を1000人年あたりで算出
event_rate <- sum(df$status) / sum(df$time) * 1000
event_rate


# Cox回帰の実行
cox_model <- coxph(Surv(time, status) ~ age + malesex + HbA1c_first + SBP_first, data = df)
summary(cox_model)

# Cox回帰の結果を整理
cox_df <- tidy(cox_model, exponentiate = TRUE, conf.int = TRUE) %>%
  mutate(
    `ハザード比 [95%CI]` = sprintf("%.2f [%.2f - %.2f]", estimate, conf.low, conf.high),
    `p値` = format.pval(p.value, digits = 3, eps = 0.001)
  ) %>%
  select(変数 = term, `ハザード比 [95%CI]`, `p値`)


# `gt` で表示
gt(cox_df) %>%
  tab_header(title = "Cox回帰の結果") 


# new user design ---------------------------------------------------------

# 必要なライブラリをロード
library(tidyverse)
library(survival)
library(survminer)

# データ読み込み
# setwd("csvのあるフォルダのpath ")
df <- read_csv("newuser.csv",na = c("", "NULL"), col_names = TRUE) #カラム名があればcol_names = FALSE
# カラム名を設定
names(df)
colnames(df) <- c("id","malesex","MI_history","MI_outcome","MI_outcome_sa","lb_diabetes","fwdays","fwdays_sa"
                  ,"age","expB_conA") #すでにカラム名があれば不要

# 除外基準
df <- df %>%
  filter(MI_history != 1) %>%  # 既往のある人を除外
  filter(!is.na(expB_conA))    # 同月に両方処方されている人を除外

# 生存時間解析の準備
df <- df %>%
  mutate(
    time = fwdays / 365.25,  # 日数を年単位に変換
    status = ifelse(MI_outcome == 1, 1, 0)) # 事象発生フラグ

# 発症率（1000人年あたり）を計算
rate_by_group <- df %>%
  group_by(expB_conA) %>%
  summarise(
    events = sum(status, na.rm = TRUE),
    person_years = sum(time, na.rm = TRUE),
    incidence_rate = (sum(status, na.rm = TRUE) / sum(time, na.rm = TRUE)) * 1000
  )
rate_by_group

# 生存時間オブジェクトの作成
surv_obj <- Surv(time = df$time, event = df$status)

# カプランマイヤー法（Kaplan-Meier曲線）
fit_surv <- survfit(surv_obj ~ expB_conA, data = df)

ggsurvplot(
  fit_surv, data = df, risk.table = TRUE, pval = TRUE,
  conf.int = TRUE, xlab = "Follow-up Time (years)",
  ylab = "Survival Probability", legend.title = "Treatment Group"
)


# ログランク検定
survdiff(surv_obj ~ expB_conA, data = df)

# コックス比例ハザード回帰モデル
cox_model <- coxph(surv_obj ~ expB_conA + age + malesex + lb_diabetes, data = df)
summary(cox_model)


# Cox回帰の結果を整理
cox_df <- tidy(cox_model, exponentiate = TRUE, conf.int = TRUE) %>%
  mutate(
    `ハザード比 [95%CI]` = sprintf("%.2f [%.2f - %.2f]", estimate, conf.low, conf.high),
    `p値` = format.pval(p.value, digits = 3, eps = 0.001)
  ) %>%
  select(変数 = term, `ハザード比 [95%CI]`, `p値`)


# `gt` で表示
gt(cox_df) %>%
  tab_header(title = "Cox回帰の結果") 

# **感度解析（fwdays_saを用いた生存時間解析）**
df <- df %>%
  mutate(
    time_sa = fwdays_sa / 365.25,  # 日数を年単位に変換
    status_sa = ifelse(MI_outcome_sa == 1, 1, 0)  # 事象発生フラグ
  )

# 生存時間オブジェクトの作成
surv_obj_sa <- Surv(time = df$time_sa, event = df$status_sa)

# カプランマイヤー法（Kaplan-Meier曲線）
fit_surv_sa <- survfit(surv_obj_sa ~ expB_conA, data = df)

ggsurvplot(
  fit_surv_sa, data = df, risk.table = TRUE, pval = TRUE,
  conf.int = TRUE, xlab = "Follow-up Time (years)",
  ylab = "Survival Probability", legend.title = "Treatment Group"
)

survdiff(surv_obj_sa ~ expB_conA, data = df)

cox_model_sa <- coxph(surv_obj_sa ~ expB_conA + age + malesex + lb_diabetes, data = df)
summary(cox_model_sa)