pt <- readRDS("../data/dialysis_cohort.rds") %>%
  # REVIEW: check if the preemptive referrals should be included in the secondary analysis
  filter(is.na(referral_date) | referral_date >= FIRSTDIAL) %>%
  mutate(referred_within_1yr = if_else(referral_date - FIRSTDIAL + 1 <= 365.25, 1, 0, missing = 0),
         referred_to_tx_outside_state = ifelse(df_state != substr(CTR_CD, 1, 2), 1, 0),
         referred_within_1yr_outside_state = case_when(referred_within_1yr == 0 ~ 0,
                                                       referred_to_tx_outside_state == 1 ~ 1,
                                                       referred_to_tx_outside_state == 0 ~ 2),
         referred_within_1yr_outside_state = factor(referred_within_1yr_outside_state, 
                                                    levels = 0:2,
                                                    labels = c("No", "To Center Outside State", "To Center Inside State"))) %>%
  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_within_1yr_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_within_1yr_outside_state)) %>%
  filter(!is.na(insurance_medevid_first)) %>%
  mutate(median_household_income = median_household_income / 10000)

var_label <- list(
  "referred_within_1yr_outside_state" = "Referred within 1 year",
  "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_within_1yr_outside_state, insurance_medevid_first) %>%
  tbl_summary(by = referred_within_1yr_outside_state,
              type = list(all_continuous() ~ "continuous2",
                          all_dichotomous() ~ "categorical"),
              percent = "row",
              digits = list(all_categorical() ~ c(0, 1)),
              statistic = all_continuous() ~ c("{mean}±{sd}", "{median} ({p25}, {p75})"),
              label = var_label
              ) %>%
  modify_spanning_header(all_stat_cols() ~ "**Referred within 1 Year**")

create_rr_tbl <- function(m) {
  m2 <- m %>% 
    tbl_regression(exponentiate = T, 
                   label = var_label) %>% 
    bold_p() %>%
    modify_table_styling(columns = c(estimate),
                         rows = reference_row == TRUE, 
                         label = "**RRR**", 
                         missing_symbol = "Ref.") %>%
    modify_footnote(estimate ~ "RRR = Relative Risk Ratio") %>%
    remove_abbreviation()
  
  outcome1_tbl <- m2$table_body %>% 
    filter(groupname_col == "To Center Inside State") 
  outcome1 <- m2
  outcome1$table_body <- outcome1_tbl
  
  outcome2_tbl <- m2$table_body %>%
    filter(groupname_col == "To Center Outside State") 
  outcome2 <- m2
  outcome2$table_body <- outcome2_tbl
  
  m2 <- tbl_merge(list(outcome1, outcome2),
                  tab_spanner = c("**To Center Inside State**", "**To Center Outside State**")) 
  
  m2$table_styling$header <- m2$table_styling$header %>%
    mutate(hide = ifelse(column %in% c("groupname_col_1", "groupname_col_2"), T, hide))

  return(m2)
}

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

m_crude <- multinom(
  referred_within_1yr_outside_state ~ adjacent_to_boundary + insurance_medevid_first, 
  data = desc_pt,
  trace = F
  )
m_adjusted <- multinom(
  as.formula(paste("referred_within_1yr_outside_state ~ adjacent_to_boundary + insurance_medevid_first + ", covars)), 
  data = desc_pt,
  trace = F
  ) 

m_interaction_crude <- multinom(
  referred_within_1yr_outside_state ~ adjacent_to_boundary * insurance_medevid_first, 
  data = desc_pt,
  trace = F
  )
m_interaction_adjusted <- multinom(
  as.formula(paste("referred_within_1yr_outside_state ~ adjacent_to_boundary * insurance_medevid_first + ", covars)), 
  data = desc_pt,
  trace = F
  )