Generate Propensity Score Matches

This file was compiled on 2022-05-27 16:34:30 by jsimmeri on argon-lc-g14-35.hpc.

library(tidyverse)
library(survival)

Data Processing

Generate Outcomes

outcomes <- read_rds("/Shared/lss_jsimmeri_backup/data/pdd/bph_rx/model_data.rds") %>%
  select(enrolid, drug, index_date, event_date, end_date) %>%
  mutate(
    event = ifelse(is.na(event_date), FALSE, event_date <= end_date),
    survival_time = ifelse(event, event_date, end_date) - index_date
  ) %>%
  select(enrolid, drug, event, survival_time) %>%
  distinct() %>%
  mutate(enrolid = as.numeric(enrolid))

Load Matching

tz_tam <- read_rds("/Shared/lss_jsimmeri_backup/data/pdd/matches/tz_tam.rds")
tz_5ari <- read_rds("/Shared/lss_jsimmeri_backup/data/pdd/matches/tz_5ari.rds")
tam_5ari <- read_rds("/Shared/lss_jsimmeri_backup/data/pdd/matches/tam_5ari.rds")

Fit Model

cox_fit <- function(outcomes, matches, treatment, reference) {
  if (is.null(matches)) {
    data <- outcomes %>%
      filter(drug %in% c(treatment, reference)) %>%
      mutate(pair_id = row_number()) %>%
      mutate(treatment = as.numeric(drug == treatment))
    type <- "unmatched"
  } else {
    data <- outcomes %>%
      inner_join(matches %>%
                   mutate(enrolid = as.numeric(enrolid)), 
                 by = "enrolid")
    type <- "matched"
  }

  survfit(Surv(survival_time, event) ~ treatment, 
          data = data) %>%
    broom::tidy() %>%
    mutate(
      label = ifelse(strata == "treatment=1", treatment, reference)
    ) %>%
    mutate(
      across(one_of("estimate", "conf.low", "conf.high"), 
             function(x) return(100 * x))
    ) %>%
    ggplot(aes(x = time / 365, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_step(aes(color = label)) + 
    geom_ribbon(aes(fill = label), alpha = 0.33) +
    theme_bw() +
    theme(legend.position = c(0.25, 0.25), 
          legend.background = element_rect(fill = NA)) +
    labs(x = "Years of Follow-Up", y = "Percent Remaining Free of Dementia",
         color = "", fill = "") 

  if (treatment == "tz/dz/az") {
    treatment <- "tz"
  }

  ggsave(glue::glue("~/projects/pdd/analysis/km_{type}_{treatment}_{reference}.png"),
         width = 16/3, height = 9/3)

  fit <- coxph(Surv(survival_time, event) ~ treatment,
          data = data, 
          cluster = pair_id, robust = TRUE) %>%
    broom::tidy() %>%
    mutate(
      treatment = treatment, 
      reference = reference,
      n_pairs = NROW(unique(data$pair_id)),
      hr = exp(estimate),
      lb = exp(estimate - 1.96 * robust.se),
      ub = exp(estimate + 1.96 * robust.se)
    ) %>%
    select(
      treatment, reference, n_pairs, hr, lb, ub, p.value
    )

  schr <- cox.zph(coxph(Surv(survival_time, event) ~ treatment, 
                        data = data,  
                        cluster = pair_id, robust = TRUE))
  schr_values <- tibble(
      x = schr$x,
      y = as.numeric(schr$y[, 1])
    ) %>%
    with(cor.test(x, y)) %>%
    broom::tidy()

  schr_plot <- tibble(
    x = schr$time / 365,
    y = as.numeric(schr$y[, 1])
  ) %>%
    ggplot(aes(x = x, y = y)) + 
    geom_point(alpha = 0.25) + 
    geom_smooth(method = "loess") + 
    theme_bw() + 
    labs(x = "Follow Up Time", y = "Schoenfeld Residual",
         subtitle = glue::glue("r = {round(schr_values$estimate, 2)}, p = {round(schr_values$p.value, 4)}"))

  ggsave(glue::glue("~/projects/pdd/analysis/schoenfeld_{type}_{treatment}_{reference}.png"),
         width = 16/3, height = 9/3)
  return(fit)
}

cox_fit(outcomes, NULL, "tz/dz/az", "tamsulosin")

# A tibble: 1 x 7
  treatment reference  n_pairs    hr    lb    ub  p.value
  <chr>     <chr>        <int> <dbl> <dbl> <dbl>    <dbl>
1 tz        tamsulosin   14184 0.708 0.638 0.787 1.22e-10

cox_fit(outcomes, tz_tam, "tz/dz/az", "tamsulosin")

# A tibble: 1 x 7
  treatment reference  n_pairs    hr    lb    ub p.value
  <chr>     <chr>        <int> <dbl> <dbl> <dbl>   <dbl>
1 tz        tamsulosin     754 0.751 0.574 0.984  0.0377

cox_fit(outcomes, NULL, "tz/dz/az", "5ari")

# A tibble: 1 x 7
  treatment reference n_pairs    hr    lb    ub p.value
  <chr>     <chr>       <int> <dbl> <dbl> <dbl>   <dbl>
1 tz        5ari         4063 0.819 0.708 0.947 0.00699

cox_fit(outcomes, tz_5ari, "tz/dz/az", "5ari")

# A tibble: 1 x 7
  treatment reference n_pairs    hr    lb    ub p.value
  <chr>     <chr>       <int> <dbl> <dbl> <dbl>   <dbl>
1 tz        5ari          475 0.833 0.610  1.14   0.253

cox_fit(outcomes, NULL, "tamsulosin", "5ari")

# A tibble: 1 x 7
  treatment  reference n_pairs    hr    lb    ub p.value
  <chr>      <chr>       <int> <dbl> <dbl> <dbl>   <dbl>
1 tamsulosin 5ari        13521  1.16  1.04  1.30 0.00775

cox_fit(outcomes, tam_5ari, "tamsulosin", "5ari")

# A tibble: 1 x 7
  treatment  reference n_pairs    hr    lb    ub p.value
  <chr>      <chr>       <int> <dbl> <dbl> <dbl>   <dbl>
1 tamsulosin 5ari          612  1.18 0.911  1.54   0.207

Time Varying

tz_tam_paired <- inner_join(outcomes, 
           tz_tam %>%
             mutate(enrolid = as.numeric(enrolid)), 
           by = "enrolid")
tz_5ari_paired <- inner_join(outcomes, 
           tz_5ari %>%
             mutate(enrolid = as.numeric(enrolid)), 
           by = "enrolid")
tam_5ari_paired <- inner_join(outcomes, 
           tam_5ari %>%
             mutate(enrolid = as.numeric(enrolid)), 
           by = "enrolid")

bind_rows(
  survfit(Surv(survival_time, event) ~ treatment, 
          data = outcomes %>% 
            filter(drug %in% c("tz/dz/az", "tamsulosin")) %>%
            mutate(treatment = as.numeric(drug == "tz/dz/az"))) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tz/Dz/Az", "Tamsulosin"),
      comparison = "Tz/Dz/Az vs Tamsulosin"
    ),
  survfit(Surv(survival_time, event) ~ treatment, 
          data = outcomes %>% 
            filter(drug %in% c("tz/dz/az", "5ari")) %>%
            mutate(treatment = as.numeric(drug == "tz/dz/az"))) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tz/Dz/Az", "5ARI"),
      comparison = "Tz/Dz/Az vs 5ARI"
    ),
  survfit(Surv(survival_time, event) ~ treatment, 
          data = outcomes %>% 
            filter(drug %in% c("5ari", "tamsulosin")) %>%
            mutate(treatment = as.numeric(drug == "tamsulosin"))) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tamsulosin", "5ARI"),
      comparison = "Tamsulosin vs 5ARI"
    )
  ) %>%
  mutate(comparison = forcats::fct_relevel(comparison, "Tz/Dz/Az vs Tamsulosin", "Tz/Dz/Az vs 5ARI")) %>%
  filter(time <= (365 * 7.5)) %>%
  ggplot(aes(x = time / 365, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_step(aes(color = drug)) + 
  geom_ribbon(aes(fill = drug), alpha = 0.33) + 
  facet_grid(cols = vars(comparison)) + 
  theme_minimal() + 
  theme(legend.position = "none") +
  labs(x = "Survival Time (Years)", y = "Fraction Remaining Dementia Free")

ggsave("~/raw.svg", width = 16/2, height = 9/2)

bind_rows(
  survfit(Surv(survival_time, event) ~ treatment, 
          data = tz_tam_paired) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tz/Dz/Az", "Tamsulosin"),
      comparison = "Tz/Dz/Az vs Tamsulosin"
    ),
  survfit(Surv(survival_time, event) ~ treatment, 
          data = tz_5ari_paired) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tz/Dz/Az", "5ARI"),
      comparison = "Tz/Dz/Az vs 5ARI"
    ),
  survfit(Surv(survival_time, event) ~ treatment, 
          data = tam_5ari_paired) %>%
    broom::tidy() %>%
    mutate(
      drug = ifelse(strata == "treatment=1", "Tamsulosin", "5ARI"),
      comparison = "Tamsulosin vs 5ARI"
    )
  ) %>%
  mutate(comparison = forcats::fct_relevel(comparison, "Tz/Dz/Az vs Tamsulosin", "Tz/Dz/Az vs 5ARI")) %>%
  filter(time <= (365 * 7.5)) %>%
  ggplot(aes(x = time / 365, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_step(aes(color = drug)) + 
  geom_ribbon(aes(fill = drug), alpha = 0.33) + 
  facet_grid(cols = vars(comparison)) + 
  theme_minimal() + 
  theme(legend.position = "none") +
  labs(x = "Survival Time (Years)", y = "Fraction Remaining Dementia Free")

ggsave("~/matched.svg", width = 16/2, height = 9/2)

survSplit(Surv(survival_time, event) ~ ., 
          data = tz_tam_paired,
          cut = c(1.5, 3) * 365,
          episode = "tgroup") %>%
  coxph(Surv(tstart, survival_time, event) ~ treatment:strata(tgroup),
        data = ., cluster = pair_id, robust = TRUE) %>%
  gtsummary::tbl_regression(exponentiate = TRUE)

Characteristic	HR¹	95% CI¹	p-value
treatment * strata(tgroup)
treatment * tgroup=1	0.54	0.36, 0.81	0.003
treatment * tgroup=2	1.16	0.68, 1.96	0.6
treatment * tgroup=3	0.91	0.54, 1.52	0.7
¹ HR = Hazard Ratio, CI = Confidence Interval

survSplit(Surv(survival_time, event) ~ ., 
          data = tz_5ari_paired,
          cut = c(1.5, 3) * 365,
          episode = "tgroup") %>%
  coxph(Surv(tstart, survival_time, event) ~ treatment:strata(tgroup),
        data = ., cluster = pair_id, robust = TRUE) %>%
  gtsummary::tbl_regression(exponentiate = TRUE)

Characteristic	HR¹	95% CI¹	p-value
treatment * strata(tgroup)
treatment * tgroup=1	0.67	0.43, 1.03	0.067
treatment * tgroup=2	1.48	0.78, 2.79	0.2
treatment * tgroup=3	0.73	0.38, 1.41	0.3
¹ HR = Hazard Ratio, CI = Confidence Interval

survSplit(Surv(survival_time, event) ~ ., 
          data = tam_5ari_paired,
          cut = c(1.5, 3) * 365,
          episode = "tgroup") %>%
  coxph(Surv(tstart, survival_time, event) ~ treatment:strata(tgroup),
        data = ., cluster = pair_id, robust = TRUE) %>%
  gtsummary::tbl_regression(exponentiate = TRUE)

Characteristic	HR¹	95% CI¹	p-value
treatment * strata(tgroup)
treatment * tgroup=1	0.99	0.68, 1.43	>0.9
treatment * tgroup=2	1.75	1.03, 2.97	0.037
treatment * tgroup=3	1.18	0.70, 2.00	0.5
¹ HR = Hazard Ratio, CI = Confidence Interval