yn_factor <- function(x) factor(x, levels = c(0, 1), labels = c("No", "Yes")) 

pt <- readRDS("../data/dialysis_cohort.rds") %>%
  filter(referral_date >= FIRSTDIAL) %>%
  mutate(referred_to_tx_outside_state = factor(ifelse(df_state != substr(CTR_CD, 1, 2), "Yes", "No"), levels = c("No", "Yes"))) %>%
  mutate_at(c("nearest_tx_outside_state", "adjacent_to_boundary", "nearest_state_boundary_within_10m", "tx_outside_state_driving_time"), yn_factor) %>%
  rowwise() %>%
  mutate(como_na = sum(is.na(COMO_CHF), is.na(COMO_ASHD), is.na(COMO_OTHCARD), is.na(COMO_CVATIA), is.na(COMO_PVD), is.na(COMO_HTN), is.na(COMO_DIABETES), is.na(COMO_COPD), is.na(COMO_CANC), is.na(COMO_TOBAC)),
         como_cnt = sum(COMO_CHF == "Y", COMO_ASHD == "Y", COMO_OTHCARD == "Y", COMO_CVATIA == "Y", COMO_PVD == "Y", COMO_HTN == "Y", COMO_DIABETES == "Y", COMO_COPD == "Y", COMO_CANC == "Y", COMO_TOBAC == "Y", na.rm = T),
         como_cnt = ifelse(como_na == 10, NA, como_cnt),
         ADI_NATRANK = ifelse(ADI_NATRANK %in% c("PH", "GQ", "PH-GQ", "QDI"), NA, as.numeric(ADI_NATRANK))) %>% 
  select(referred_to_tx_outside_state, insurance_medevid_first, adjacent_to_boundary, male, INC_AGE, RACE, HISPANIC, EMPCUR, NEPHCARE, TRCERT, COMO_INST, como_cnt, inambulate, median_household_income, public_transportation_pct, ADI_NATRANK) 

desc_pt <- pt %>% 
  filter(!is.na(referred_to_tx_outside_state)) %>%
  filter(!is.na(insurance_medevid_first)) %>%
  mutate(median_household_income = median_household_income / 10000)

var_label <- list(
  "referred_to_tx_outside_state" = "Referred to transplant outside state",
  "adjacent_to_boundary" = "Located in the county adjacent to state boundary",
  "insurance_medevid_first" = "Insurance type",
  "male" = "Sex",
  "INC_AGE" = "Age at Dialysis Start",
  "RACE" = "Race", 
  "HISPANIC" = "Ethnicity", 
  "EMPCUR" = "Employment Status", 
  "NEPHCARE" = "Prior Nephrology Care", 
  "TRCERT" = "Patient Completing Home Dialysis Training", 
  "COMO_INST" = "Institutionalized", 
  "como_cnt" = "Number of Comorbidities",
  "inambulate" = "Inability to Ambulate",
  "median_household_income" = "Median Household Income (×$10k)",
  "public_transportation_pct" = "Public transportation to Work (%)", 
  "ADI_NATRANK" = "ADI National Rank"
  )

desc_pt %>%
  tbl_summary(by = adjacent_to_boundary,
              type = list(all_continuous() ~ "continuous2",
                          all_dichotomous() ~ "categorical"),
              statistic = all_continuous() ~ c("{mean}±{sd}", "{median} ({p25}, {p75})"),
              label = var_label
              ) %>%
  modify_spanning_header(all_stat_cols() ~ "**Located in the county adjacent to state boundary**") %>%
  add_overall()

desc_pt %>%
  select(referred_to_tx_outside_state, insurance_medevid_first) %>%
  tbl_summary(by = referred_to_tx_outside_state,
              type = list(all_continuous() ~ "continuous2",
                          all_dichotomous() ~ "categorical"),
              percent = "row",
              digits = list(all_continuous() ~ 1),
              statistic = all_continuous() ~ c("{mean}±{sd}", "{median} ({p25}, {p75})"),
              label = var_label
              ) %>%
  modify_spanning_header(all_stat_cols() ~ "**Referred to transplant outside state**") 

or_tbl <- function(m) {
  m2 <- m %>% 
    tbl_regression(exponentiate = T, 
                   label = var_label) %>% 
    bold_p() %>%
    modify_table_styling(columns = c(estimate),
                         rows = reference_row == TRUE, 
                         missing_symbol = "Ref.")
  return(m2)
}

clean_tbl <- function(tbl, varname) {
  varname2 <- str_to_title(str_replace_all(varname, "_", " "))
  tbl_clean <- tbl %>% 
    rename_with(~ "At", all_of(varname)) %>%
    mutate_at(vars(AME, lower, upper, p), ~ round(., 5)) %>%
    mutate(Insurance = paste0(substr(factor, nchar("insurance_medevid_first") + 1, nchar(factor)), " vs. TM w/o MDCD"), 
           Insurance = factor(Insurance, levels = paste0(c("TM w/ MDCD", "MA w/o MDCD", "MA w/ MDCD", "Employer", "Other", "None"), " vs. TM w/o MDCD")),
           At = factor(At, levels = c("No", "Yes"), labels = c(paste0(varname2, ": No"), paste0(varname2, ": Yes"))),
           CI = paste0("(", lower, ", ", upper, ")")) %>%
    select(At, Insurance, AME, CI, p) %>%
    arrange(Insurance)
  return(tbl_clean)
}

covars <- "ADI_NATRANK + male + INC_AGE + RACE + HISPANIC + EMPCUR + NEPHCARE + TRCERT + COMO_INST + como_cnt + inambulate + median_household_income + public_transportation_pct"

# whether the dialysis facility is located in the county adjacent to state boundary
m1_c <- glm(referred_to_tx_outside_state ~ adjacent_to_boundary + insurance_medevid_first, 
            data = desc_pt, family = binomial) 
m1_a <- glm(as.formula(paste("referred_to_tx_outside_state ~ adjacent_to_boundary + insurance_medevid_first + ", covars)), 
            data = desc_pt, family = binomial) 
m1_in_c <- glm(referred_to_tx_outside_state ~ adjacent_to_boundary * insurance_medevid_first, 
               data = desc_pt, family = binomial)
m1_in_a <- glm(as.formula(paste("referred_to_tx_outside_state ~ adjacent_to_boundary * insurance_medevid_first + ", covars)), , 
               data = desc_pt, family = binomial)