suppressPackageStartupMessages({
  library(readr)
  library(dplyr)
  library(tidyr)
  library(MatchIt)
  library(WeightIt)
  library(survey)
  library(survival)
})

out_dir <- "/Users/Zhuanz/Documents/ChatGPT/过度通气MIMICIV/outputs"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

dat <- read_csv(
  file.path(out_dir, "analysis_with_hv_definitions.csv"),
  show_col_types = FALSE
)

dat <- dat %>%
  mutate(
    death = as.integer(death),
    age = as.numeric(age),
    gcs = as.numeric(gcs),
    platelet = as.numeric(platelet),
    creatinine = as.numeric(creatinine),
    sbp = as.numeric(sbp),
    dbp = as.numeric(dbp),
    spo2 = as.numeric(spo2),
    hemoglobin = as.numeric(hemoglobin),
    invasive = as.integer(invasive),
    hv_vent_primary = as.integer(hv_primary & invasive == 1),
    hv_vent_strict = as.integer(hv_strict & invasive == 1)
  )

base_cov <- c("age", "gcs", "platelet", "creatinine")
full_cov <- c(
  "age", "gcs", "platelet", "creatinine",
  "sbp", "dbp", "spo2", "hemoglobin", "invasive"
)

or_ci <- function(fit, coef_name = "group") {
  or <- unname(exp(coef(fit)[coef_name]))
  ci <- unname(exp(confint(fit))[coef_name, ])
  c(OR = or, lower = ci[[1]], upper = ci[[2]])
}

run_definition <- function(data, groupvar, label = groupvar, restrict_vent = FALSE) {
  d <- data
  if (restrict_vent) {
    d <- d %>% filter(invasive == 1)
  }

  d <- d %>%
    filter(!is.na(.data[[groupvar]])) %>%
    mutate(group = as.integer(.data[[groupvar]])) %>%
    filter(
      complete.cases(
        select(., death, group, all_of(full_cov))
      )
    )

  if (sum(d$group) < 5) {
    return(NULL)
  }

  formula_base <- as.formula(paste("group ~", paste(base_cov, collapse = " + ")))

  set.seed(123)
  m <- matchit(
    formula_base,
    data = d,
    method = "nearest",
    ratio = 1,
    replace = FALSE,
    caliper = 0.2
  )
  matched <- match.data(m)

  fit_unmatched <- glm(
    as.formula(paste("death ~ group +", paste(full_cov, collapse = " + "))),
    data = d,
    family = binomial()
  )

  fit_matched <- glm(
    as.formula(paste("death ~ group +", paste(full_cov, collapse = " + "))),
    data = matched,
    family = binomial()
  )

  fit_conditional <- clogit(
    as.formula(paste(
      "death ~ group +", paste(base_cov, collapse = " + "),
      "+ strata(subclass)"
    )),
    data = matched
  )

  set.seed(123)
  w <- weightit(
    formula_base,
    data = d,
    method = "ps",
    estimand = "ATE"
  )
  d$w <- w$weights
  design <- svydesign(ids = ~1, weights = ~w, data = d)
  fit_iptw <- svyglm(
    as.formula(paste("death ~ group +", paste(base_cov, collapse = " + "))),
    design = design,
    family = binomial()
  )

  unmatched_ci <- or_ci(fit_unmatched)
  matched_ci <- or_ci(fit_matched)
  conditional_ci <- or_ci(fit_conditional)

  tibble(
    analysis = label,
    n_total = nrow(d),
    n_hv = sum(d$group),
    n_pairs = sum(matched$group),
    death_hv = mean(matched$death[matched$group == 1]),
    death_control = mean(matched$death[matched$group == 0]),
    unmatched_or = unmatched_ci[["OR"]],
    unmatched_lower = unmatched_ci[["lower"]],
    unmatched_upper = unmatched_ci[["upper"]],
    matched_or = matched_ci[["OR"]],
    matched_lower = matched_ci[["lower"]],
    matched_upper = matched_ci[["upper"]],
    conditional_or = conditional_ci[["OR"]],
    conditional_lower = conditional_ci[["lower"]],
    conditional_upper = conditional_ci[["upper"]],
    iptw_or = unname(exp(coef(fit_iptw)[["group"]])),
    iptw_lower = unname(exp(confint(fit_iptw))[["group", "2.5 %"]]),
    iptw_upper = unname(exp(confint(fit_iptw))[["group", "97.5 %"]])
  )
}

result_rows <- list(
  run_definition(dat, "hv_original", "Original pCO2-only"),
  run_definition(dat, "hv_primary", "Primary"),
  run_definition(dat, "hv_strict", "Strict"),
  run_definition(dat, "hv_very_strict", "Very strict"),
  run_definition(dat, "hv_vent_primary", "Ventilated: primary"),
  run_definition(dat, "hv_vent_strict", "Ventilated: strict"),
  run_definition(dat, "hv_primary", "Restricted to MV: primary", restrict_vent = TRUE),
  run_definition(dat, "hv_strict", "Restricted to MV: strict", restrict_vent = TRUE)
)

results <- bind_rows(result_rows)
write_csv(results, file.path(out_dir, "final_hv_results_v2.csv"))
print(results)