chapter_two_scripts.R

library(tidyverse)
library(lubridate)
library(survival)
library(survminer)

#############################
### Making A Life Table In R!

# Step 1: Get Data
tmm <-read.csv(
  "https://raw.githubusercontent.com/ribsy/mdata/main/tmm_start_small.csv", stringsAsFactors = FALSE)

# Step 2: Tidy Data 
tmm <- tmm %>% mutate(status = ifelse(last.seen > "2020-12-31", 0, 1),
                      time.diff = ymd(last.seen) - ymd(first.seen) + 1)  %>%
  # filter(first.seen <= ymd("2020-10-01") &
  #       group == 1 & idTeam %in% 1 ) %>%
  select(time.diff, status)


# Step 3: Execute Survival Functions
t_surv <- survfit(data = tmm,
                  Surv(time.diff, status) ~ 1)

# Step 4a: Make Life Table
ltable <-
  tibble(day = t_surv$time) %>% 
  mutate(time_int     = str_c("[", day, ", ", day + 1, ")"), 
         n_opened     = t_surv$n.risk, 
         n_closed     = t_surv$n.event) %>% 
  mutate(n_censored   = n_opened - n_closed - lead(n_opened, default = 0),
         hazard_fun   = n_closed / n_opened,
         survivor_fun   = t_surv$surv,
         cum_hazard   = t_surv$cumhaz)

# Step 4b: Update Life Table - Add First Row
ltable <- bind_rows(tibble(day   = 0, 
                           time_int     = "[0, 1)", 
                           n_opened     = t_surv$n.risk[1], 
                           n_closed     = 0, 
                           n_censored   = 0, 
                           hazard_fun   = 0, 
                           survivor_fun   = 1,
                           cum_hazard   = 0), ltable)


##########################
### The Survival Functions

# The Survival Function Code Snippets 
t_surv <- survfit(data = tmm, Surv(time.diff, status) ~ 1)
t_surv

Surv(tmm$time.diff, tmm$status)[65:80]

Surv(tmm$time.diff, tmm$status)[(nrow(tmm)-15):nrow(tmm)]

############################
### Basic Life Table Metrics
max(ltable$day[(1 - ltable$survivor_fun) <= .5])

### Example with one Team and one Group
# Another downloadable data set with vuln severity
tmm <-read.csv("https://raw.githubusercontent.com/ribsy/mdata/main/tmm_start.csv")

# Step 2: Tidy Data 
tmm <- tmm %>% 
  
  # Add two new columns
  mutate(status = ifelse(last.seen > "2020-12-31", 0, 1),
         time.diff = ymd(last.seen) - ymd(first.seen) + 1)  %>%
  
  # Subset data – look for crit and high only
  filter(first.seen <= ymd("2020-10-01") & 
           group == 1 & 
           idTeam == 1 ) %>%
  
  # Remove these two columns
  select(time.diff,status,severity)

t_surv <- survfit(data = tmm,
                  Surv(time.diff, status) ~ severity)

ggsurvplot(t_surv,
           conf.int = TRUE,
           risk.table.col = "strata", # Change risk table color by groups
           ggtheme = theme_light(), # Change ggplot2 theme
           palette = c("coral3", "darkorange","deepskyblue4"),
           risk.table = "abs_pct",
           risk.table.y.text.col = T,
           fun = "event")