pcicabg-death.Rmd

---
title: "Death after PCI or CABG"
author: "Joonatan Palmu"
date: "`r format(Sys.time(), '%d.%m.%Y')`"
output: pdf_document
---
        
```{r setup, include=FALSE}
knitr::opts_chunk$set(include = TRUE, xecho = TRUE, message = FALSE, results='asis', cache=FALSE, warning=FALSE)
knitr::opts_chunk$set(cache.path = 'cache/', output.dir="cache/",
                      file.path = 'cache/', fig.path = 'cache/')
```

# Libraries

<details>
  <summary>Open/Close</summary>

```{r libraries, cache = FALSE}
library("haven")
library("magrittr")
library("dplyr")
library("survival")
library("knitr")
library("tibble")
library("tidyr")
library("broom")
library("paf")
library("pander")
library("survminer")
library("AF")
library("tableone")
library("purrr")
library("readr")
```

</details>

# Session info
 
 <details>
  <summary>Open/Close</summary>
 
```{r Session info, echo = FALSE}
pander(sessionInfo(), compact = TRUE)
```

</details>

# Helper functions

<details>
  <summary>Open/Close</summary>


```{r helper functions}
myread <- function(file, variables) {
  df <- read_sas(file)
  names(df) <- names(df) %>% toupper
  if (missing(variables))
      df
  else
      df %>% dplyr::select(variables)
}

kable <- function(data) {
  knitr::kable(data, booktabs = TRUE, digits = 2)
}

yeardiff <- function(start, end, days = 365.2425) {
    as.numeric(difftime(as.Date(strptime(start, format = "%d%m%Y")),
                        as.Date(strptime(end, format = "%d%m%Y")),
                        unit="days"))/days
}

mycharacteristics <- function(df, strata) {
    ret <- CreateTableOne(data = df,
                   strata = strata,
                   vars = c("CABG", "BL_AGE", "SUKUP", "BMI", "SYS",
                            "HDL", "LDL", "Q57", "SMOKER", "DIABETES",
                            "PREVAL_PAD", "PREVAL_COPD",
                            "PREVAL_HFAIL_STRICT"),
                   factorVars = c("SUKUP", "Q57", "SMOKER", "DIABETES",
                                  "CABG", "PREVAL_PAD",
                                  "PREVAL_COPD", "PREVAL_HFAIL_STRICT")
                   )
    capture.output(x <- print(ret, catDigits = 1, contDigits = 1))
    return(x)
}

pub.p <- function(p, Nfdr = FALSE) {
    p <- as.numeric(p)
    if (Nfdr) p <- p.adjust(p, method="BH", n = Nfdr)
    ifelse(p < 0.01, ifelse(p<0.001, "<0.001", sprintf("%.3f", p)), sprintf("%.3f", p))
}

c2l <- function(...) {
    l <- as.list(c(...))
    names(l) <- c(...)
    l
}
```

</details>

## Data import

<details>
  <summary>Open/Close</summary>

Reading follow up data and picking variables. If DEATH is reported
multiple times, take the one with smalles AGEDIFF.

```{r load followup data}
df.followup <- read_tsv("dset/followup/Release_FR_ALL_FU17Default_predefined_endpoints_2019-03-26_v1_1.txt",
                        col_types = cols(.default = col_double(), UID = col_character())) %>%
    filter(grepl("FR02_|FR07_|FR12_", UID)) %>%
    select("UID", "DEATH", "DEATH_AGE", "HFAIL_STRICT",
           "HFAIL_STRICT_AGE", "MI","MI_AGE", "STR", "STR_AGE", "PAD",
           "PAD_AGE", "COPD", "COPD_AGE")
```

Reading fr02, fr07, and fr12. Diabetes needs to be defined different ways between exams.

```{r fr data}
variables <- c("FR_ID", "VUOSI", "SUKUP", "BMI", "Q69", "HDL",
               "LDL", "SYS1", "SYS2", "SYS3", "DIAS1", "DIAS2",
               "DIAS3", "Q36", "Q57", "FR02_25", "FR02_25V",
               "FR02_24", "FR02_24V", "TULOAIKA", "Q67", "UID", "SYNTAIKA", "TUTKPVM")

fr02 <- myread('dset/FR02.sas7bdat', c(variables, "Q34")) %>%
    dplyr::mutate(DIABETES = as.factor(ifelse(Q34 == 3, 1, 0)),
                  COHORT = as.factor(1)) %>%
    dplyr::select(-Q34)

fr07 <- myread('dset/FR07.sas7bdat', c(variables, "FR07_38")) %>%
    dplyr::mutate(COHORT = as.factor(2),
                  DIABETES = as.factor(case_when(FR07_38 == 6 ~ 0,
                                                 FR07_38 == 5 ~ 1,
                                                 FR07_38 == 4 ~ 1,
                                                 FR07_38 == 3 ~ 1,
                                                 FR07_38 == 2 ~ 0,
                                                 FR07_38 == 1 ~ 0)),
                  UID = gsub(".$", "", UID)) %>%
    select(-FR07_38)

fr12 <- myread('dset/FR12.sas7bdat',
               c(variables, "FR12_20A", "FR12_20B", "FR12_20C", "FR12_20D",
                 "FR12_20E", "FR12_20F")) %>%
    dplyr::mutate(COHORT = as.factor(3),
                  DIABETES = as.factor(case_when(FR12_20F == 1 ~ 0,
                                       FR12_20E == 1 ~ 1,
                                       FR12_20D == 1 ~ 1,
                                       FR12_20C == 1 ~ 1,
                                       FR12_20B == 1 ~ 0,
                                       FR12_20A == 1 ~ 0))) %>%
    select(-FR12_20A, -FR12_20B, -FR12_20C, -FR12_20D, -FR12_20E, -FR12_20F)
```

Combining samples and for participants taking part in multiple studies, only the last one is included.

```{r merge fr}
dup.fr02 <- pull(fr02, FR_ID) %in% pull(fr07, FR_ID) |
    pull(fr02, FR_ID) %in% pull(fr12, FR_ID) 
dup.fr07 <- pull(fr07, FR_ID) %in% pull(fr12, FR_ID) 
df.unique <- rbind(fr12, fr07[!dup.fr07,], fr02[!dup.fr02,])
```

Merging sample data to followup data.

```{r merge without genetic}
df.merged <- right_join(df.unique,
                   df.followup,
                   by = "UID")
```

Taking the subset of participants that have PCI or CABG in history. Defining variables.

```{r mutations}
df.raw <- df.merged %>%
    dplyr::filter(FR02_25 == 2 | FR02_24 == 2) %>%
    dplyr::mutate(VUOSI = as.factor(VUOSI),
                  SYS = (SYS1 + SYS2 + SYS3)/3,
                  DIAS = (DIAS1 + DIAS2 + DIAS3)/3,
                  SMOKER = as.factor(case_when(Q67 == 1 ~ 0,
                                     Q69 == 1 ~ 1,
                                     Q69 == 2 ~ 1,
                                     Q69 == 3 ~ 0)),
                  INTERVENTION_TIME = min(FR02_24V, FR02_25V, na.rm = TRUE),
                  DEATH = as.integer(DEATH),
                  DIABETES = as.factor(DIABETES),
                  Q57 = as.factor(case_when(Q57 == 1 ~ 1,
                                     Q57 == 2 ~ 2,
                                     Q57 == 3 ~ 2,
                                     Q57 == 4 ~ 2)),
                  SUKUP = as.factor(SUKUP),
                  TULOAIKA = ifelse(TULOAIKA == "", NA, TULOAIKA),
                  COHORT = factor(COHORT, levels = c(1, 2, 3)),
                  CABG = as.factor(case_when(FR02_24 == 2 ~ 1,
                                             FR02_25 == 2 ~ 0)),
                  BL_AGE = yeardiff(TUTKPVM, SYNTAIKA),
                  DEATH_AGEDIFF = DEATH_AGE - BL_AGE,
                  PREVAL_PAD = as.factor(ifelse(PAD_AGE< BL_AGE & PAD == 1, 1 , 0)),
                  PREVAL_COPD = as.factor(ifelse(COPD_AGE < BL_AGE & COPD == 1, 1 , 0)),
                  PREVAL_HFAIL_STRICT = as.factor(ifelse(HFAIL_STRICT_AGE< BL_AGE &
                                                         HFAIL_STRICT == 1, 1 , 0)),
                  cohort = substr(UID, 1, 4)) %>% 
    select(-SYS1, -SYS2, -SYS3, -DIAS1, -DIAS2, -DIAS3, -Q36, -Q67, -Q69,
           -FR02_25V, -FR02_24V, -FR02_25, -FR02_24, -TULOAIKA)
```

Table of sample size and missing values

```{r number of nas, echo = FALSE}
df.raw %>%
    select_if(is.numeric) %>%
    gather(key, value) %>%
    group_by(key) %>%
    summarise_all(funs(N = n(),
                       Mean = mean(as.numeric(value), na.rm = TRUE),
                       SD = sd(as.numeric(value), na.rm = TRUE),
                       NAs = sum(is.na(value))))  %>%
    kable()
```

Excluding all participants with NAs.

```{r trimming}
df <- df.raw %>%
    na.omit()
```

Table of sample size after omiting NAs

```{r number of nas after na.omit, echo = FALSE}
df %>%
    select_if(is.numeric) %>%
    gather(key, value) %>%
    group_by(key) %>%
    summarise_all(funs(N = n(),
                       Mean = mean(as.numeric(value), na.rm = TRUE),
                       SD = sd(as.numeric(value), na.rm = TRUE),
                       NAs = sum(is.na(value))))  %>%
    kable()
```

</details>

# Results

## Characteristics

```{r characteristics, message = FALSE, echo = FALSE}
cbind(mycharacteristics(df),
      mycharacteristics(df, strata = "COHORT")) %>% kable
```


## Followup time

Total mean follow-up time is `r round(mean(df$DEATH_AGEDIFF), 2)` and median
`r round(median(df$DEATH_AGEDIFF), 2)`.

```{r number of events, echo = FALSE}
df %>%
    group_by(COHORT) %>%
    summarize(n = n(),
              dth = sum(DEATH == 1),
              alive = sum(DEATH == 0),
              percent = dth/alive,
              meantime = mean(DEATH_AGEDIFF)) %>%
    kable
```

## Death

Cox model for Death. Values reported as hazard ratios (exp(estimate)).

```{r cox death}
res.cox <- coxph(Surv(DEATH_AGEDIFF, DEATH) ~ 
                     BL_AGE + SUKUP + BMI + SYS + HDL + LDL + Q57 +
                     SMOKER + DIABETES + PREVAL_PAD + PREVAL_COPD +
                     PREVAL_HFAIL_STRICT + strata(VUOSI),
                 ties = "breslow",
                 data = df)

res.cox %>%
    broom::tidy(conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE) %>%
    arrange(p.value) %>%
    mutate(p.value = pub.p(p.value),
           betaci = sprintf("%.2f (%.2f−%.2f)", estimate, conf.low, conf.high)) %>%
    mutate_if(is.numeric, round, 2) %>%
    kable
```


### Assumptions

<details>
  <summary>Open/Close</summary>

Proportional hazard assumption

```{r cox assumption}
cox.zph(res.cox) %>%
    extract2("table") %>%
    kable
```

Influential observations

```{r influential observations}
ggcoxdiagnostics(res.cox,
                 type = "deviance",
                 linear.predictions = FALSE,
                 ggtheme = theme_bw())
```

Non linearity

```{r non-linearity}
survminer::ggcoxfunctional(formula = Surv(DEATH_AGEDIFF, DEATH) ~ 
                               BL_AGE + BMI + SYS, data = df)

```

</details>

### Testing effect of squared BMI

Testing if squared BMI is significant in previous model

```{r testing bmi2}
coxph(Surv(DEATH_AGEDIFF, DEATH) ~ 
                     BL_AGE + SUKUP + BMI + I(BMI**2) + SYS + HDL + LDL + Q57 +
                     SMOKER + DIABETES + PREVAL_PAD + PREVAL_COPD +
                     PREVAL_HFAIL_STRICT + strata(VUOSI),
                 ties = "breslow",
      data = df) %>%
    tidy %>%
    filter(grepl("BMI", term)) %>%
    kable
```

### Interactions with time

```{r continuous variable}
coxph(Surv(DEATH_AGEDIFF, DEATH) ~ 
                        BL_AGE + SUKUP + BMI + SYS + HDL + LDL + Q57 +
                        DIABETES + PREVAL_PAD + PREVAL_COPD +
                        PREVAL_HFAIL_STRICT + SMOKER*as.numeric(VUOSI),
                    ties = "breslow",
      data = df) %>%
    tidy %>%
    filter(grepl("SMOKER|VUOSI", term)) %>%
    kable
```

Checking if associations agaist outcomes has changed over time



```{r interaction with time}
vars <- c("BL_AGE", "SUKUP", "BMI", "SYS", "HDL", "LDL", "Q57", "SMOKER",
          "DIABETES", "PREVAL_PAD", "PREVAL_COPD",
          "PREVAL_HFAIL_STRICT")

res.strata <- coxph(Surv(DEATH_AGEDIFF, DEATH) ~ 
                        BL_AGE + SUKUP + BMI + SYS + HDL + LDL + Q57 +
                        SMOKER + DIABETES + PREVAL_PAD + PREVAL_COPD +
                        PREVAL_HFAIL_STRICT + strata(VUOSI),
                    ties = "breslow",
                    data = df)

ret.int <- lapply(c2l(vars), function(var, vars, compare) {
    fo <- sprintf("Surv(DEATH_AGEDIFF, DEATH) ~  %s + strata(VUOSI) + %s:VUOSI",
                  paste0(vars, collapse = " + "),
                  var)
    res.interaction <- coxph(as.formula(fo), ties = "breslow", data = df)
    anova(res.interaction, compare) %>% tidy %>% mutate(fo = fo)
}, vars = vars, compare = res.strata)

ret.int %>%
    map_df(., ~as.data.frame(.x), .id = "model") %>%
    filter(!is.na(p.value)) %>%
    select(model, fo, p.value) %>%
    arrange(p.value) %>%
    kable
```

## Attributable fraction function based on a Cox Proportional Hazard

[Manual for AF package](https://cran.r-project.org/web/packages/AF/AF.pdf)

```{r paf model}
df.af <- as.data.frame(df) %>%
    mutate(SMOKER = ifelse(SMOKER == 1, 1, 0))

res.af <- coxph(Surv(DEATH_AGEDIFF, DEATH) ~ BL_AGE + SUKUP + BMI +
                    SYS + HDL + LDL + Q57 + SMOKER + DIABETES +
                    PREVAL_PAD + PREVAL_COPD + PREVAL_HFAIL_STRICT +
                    VUOSI,
                ties = "breslow",
                data = df.af)
```

PAF for smoking

```{r paf smoker, message = FALSE}
AFcoxph(res.af, df.af, exposure = "SMOKER", time = c(5)) %>% summary %>%
    kableone
```