---
title: 'op_infections_ccasanet'
output:
  html_document:
    df_print: paged
  html_notebook: default
  word_document: default
  pdf_document: default
params: 
  printcode: FALSE
  date: !r lubridate::today()
always_allow_html: yes
---


#library(rms);library(tidyverse);library(lubridate);library(survival);library(tinytex);library(markdown);
#library(knitr);library(timereg);library(ff);library(pryr);library(survminer);library(cmprsk);library(stringr);
#library(kableExtra);library(flextable);library(DiagrammeR);library(rticles);library(agricolae);library(patchwork)

```{r Packages and databases are loaded here, include = T, echo = F, message=F, warning=F}
#Packages we'll use
library(rms);library(tidyverse);library(lubridate);library(survival);library(survminer);library(agricolae);library(patchwork);library(cmprsk)


#Databases we used
art <- read_csv("art.csv")
basic <- read_csv("basic.csv")
ce <- read_csv("ce.csv")
ce_cancer <- read_csv("ce_cancer.csv")
follow <- read_csv("follow.csv")
cd4 <- read_csv("lab_cd4.csv")
rna <- read_csv("lab_rna.csv")
visit <- read_csv("visit.csv")

```


```{r Here we filtered basic by patient age, date of enrolment, and presence in art&visit dataframes, include = T, echo = F}
basic_1 <-  basic %>%
  mutate(age_years = as.numeric(round((enrol_d-birth_d)/365,0))) %>%
  filter(age_years >= 18 & (year(enrol_d) >= "2001" & year(enrol_d) <= "2015"))

basic_filtered <- basic_1 %>%
  filter(patient_id %in% art$patient_id) %>%
  select(patient_id, birth_d, enrol_d, site, age_years, male, center)


#Number of patients at baseline
n_baseline <- n_distinct(basic$patient_id)

#Number of patients enroled between 2001 and 2015
n_enrol_year <- n_distinct(basic_1$patient_id)


#Number of patients with info on both the art database and the visit database
n_info_art_visit <- n_distinct(basic_filtered$patient_id)
```


```{r Haiti is excluded here because of its limited reporting of clinical events, include = T, echo = F}
ce_filtered <- ce %>%
  filter(patient_id %in% basic_filtered$patient_id) 
  

countries <- as.character(levels(as.factor(ce_filtered$site)))  
list1 <- vector("list", length(countries))
names(list1) <- countries

for(i in seq_along(countries)){
ce_1 <- ce_filtered %>%
  filter(site == countries[i])
list1[[i]] <- as.character(levels(as.factor(ce_1$ce_id)))
}

basic_filtered_2 <- basic_filtered %>%
  filter(site != "haiti")


#Number of patients from countries other than Haiti
n_not_haiti <- n_distinct(basic_filtered_2$patient_id);n_not_haiti
```


```{r Removing non-naive patients and those who didn´t have HAART as their first ART, include = T, echo =F}
art_filtered <- art %>%
  filter(str_count(art$art_id, ",") > 1) %>% #Filtering out those patients that didn't have HAART
  mutate(has_nnrti = if_else(str_detect(art_id, "EFV")|
                               str_detect(art_id, "ETR") |
                               str_detect(art_id, "NVP") |
                               str_detect(art_id, "RPV"), T, F),
         has_pi = if_else(str_detect(art_id, "LPV") | 
                          str_detect(art_id, "ATZ") |
                          str_detect(art_id, "RTV") |
                          str_detect(art_id, "IDV") |
                          str_detect(art_id, "SQV"), T, F),
         has_ii = if_else(str_detect(art_id, "RAL")|
                            str_detect(art_id, "EVG")|
                            str_detect(art_id, "DTG"), T, F)) %>% 
 filter(patient_id %in% basic_filtered_2$patient_id &
           art_class == "HAART") %>%
  group_by(patient_id) %>%
  arrange(art_sd, .by_group = T) %>%
  slice(1) %>%
  select(patient_id, art_id:ii2, art_class, art_sd, art_ed, has_nnrti, has_pi, has_ii)

basic_art <- basic_filtered_2 %>%
  left_join(art_filtered, by = "patient_id") %>%
  mutate(naive_enrol = if_else(art_sd >= enrol_d, T, F)) %>%
  filter(naive_enrol == T) #Joining dataframes and excluding non-naive patients

#Number of naive patients up to this point
n_only_naive <- n_distinct(basic_art$patient_id)
```


```{r Defining early ades and OIs DATE FILTER, include = T, echo = F}
ce_filtered2 <- ce_filtered %>%
  mutate(ade_mycosis = if_else(str_detect(ce_id, "candid")|
                                 str_detect(ce_id, "cryptoco")|
                                 str_detect(ce_id, "histo")|
                                 str_detect(ce_id, "pneumocy"), T, F),
         ade_crypto = if_else(str_detect(ce_id, "cryptoco"), T, F),
         is_ade = if_else(startsWith(ce_id, "ade"), T, F),
         is_oi = if_else(startsWith(ce_id, "ade") & !str_detect(ce_id, "cancer") &
                           !str_detect(ce_id, "wasting") & !str_detect(ce_id, "peds"), T, F)) %>%
  select(patient_id, ce_id, ce_d, ade_mycosis, ade_crypto, is_ade, is_oi) 


basic_art_ce <- basic_art %>%
  left_join(ce_filtered2, by = "patient_id") %>% 
  distinct()%>%
  group_by(patient_id) %>% 
  mutate(early_ade_there = ifelse(is_ade == T & (ce_d >= art_sd & ce_d <= art_sd + 180), T, F),
         early_ade = if_else(any(early_ade_there == T), T, F),
         which_early_ade = ce_id,
         early_ade_date = ce_d) #DF with info on clinical events and ART
```


```{r Including information on CD4 count at ART start DATE FILTER, include = F, echo = F, warning = F}
#The whole chunk must be run to correctly obtain the dataframe with CD4 at ART information
cd4_around_art <- basic_art_ce %>%
  left_join(select(cd4, cd4_d, cd4_v, patient_id, cd4_per), by = "patient_id") %>% 
  mutate(cd4_diff_absolute = if_else(cd4_d - art_sd >= 0 , cd4_d - art_sd , -(cd4_d - art_sd)),
         cd4_diff = art_sd - cd4_d) %>% 
  group_by(patient_id) %>% 
  arrange(cd4_diff_absolute, .by_group = T) %>% 
  mutate(low_cd4_prior = if_else(between(cd4_diff, -90, 90) & cd4_v < 200, T, F)) %>% 
  filter(between(cd4_diff, -90, 90)) %>% 
  slice(1)
         
cd4_last <- group_by(cd4, patient_id) %>% 
  arrange(desc(cd4_d), .by_group = T) %>% 
  mutate(last_date = cd4_d) %>% 
  slice(1) %>% 
  select(patient_id, last_date)

basic_all_cd4 <- basic_art_ce %>% 
  left_join(select(cd4_around_art, patient_id, cd4_d:low_cd4_prior), by = "patient_id") %>%
  left_join(cd4_last, by = "patient_id") %>% 
  slice(1)
```


```{r Adding survival time to oi or death or LTFU, include = T, echo = F, warning = F}

visit_follow_mod <- visit %>% 
              group_by(center, patient_id) %>% 
              arrange(desc(visit_d), .by_group = T) %>% 
              slice(1) %>% 
              ungroup() %>% 
  group_by(center) %>% 
  mutate(center_last_date = max(visit_d)) %>% 
  ungroup() %>% 
  mutate(date_visit = visit_d,
         time_to_last_visit = center_last_date - visit_d) %>% 
  group_by(patient_id) %>% 
  mutate(last_date = max(visit_d),
         last_visit = max(visit_d)) %>% 
  ungroup() %>% 
  select(patient_id, last_date, site, center, center_last_date, last_visit) %>% 
  full_join(follow, by = "patient_id") %>% 
  mutate(center= if_else(!is.na(center.x), center.x, center.y)) %>% 
  select(-center.x, -center.y, -site.y, -site.x) %>%
  group_by(center) %>% 
  mutate(center_last_date = if_else(center != "ar-ss" &!is.na(center_last_date), 
                                    center_last_date,
                                    max(l_alive_d))) %>% 
  ungroup()
  

rna_general <- group_by(rna, patient_id) %>% 
  arrange(desc(rna_d), .by_group=T) %>% 
  slice(1) %>% 
  mutate(last_date = rna_d) %>% 
  select(patient_id, last_date, center)

basic_preliminary_prelim <- basic_all_cd4 %>% 
  left_join(select(visit_follow_mod, patient_id, death_d, l_alive_d, drop_y, drop_d, last_date, last_visit, center_last_date),
            by = "patient_id") %>%
  left_join(ce %>% mutate(op_inf  = if_else(startsWith(ce_id, "ade") &
                                              !str_detect(ce_id, "cancer") &
                                              !str_detect(ce_id, "wasting") & !str_detect(ce_id, "peds"), T, F)) %>% 
              filter(op_inf==T), by = "patient_id") %>% #Aquí estoy haciendo que se queden puras oportunistas
  left_join(rna_general, by = "patient_id") %>% 
  group_by(patient_id) %>% 
  arrange(desc(last_date), .by_group = T) %>% 
  mutate(last_alive = if_else(l_alive_d >=last_date, l_alive_d, last_date),
         ltfu = if_else(center_last_date - last_alive >365, T, F),
         ltfu_date = as.Date(ifelse(ltfu == T, last_alive, NA), origin="1970-01-01"),#AGREGAR EL CONTRASTE DE LAS FECHAS DE LOS EVENTOS
         early_ade = if_else(is.na(ce_id.y), F, if_else(any(startsWith(ce_id.y, "ade") & 
                                                              between(difftime(ce_d.y, art_sd), 0, 180)), T, F)),
         last_alive_dif = as.numeric(last_alive - art_sd),
         op_inf_dif = ifelse(as.numeric(ce_d.y-art_sd)>=0, as.numeric(ce_d.y-art_sd), NA),
         op_inf = if_else(!is.na(op_inf_dif), T, F),
         ltfu_dif = as.numeric(ltfu_date-art_sd),
         death_dif = as.numeric(death_d-art_sd)) %>% 
  filter(ltfu_dif >= 180 | death_dif >=180 | op_inf_dif >= 180 | last_alive_dif >=180) %>% 
  mutate(event_oi = factor(if_else(op_inf == T & op_inf_dif >= 180, 
                           1,
                                   if_else(!is.na(death_dif) & (is.na(ltfu_dif) | death_dif <= ltfu_dif), 
                                           2, 
                                           if_else(!is.na(ltfu_dif), 
                                                   3, 0))), levels = c("0", "1", "2", "3")), 
         event_oi_surv = if_else(event_oi == 0, difftime(last_alive, art_sd, units = "days"),
                                 if_else(event_oi == 1, difftime(ce_d.y, art_sd, units = "days"),
                                         if_else(event_oi == 2, difftime(death_d, art_sd, units = "days"),
                                                 difftime(ltfu_date, art_sd, units = "days")))),
         event_date = if_else(event_oi == 1, ce_d.y, 
                              if_else(event_oi == 2, death_d,
                                      if_else(event_oi == 3, ltfu_date, last_alive)))) %>% 
  arrange(event_oi_surv, .by_group = T) %>% 
  slice(1) %>% 
  mutate(event_oi_surv = event_oi_surv-180)

#Here we filter out patients that died or were ltfu before completing six months of treatment (remembering that our time zero is actually 180 days)

#Number of patients in care 180 days after initiating ART
n_distinct(basic_preliminary_prelim$patient_id)
```


```{r Viral load information , include = F, echo = F, warnings = F, message=F}
#Here we add information on viral load both at ART start and at event. To correctly obtain the database 'basic_intermediate' (which contains this information) this whole chunk must be ran


rna_at_art <- basic_preliminary_prelim %>% 
  left_join(select(rna, rna_d, rna_v, patient_id), by = "patient_id") %>% 
  ungroup() %>% 
  mutate(rna_diff_art = art_sd-rna_d,
         rna_diff_art_absolute = if_else(art_sd - rna_d >= 0, art_sd - rna_d, -(art_sd - rna_d)),
         rna_at_art = if_else(rna_v < 0, 0, rna_v)) %>% 
  group_by(patient_id) %>% 
  arrange(rna_diff_art_absolute, .by_group = T) %>% 
  filter(rna_diff_art >= -7 & rna_diff_art <= 90) %>% 
  slice(1) %>% 
  select(patient_id, rna_diff_art:rna_at_art) %>% 
  mutate(rna_at_art_log_ten = log10(rna_at_art))


rna_at_6_months<- basic_preliminary_prelim %>% 
  left_join(select(rna, rna_d, rna_v, patient_id), by = "patient_id") %>% 
  ungroup() %>% 
  mutate(rna_diff_6m = (art_sd + 180) - rna_d,
         rna_diff_6m_absolute = if_else((art_sd + 180) - rna_d >= 0, (art_sd + 180) - rna_d, 
                                        -((art_sd+ 180) - rna_d)),
         rna_at_6m = rna_v) %>% 
  group_by(patient_id) %>% 
  arrange(rna_diff_6m_absolute, .by_group = T) %>% 
  filter(rna_diff_6m <= 7, rna_diff_6m >= -7) %>% 
  slice(1) %>% 
  mutate(rna_at_6_months_det = if_else(rna_at_6m > 200, T, F)) %>% 
  select(patient_id, rna_diff_6m:rna_at_6m, rna_at_6_months_det) %>% 
  mutate(rna_at_6_months_log_ten = log10(rna_at_6m))



#Here we obtain rna load near event for those patients who have at least a value (excluding those who have no viral load within the range, the range being at least 90 days posterior to the event or 180 days prior)
rna_at_event_range <- basic_preliminary_prelim %>% 
  left_join(select(rna, rna_d, rna_v, patient_id), by = "patient_id") %>% 
  ungroup() %>% 
  mutate(rna_diff_event = event_date-rna_d,
         rna_diff_event_absolute = if_else(event_date - rna_d >= 0, event_date - rna_d, -(event_date - rna_d)),
         rna_at_event = rna_v) %>% 
  group_by(patient_id) %>% 
  arrange(rna_diff_event_absolute, .by_group = T) %>% 
  filter(if_else(event_oi != 0, rna_diff_event >= -7 & rna_diff_event <= 7, !is.na(rna_diff_event))) %>% 
  slice(1) %>% 
  mutate(in_range_rna = T) %>% 
  select(patient_id, rna_diff_event:rna_at_event, in_range_rna) 
  

#Here we include rna load for those patients that have no rna load inside the specified range. The variable 'in_range' identifies such patients
rna_at_event_out_of_range <- basic_preliminary_prelim %>% 
  left_join(select(rna, rna_d, rna_v, patient_id), by = "patient_id") %>% 
  ungroup() %>% 
  filter(!(patient_id %in% rna_at_event_range$patient_id)) %>% 
  mutate(rna_diff_event = event_date-rna_d,
         rna_diff_event_absolute = if_else(event_date - rna_d >= 0, event_date - rna_d, -(event_date - rna_d)),
         rna_at_event = rna_v) %>% 
  group_by(patient_id) %>% 
  arrange(rna_diff_event_absolute, .by_group = T) %>% 
  slice(1) %>% 
  mutate(in_range_rna = F) %>% 
  select(patient_id, rna_diff_event:rna_at_event, in_range_rna)

rna_at_event_general <- rbind(rna_at_event_range, rna_at_event_out_of_range)


basic_intermediate <- basic_preliminary_prelim %>% 
  left_join(rna_at_art, by = "patient_id") %>% 
  left_join(rna_at_event_general, by = "patient_id") %>% 
  left_join(rna_at_6_months, by = "patient_id") 
```


```{r Adding mode of transmition:prior_aids:age:event_oi into final dataframe basic_cox, include = T, echo =F}
basic <- basic %>% 
  mutate(real_mode = if_else(str_detect(mode, "Homo") |
                             str_detect(mode, "Bisexual"), "Homosexual", 
                             if_else(str_detect(mode, "Hetero"), "hetero", 
                                             if_else(str_detect(mode, "Unknown"),
                                                     "unknown", "other"))))

#Adding info about transmition mode
basic_intermediate_1 <- basic_intermediate %>% 
  left_join(select(basic, patient_id, real_mode), by = "patient_id")

#Adding info about ADEs at or near ART sd
#Here Im only defining previous aids defining event by any clinical event date
ce_mod <- ce %>% 
  right_join(select(basic_intermediate, patient_id, art_sd), by = "patient_id") %>%
  group_by(patient_id) %>% 
  mutate(ade_prev_30 = if_else(any(art_sd-ce_d <= 30 & art_sd-ce_d >=0), T, F),
         ade_prev_any = if_else(any(art_sd-ce_d >= 0), T, F))


basic_almost <- basic_intermediate_1 %>% 
  left_join(select(ce_mod, patient_id, ade_prev_30, ade_prev_any), by = "patient_id")%>%
  distinct() %>% 
  left_join(select(basic, patient_id, aids_enrol_y), by = "patient_id") %>% 
  mutate(art_sd_year = year(art_sd),
         prior_aids = ifelse(cd4_v <200 | 
                               ade_prev_any == T | aids_enrol_y == 1, T, F),
         age_40 = if_else(age_years >= 40, T, F),
         event_oi_mod = if_else(event_oi == 0, "None", 
                                if_else(event_oi == 1, "LOI", 
                                        if_else(event_oi == 2, "Death", "LTFU"))),
         unknown_prior_aids = ifelse(is.na(prior_aids), T, F)) %>% 
  select(-site.y, -ce_id.x, -ce_d.x) %>% 
  rename(site = site.x, ce_id = ce_id.y, ce_d = ce_d.y)

```


```{r Adding information on CD4 count at event, include = T, echo = F}

cd4_at_event <- basic_almost %>% 
  rename(cd4_count_at_art = cd4_v,
         cd4_diff_at_art_absolute = cd4_diff_absolute,
         cd4_diff_at_art = cd4_diff,
         cd4_date_at_art = cd4_d) %>%
  left_join(select(cd4, cd4_d, cd4_v, patient_id), by = "patient_id") %>% 
  mutate(cd4_diff_absolute_event = if_else(event_date - cd4_d > 0 , event_date - cd4_d, -(event_date - cd4_d)),
         cd4_diff_at_event = event_date - cd4_d) %>% 
  rename(cd4_count_at_event = cd4_v,
         cd4_date_at_event = cd4_d) %>% 
  group_by(patient_id) %>% 
  filter(if_else(event_oi != 0, cd4_diff_at_event >= -30 | cd4_diff_at_event <= 30, cd4_diff_at_event >0)) %>% 
  arrange(cd4_diff_absolute_event, .by_group = T) %>% 
  slice(1) %>% 
  select(patient_id, cd4_date_at_event:cd4_diff_at_event)


basic_cox <- basic_almost %>% 
rename(cd4_count_at_art = cd4_v,
         cd4_diff_at_art_absolute = cd4_diff_absolute,
         cd4_diff_at_art = cd4_diff,
         cd4_date_at_art = cd4_d) %>%
  left_join(cd4_at_event, by = "patient_id") %>% #Here I'm adding info on CD4 count at event
  left_join(filter(art, patient_id %in% basic_almost$patient_id) %>% group_by(patient_id) %>% arrange(.by_group = T, art_sd) %>% 
              slice(1) %>% select(patient_id, art_rs), by = "patient_id") %>% 
  mutate(rna_at_event_det = if_else(rna_at_event>200, T, F),
         art_susp_ever = if_else(is.na(art_ed) | event_date < art_ed, F, T),#Here I'm adding information on treatment suspension
         art_susp_early = if_else(is.na(art_ed) | art_sd + 180 > art_ed, F, T),#OF NOTE, if someone has no end date for ART then I'm assuming
         art_susp_ever_f = if_else(art_susp_ever == T & any(str_detect(art_rs, c("failure" , "toxi" , "acidosis" ,
                                                                     "Side effects" , "Hypersens"))), T, F),#They're still on treatment
         art_susp_early_f = if_else(art_susp_early == T & any(str_detect(art_rs, c("failure" , "toxi" , "acidosis" ,
                                                                     "Side effects" , "Hypersens"))), T, F),
         event_bi = if_else(event_oi == 1, 1, 0),#This variable is for use in obtaining the p values on table_1
         death = if_else(!is.na(death_d), T, F),
         rna_at_art_in = if_else(rna_at_art < 200, T, F),
         real_mode_num = if_else(real_mode == "hetero", 0, 
                                 if_else(real_mode == "Homosexual", 1, 
                                         if_else(real_mode == "other", 2,3))),
         multiple_art = if_else((has_nnrti + has_pi + has_ii)>=2, T, F),
         real_mode = factor(real_mode, levels = c("Homosexual", "hetero", "other", "unknown"))) %>% 
  filter(rna_at_art_in == F|is.na(rna_at_art_in)) %>% 
  mutate(haart_type = if_else(has_nnrti == T, "NNRTI", 
                              if_else(has_pi == T, "PI", "II")),
         haart_type = factor(haart_type, levels = c("NNRTI","PI","II")),
         haart_type_numeric = if_else(haart_type == "NNRTI", 1, 
                                      if_else(haart_type == "PI",2, 3)))

```


```{r are the early ADE the same as the ones in the outcome, include = F}

basic_cox_only_loi <- filter(basic_cox, event_bi == 1) 
  

basic_art_ce_early_ade_not <- filter(basic_art_ce, early_ade ==T & patient_id %in% basic_cox_only_loi$patient_id) %>% 
  anti_join(select(basic_cox_only_loi, patient_id, ce_id), by = c( "patient_id", "ce_id"))


basic_art_ce_early_ade <- filter(basic_art_ce, early_ade ==T  
                                   & patient_id %in% basic_cox_only_loi$patient_id &
                                   !(patient_id %in% basic_art_ce_early_ade_not$patient_id)) %>% 
  left_join(basic_cox %>% rename(LOI_name = ce_id, LOI_date = ce_d) %>% 
              select(patient_id, LOI_date, LOI_name), by = "patient_id") %>% 
  select(patient_id, ce_d, ce_id, LOI_date, LOI_name, everything()) %>% 
  filter(LOI_date-ce_d <=365)
  
basic_intermediate <- filter(basic_cox, patient_id %in%basic_art_ce_early_ade$patient_id) %>% 
  mutate(event_bi = 0,
         event_oi = factor(if_else(!is.na(death_d),2,
                              if_else(center_last_date-l_alive_d >=365,3,0))),
         event_oi_mod = case_when(event_oi == 2~ "Death",
                                  event_oi==3~"LTFU",
                                  event_oi==0~"None"),
         event_oi_surv = l_alive_d-(art_sd+180)
  )
                              
                              
basic_cox <- rbind(filter(basic_cox, !(patient_id %in%basic_intermediate$patient_id)),basic_intermediate)
```


```{r Opportunistic infection count to include in results section, include = T, echo = F}
oi_count <- basic_cox %>% 
  ungroup() %>% 
  filter(startsWith(ce_id, "ade") & event_oi == 1) %>%
  count(ce_id) %>% 
  arrange(desc(n)) %>% 
  mutate(porcentaje = str_c(round((n/sum(n)*100),1), "%"),
         per = round((n/sum(n)*100),1))

medianas <- basic_cox %>% 
  ungroup() %>% 
  group_by(event_oi) %>% 
  summarise(mediana = median(event_oi_surv))
```


```{r Table 1,  warnings = F, include = F, echo = F}
#Vectors to include in the table
generalisimo <- function(x){
  c(n_distinct(x[,1]),  
    round(median(x$event_oi_surv[x$event_oi != 0])/365,1),
    NA, 
    sum(x$site == "argentina"), 
    sum(x$site == "brazil"), 
    sum(x$site == "chile"), 
    sum(x$site == "honduras"),
    sum(x$site == "mexico"),
    sum(x$site == "peru"),
    median(x$age_years),
    NA, 
    sum(x$male != 1), 
    sum(x$male), 
    NA, 
    sum(x$real_mode == "Homosexual"), 
    sum(x$real_mode == "hetero"),
    sum(x$real_mode == "other"),
    sum(x$real_mode == "unknown"),
    NA, 
    sum(x$has_nnrti), 
    sum(x$has_ii),
    sum(x$has_pi),
    median(x$art_sd_year),
    sum(x$prior_aids),
    sum(x$early_ade),
    sum(x$event_oi_mod == "LTFU", na.rm = T),
    sum(x$event_oi_mod == "Death", na.rm = T),
    NA,
    sum(x$rna_at_art < 200 , na.rm = T),
    sum(x$rna_at_6_months_det == F , na.rm = T),
    sum(x$rna_at_event_det == F, na.rm = T), 
    median(x$cd4_count_at_art, na.rm = T),
    median(x$cd4_count_at_event, na.rm = T),
    sum(x$art_susp_ever_f, na.rm = T),
    sum(x$art_susp_early_f, na.rm = T))
}


variables <- c("Count", "Median time of follow up/to event (years)", "Site", "Argentina", "Brazil", "Chile", "Honduras", "Mexico", "Peru", "Age", "Gender", 
                                      "Female", "Male",
                                      "Mode of acquisition", "Homosexual", "Heterosexual", "Other", "Unknown",
                                      "ART includes:", "Non nucleoside retrotranscriptase inhibitor", "Integrase inhibitor",
                                      "Protease inhibitor", "Year ART start", "AIDS diagnosis prior to ART start",
                                      "Early ADE (< 6 months)", "LTFU", "Died at some point", "Viral load (percentage undetectable)", 
               str_c("At ART start ","(", "n = ", sum(!is.na(basic_cox$rna_at_art)),")"),
               str_c("At 6 months ","(", "n = ", sum(!is.na(basic_cox$rna_at_6_months_det)),")"),
               str_c("At event ","(", "n = ", sum(!is.na(basic_cox$rna_at_event_det)),")"),
                                       "CD4 at ART start", "CD4 at event", "Suspended/changed ART prior to event",
               "Suspended/changed ART during the first 6 months of treatment")


porcentajes <- function(x){
  if_else(variables == "Age", 
          str_c(" ","(", fivenum(x$age_years, na.rm = T)[2], "-", 
                fivenum(x$age_years, na.rm = T)[4], ")"),
          if_else(variables == "Year ART start",
                  str_c(" ","(", fivenum(x$art_sd_year, na.rm = T)[2], "-", 
                        fivenum(x$art_sd_year, na.rm = T)[4], ")"), 
                  if_else(variables == "Median time of follow up/to event (years)",
                          str_c(" ","(", as.numeric(round(fivenum(x$event_oi_surv[x$event_oi != 0], na.rm = T)
                                               /365,1)[2]), "-", 
                                as.numeric(round(fivenum(x$event_oi_surv[x$event_oi != 0], na.rm = T)/365,1)[4])
                                , ")"),
                          if_else(variables == "CD4 at ART start",
                                  str_c(" ","(", fivenum(x$cd4_count_at_art, na.rm = T)[2], "-", 
                                        fivenum(x$cd4_count_at_art, na.rm = T)[4], 
                                        ", n = ", sum(!is.na(x$cd4_count_at_art)),
                                        ")"),
                                  
                                  if_else(variables == "CD4 at event",
                                          str_c(" ","(", fivenum(x$cd4_count_at_event, na.rm = T)[2], "-", 
                                                fivenum(x$cd4_count_at_event, na.rm = T)[4], 
                                                ", n = ", sum(!is.na(x$cd4_count_at_event)),
                                                ")"),
                                          str_c(" ","(",round(generalisimo(x)/n_distinct(x[,1])*100,1),"%",")"))))))
}



variables_p <-c(NA, "event_oi_surv", "site", NA, NA, NA, NA, NA , NA, "age_years", "male", NA, NA,
                                      "real_mode", NA, NA, NA, NA, NA, "has_nnrti", "has_ii","has_pi", "art_sd_year", "prior_aids",
                                      "early_ade", "drop_y", "death", NA, "rna_at_art", "rna_at_6_months_det", "rna_at_event_det",
                                       "cd4_count_at_art", "cd4_count_at_event", "art_susp_ever_f", "art_susp_early_f")


variables_p_2 <- variables_p[!is.na(variables_p)]
p_values_only <- variables_p
p_values_list <- variables_p
median_tabla_1 <- 1

#In this for loop we obtain results of statistical significance tests for table_1
for(i in 1:length(variables_p)){
  if(!is.na(variables_p[i]) & variables_p[i] == "event_oi_surv"){
    median_tabla_1 <- Median.test(y = basic_cox$event_oi_surv[basic_cox$event_oi != 0],
                        trt = basic_cox$event_oi[basic_cox$event_oi != 0])$statistics$p.chisq[[1]]
    p_values_list[i] <- median_tabla_1
  }  else  {
    if(!is.na(variables_p[i]) & (variables_p[i] == "age_years" |
       variables_p[i] == "art_sd_year" |
       variables_p[i] == "cd4_count_at_art" |
       variables_p[i] =="cd4_count_at_event")){
      p_values_list[i] <- wilcox.test(basic_cox[[variables_p[i]]] ~ event_bi, data = basic_cox)$p.value
    } else {
      if(!is.na(variables_p[i])){
      p_values_list[i] <- chisq.test(basic_cox[[variables_p[i]]], basic_cox$event_bi)$p.value
      
      } else {
      p_values_list[i] <- NA  
      }
    }
  }
}

p_values_list <- as.numeric(p_values_list)

p_values_only <- if_else(!is.na(p_values_list) & p_values_list < 0.001, "< 0.001", 
                         as.character(round(p_values_list, 3)))

#Here Im including two columns, one for all the cohort except those that developed an event and another one for those that did develop an event. 

table_1_df <- data.frame(Variable = variables,
                         'Did not develop OI' = str_c(generalisimo(basic_cox %>% filter(event_oi != 1)), " ",
                         porcentajes(basic_cox %>% filter(event_oi != 1))),
                         'Developed OI' = str_c(generalisimo(basic_cox %>% filter(event_oi == 1))," ", 
                                                porcentajes(basic_cox %>% filter(event_oi == 1))),
                         p_values = p_values_only)

View(table_1_df)

#Patients with prior aids
table(basic_cox$prior_aids, basic_cox$event_bi)
7192/(7192+612)
810/(810+80)
chisq.test(basic_cox$prior_aids, basic_cox$event_bi)
```


```{r univariated cox models for continuous variables, include = T, echo = F}

#Function to obtain cox HRs directly 
cox_univariado <- function(y){
  summary(coxph(Surv(event_oi_surv, event_oi == "1") ~ y + strata(basic_cox$site), data = basic_cox))
}

#univariated cox w/splines for CD4 count
uni_cd4 <- coxph(Surv(event_oi_surv, event_oi == "1") ~ rcs(cd4_count_at_art,3) + strata(basic_cox$site), data = basic_cox)

coxcd4_uv <- termplot(uni_cd4, se=TRUE, plot = F)$cd4_count_at_art
centercd4_uv <- with(coxcd4_uv, y[x == 300])
cicd4_uv <- coxcd4_uv$y + outer(coxcd4_uv$se, c(0, -1.96, 1.96), '*')

cd4u_df <- data.frame(cd4_count = coxcd4_uv$x, 
                     y = coxcd4_uv$y,
                     se = coxcd4_uv$se,
                     lower_ci = cicd4_uv[,2],
                     upper_ci = cicd4_uv[,3],
                     mean_ci = cicd4_uv[,1])

cd4_univar <- cd4u_df[with(cd4u_df, cd4_count == 100 | cd4_count == 200 | cd4_count == 300 | 
              cd4_count == 400 |cd4_count== 500),] %>% 
  mutate(hr = round(exp(y-centercd4_uv),2),
         upper_hr = round(exp(upper_ci-centercd4_uv),2),
         lower_hr = round(exp(lower_ci-centercd4_uv),2),
         full_ci = str_c(lower_hr, "-",upper_hr))

#univariated cox w/splines for age at art start
uni_age <- coxph(Surv(event_oi_surv, event_oi == "1") ~ rcs(age_years,3) + strata(basic_cox$site), data = basic_cox)

coxage_uv <- termplot(uni_age, se=TRUE, plot = F)$age_years
centerage_uv <- with(coxage_uv, y[x == 40])
ciage_uv <- coxage_uv$y + outer(coxage_uv$se, c(0, -1.96, 1.96), '*')

age_df <- data.frame(age = coxage_uv$x, 
                     y = coxage_uv$y,
                     se = coxage_uv$se,
                     lower_ci = ciage_uv[,2],
                     upper_ci = ciage_uv[,3],
                     mean_ci = ciage_uv[,1])

age_univar <- age_df[with(age_df, age == 20 | age == 30 | age == 40 | 
              age == 50 | age == 60),] %>% 
  mutate(hr = round(exp(y-centerage_uv),2),
         upper_hr = round(exp(upper_ci-centerage_uv),2),
         lower_hr = round(exp(lower_ci-centerage_uv),2),
         full_ci = str_c(lower_hr, "-",upper_hr))



##univariated cox w/splines for year of art start
uni_art_sd <- coxph(Surv(event_oi_surv, event_oi == "1") ~ rcs(art_sd_year,3) + strata(basic_cox$site), data = basic_cox)

coxart_sd_uv <- termplot(uni_art_sd, se=TRUE, plot = F)$art_sd_year
center_art_sd_uv <- with(coxart_sd_uv, y[x == 2010])
ci_art_sd_uv <- coxart_sd_uv$y + outer(coxart_sd_uv$se, c(0, -1.96, 1.96), '*')

art_sd_df <- data.frame(art_sd = coxart_sd_uv$x, 
                     y = coxart_sd_uv$y,
                     se = coxart_sd_uv$se,
                     lower_ci = ci_art_sd_uv[,2],
                     upper_ci = ci_art_sd_uv[,3],
                     mean_ci = ci_art_sd_uv[,1])

art_sd_univar <- art_sd_df[with(art_sd_df, art_sd == 2005 | art_sd == 2010 | art_sd == 2015),] %>% 
  mutate(hr = round(exp(y-center_art_sd_uv),2),
         upper_hr = round(exp(upper_ci-center_art_sd_uv),2),
         lower_hr = round(exp(lower_ci-center_art_sd_uv),2),
         full_ci = str_c(lower_hr, "-",upper_hr))


```


```{r cox multivariated models, include = T, echo = F}
cox_mv <- coxph(Surv(event_oi_surv, event_oi == 1)~  rcs(age_years, 3) + male + haart_type + early_ade + real_mode + art_sd_year  + art_susp_early_f + rcs(cd4_count_at_art, 3) + strata(site), data = basic_cox)

#Sex-> age at art start, cd4 count at art, early art switch
cox_mv_sex <- coxph(Surv(event_oi_surv, event_oi == 1)~  male + rcs(age_years, 3) + rcs(cd4_count_at_art, 3) + art_susp_early_f + strata(site), data = basic_cox) 

#Initial ART scheme-> age at art start, cd4 at art start, year of art start
cox_mv_art <- coxph(Surv(event_oi_surv, event_oi == 1)~  haart_type + rcs(age_years, 3) + art_sd_year  + rcs(cd4_count_at_art, 3) + strata(site), data = basic_cox) 

#Age of art start-> art start year, sex
cox_mv_age <- coxph(Surv(event_oi_surv, event_oi == 1)~  rcs(age_years, 3) + male + art_sd_year + strata(site), data = basic_cox) 

#ADE during the first 6 months of art start-> age of ART start, cd4 at art start, ART scheme
cox_mv_early_ade <- coxph(Surv(event_oi_surv, event_oi == 1)~ early_ade + rcs(age_years, 3) + haart_type + rcs(cd4_count_at_art, 3) + strata(site), data = basic_cox) 

#Mode of transmission-> no adjustment
cox_mv_mode <- coxph(Surv(event_oi_surv, event_oi == 1)~  real_mode + strata(site), data = basic_cox) 

#Year of art start_>No adjustment
cox_mv_year <- coxph(Surv(event_oi_surv, event_oi == 1)~  art_sd_year + strata(site), data = basic_cox) 

#CD4 count at art start-> age at ART start, year of ART start
cox_mv_cd4 <- coxph(Surv(event_oi_surv, event_oi == 1)~  rcs(cd4_count_at_art, 3) + rcs(age_years, 3) + art_sd_year + strata(site), data = basic_cox) 

#ART switch prior to event-> initial ART scheme, CD4 count at ART start
cox_mv_switch <- coxph(Surv(event_oi_surv, event_oi == 1)~  art_susp_early_f + haart_type+rcs(cd4_count_at_art, 3) + strata(site), data = basic_cox) 

```


```{r Graphs of continuous variables adjusted with splines in the multivariated model, echo = F}
#Cd4 count
coxcd4 <- termplot(cox_mv_cd4, se=TRUE, plot = F)$cd4_count_at_art
centercd4 <- with(coxcd4, y[x==301])
cicd4 <- coxcd4$y + outer(coxcd4$se, c(0, -1.96, 1.96), '*')

cd4_df <- data.frame(cd4_count = coxcd4$x, 
                     y = coxcd4$y,
                     se = coxcd4$se,
                     lower_ci = cicd4[,2],
                     upper_ci = cicd4[,3],
                     mean_ci = cicd4[,1])

cd4_spline_plot <- ggplot(cd4_df, aes(x = cd4_count, y = exp(mean_ci - centercd4))) +
  geom_hline(yintercept = c(.5,1,1.5,2,2.5), colour = "gray") +
  geom_smooth(se = F, colour = "red") +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(upper_ci - centercd4)), size = .1) +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(lower_ci - centercd4)), size = .1) +
  ylab("HR of LOI") +
  xlab("CD4 count at ART start (cells/mm3)")+
  theme_light()+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12)) +
  scale_y_continuous(breaks = (seq(0,3,.5))) +
  scale_x_continuous(breaks = (seq(0,1000,100)),
                     limits = c(0, 1000)) 
  


#Age at art start
cox_age <- termplot(cox_mv_age, se=TRUE, plot = F)$age_years
center_age <- with(cox_age, y[x==40])
ci_age <- cox_age$y + outer(cox_age$se, c(0, -1.96, 1.96), '*')

age_df <- data.frame(age = cox_age$x, 
                     y = cox_age$y,
                     se = cox_age$se,
                     lower_ci = ci_age[,2],
                     upper_ci = ci_age[,3],
                     mean_ci = ci_age[,1])

age_spline_plot <- ggplot(age_df, aes(x = age, y = exp(y - center_age))) +
  geom_hline(yintercept = seq(from = 0.5, to = 5, by = .5), colour = "gray") +
  geom_smooth(se = F, colour = "red") +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(upper_ci - center_age)), size = .1) +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(lower_ci - center_age)), size = .1) +
  ylab("HR of LOI") +
  xlab("Age at ART start (years)")+
  theme_light()+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12))+
  scale_x_continuous(breaks = seq(20,70,10),
                     limits = c(20, 70))+
  scale_y_continuous(limits = c(0, 2.5))


#Code for year of art start
cox_art <- termplot(cox_mv_year, se=TRUE, plot = F)$art_sd_year
center_art <- with(cox_art, y[x==2010])
ci_art <- cox_art$y + outer(cox_art$se, c(0, -1.96, 1.96), '*')

art_df <- data.frame(art_sd = cox_art$x, 
                     y = cox_art$y,
                     se = cox_art$se,
                     lower_ci = ci_art[,2],
                     upper_ci = ci_art[,3],
                     mean_ci = ci_art[,1])

art_spline_plot <- ggplot(art_df, aes(x = art_sd, y = exp(y - center_art))) +
  geom_hline(yintercept = seq(from = 0.5, to = 5, by = .5), colour = "gray") +
  geom_smooth(se = F, colour = "red") +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(upper_ci - center_art)), size = .1) +
  geom_smooth(se = F, colour = "red", linetype = "dashed", aes(y = exp(lower_ci - center_art)), size = .1) +
  ylab("HR of LOI") +
  xlab("Calendar year of ART start (years)")+
  theme_light()+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12))+
  scale_x_continuous(breaks = seq(2005,2015,5),
                     limits = c(2005, 2015))+
  scale_y_continuous(limits = c(0, 2.5))

#Joined graph of cd4 and age at art start
joined_plots <- cd4_spline_plot / 
  age_spline_plot /
  art_spline_plot +
  plot_annotation(tag_levels ="A") 
  
ggsave("joined_plots.jpeg", joined_plots, width = 100, height = 200, units = "mm")
```


```{r univariated cox models for categorical variables, include = F}
#This list does not include continuous variables adjusted with splines (which were defined in a previous chunk)
lista_cox_univ <- list(male = cox_univariado(y = basic_cox$male), 
                      early_ade = cox_univariado(basic_cox$early_ade), 
                      haart_type = cox_univariado(basic_cox$haart_type),
                      real_mode = cox_univariado(basic_cox$real_mode),
                      art_susp_ever_f = cox_univariado(basic_cox$art_susp_ever_f)
)


cox_u_df <- data.frame(variables = c("male", "has_nnrti", "has_pi", "has_ii", "age_years20",
                                   "age_years30", "age_years40",  "age_years50",
                                   "age_years60","early_ade", "real_modehomosexual", "real_modehetero", 
                                   "real_modeother", "real_modeunknown", "art_sd_year2005",
                                   "art_sd_year2010","art_sd_year2015",
                                   "cd4_100", "cd4_200", "cd4_300", "cd4_400", "cd4_500","art_susp_ever_f"),
                     hr = round(c(lista_cox_univ[["male"]]$conf.int[1],
                                  1,
                            lista_cox_univ[["haart_type"]]$conf.int["yPI", 1],
                            lista_cox_univ[["haart_type"]]$conf.int["yII", 2],
                            age_univar$hr,
                            lista_cox_univ[["early_ade"]]$conf.int[1],
                            1,
                            lista_cox_univ[["real_mode"]]$conf.int["yhetero",1],
                            lista_cox_univ[["real_mode"]]$conf.int["yother",1],
                            lista_cox_univ[["real_mode"]]$conf.int["yunknown",1],
                            art_sd_univar$hr,
                            cd4_univar$hr,
                            lista_cox_univ[["art_susp_ever_f"]]$conf.int[1]),2),
                     ci = c(str_c(round(lista_cox_univ[["male"]]$conf.int[3],2),"-",
                                        round(lista_cox_univ[["male"]]$conf.int[4],2)),
                            str_c("Reference"),
                            str_c(round(lista_cox_univ[["haart_type"]]$conf.int["yPI", 3],2),"-",
                                  round(lista_cox_univ[["haart_type"]]$conf.int["yPI", 4],2)),
                            str_c(round(lista_cox_univ[["haart_type"]]$conf.int["yII", 3],2),"-",
                                  round(lista_cox_univ[["haart_type"]]$conf.int["yII", 4],2)),
                            age_univar$full_ci,
                            str_c(round(lista_cox_univ[["early_ade"]]$conf.int[3],2),"-",
                                  round(lista_cox_univ[["early_ade"]]$conf.int[4],2)),
                            "Reference",
                            str_c(round(lista_cox_univ[["real_mode"]]$conf.int["yhetero",3],2),"-",
                                  round(lista_cox_univ[["real_mode"]]$conf.int["yhetero",4],2)),
                            str_c(round(lista_cox_univ[["real_mode"]]$conf.int["yother",3],2),"-",
                                  round(lista_cox_univ[["real_mode"]]$conf.int["yother",4],2)),
                            str_c(round(lista_cox_univ[["real_mode"]]$conf.int["yunknown",3],2),"-",
                                  round(lista_cox_univ[["real_mode"]]$conf.int["yunknown",4],2)),
                           art_sd_univar$full_ci,
                           cd4_univar$full_ci,
                            str_c(round(lista_cox_univ[["art_susp_ever_f"]]$conf.int[3],2),"-",
                                  round(lista_cox_univ[["art_susp_ever_f"]]$conf.int[4],2))))
           
cox_u_df[] <- lapply(cox_u_df, as.character)

```


```{r Multivariated cox model info for table 2, include = F}
#HRs for CD4s
cd4_hrs <- cd4_df[with(cd4_df, cd4_count == 100 | cd4_count == 200 | cd4_count == 301 | 
                         cd4_count == 400 |cd4_count== 499),] %>% 
  mutate(hr = exp(y-centercd4),
         upper_hr = exp(upper_ci-centercd4),
         lower_hr = exp(lower_ci-centercd4))

#HRs for age
age_hrs <- age_df[with(age_df, age == 20 | age == 30 | age == 40 | age == 50 | age == 60),] %>% 
  mutate(hr = exp(y-center_age),
                                                                                                                             upper_hr = exp(upper_ci-center_age),
                                                                                                                             lower_hr = exp(lower_ci-center_age))

#HRs for art start year
artsd_hrs <- art_df[with(art_df, art_sd == 2005 | art_sd == 2010 | art_sd == 2015),] %>% 
  mutate(hr = exp(y-center_art),
         upper_hr = exp(upper_ci-center_art),
         lower_hr = exp(lower_ci-center_art))

#Information for the data frame below will come from cox_mv for most of the variables. For cd4_count_at_art, age_years and art_sd_year it will be from the individual data frames I created which hold their respective HRs and confidence intervals.

cox_mv_df <- data.frame(variables = c("male", "has_nnrti", "has_pi", "has_ii", "age_years20",
                                   "age_years30", "age_years40",  "age_years50",
                                   "age_years60","early_ade", "real_modehomosexual","real_modehetero", 
                                   "real_modeother", "real_modeunknown", "art_sd_year2005",
                                   "art_sd_year2010","art_sd_year2015",
                                   "cd4_100", "cd4_200", "cd4_300", "cd4_400", "cd4_500","art_susp_early_f"),
                     hr = round(c(summary(cox_mv_sex)$conf.int["male",][1],
                            1,
                            summary(cox_mv_art)$conf.int["haart_typePI",][1],
                            summary(cox_mv_art)$conf.int["haart_typeII",][1],
                            age_hrs$hr,
                            summary(cox_mv_early_ade)$conf.int["early_adeTRUE",][1],
                            1,
                            summary(cox_mv_mode)$conf.int["real_modehetero",][1],
                            summary(cox_mv_mode)$conf.int["real_modeother",][1],
                            summary(cox_mv_mode)$conf.int["real_modeunknown",][1],
                            artsd_hrs$hr,
                            cd4_hrs$hr,
                            summary(cox_mv_switch)$conf.int["art_susp_early_fTRUE",][1]),2),
                     ci = c(str_c(round(summary(cox_mv_sex)$conf.int["male",][3],2),"-",
                                        round(summary(cox_mv_sex)$conf.int["male",][4],2)),
                            str_c("Reference"),
                            str_c(round(summary(cox_mv_art)$conf.int["haart_typePI",][3],2),"-",
                                  round(summary(cox_mv_art)$conf.int["haart_typePI",][4],2)),
                            str_c(round(summary(cox_mv_art)$conf.int["haart_typeII",][3],2),"-",
                                  round(summary(cox_mv_art)$conf.int["haart_typeII",][4],2)),
                            str_c(round(age_hrs$lower_hr,2), "-", round(age_hrs$upper_hr,2)),
                            str_c(round(summary(cox_mv_early_ade)$conf.int["early_adeTRUE",][3],2),"-",
                                  round(summary(cox_mv_early_ade)$conf.int["early_adeTRUE",][4],2)),
                            "Reference",
                            str_c(round(summary(cox_mv_mode)$conf.int["real_modehetero",][3],2),"-",
                                  round(summary(cox_mv_mode)$conf.int["real_modehetero",][4],2)),
                            str_c(round(summary(cox_mv_mode)$conf.int["real_modeother",][3],2),"-",
                                  round(summary(cox_mv_mode)$conf.int["real_modeother",][4],2)),
                            str_c(round(summary(cox_mv_mode)$conf.int["real_modeunknown",][3],2),"-",
                                  round(summary(cox_mv_mode)$conf.int["real_modeunknown",][4],2)),
                            str_c(round(artsd_hrs$lower_hr,2), "-", round(artsd_hrs$upper_hr, 2)),
                            str_c(round(cd4_hrs$lower_hr,2), "-", round(cd4_hrs$upper_hr,2)),
                            str_c(round(summary(cox_mv_switch)$conf.int["art_susp_early_fTRUE",][3],2),"-",
                                  round(summary(cox_mv_switch)$conf.int["art_susp_early_fTRUE",][4],2))))
           
cox_mv_df[]  <- lapply(cox_mv_df, as.character)
```


```{r Competing risk analysis for table 2, include = F}
#All variables coded numerically
cov1 <- model.matrix(~
                       male + early_ade + age_years + haart_type +
                      real_mode_num + art_sd_year +  art_susp_early_f+
                       cd4_count_at_art,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 


cov_ej <- model.matrix(~
                       male + early_ade + age_years + haart_type +
                      real_mode_num + art_sd_year +  art_susp_early_f+
                       rcs(cd4_count_at_art,3),
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1)


hr_sub_ej <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov_ej)

#Sex-> age at art start, cd4 count at art, early art switch
cov1_sex <- model.matrix(~
                       male + age_years + haart_type +art_susp_early_f+
                       cd4_count_at_art,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_sex <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_sex)

#Initial ART scheme-> age at art start, cd4 at art start, year of art start
cov1_art <- model.matrix(~
                       haart_type + age_years  +
                      art_sd_year +  
                       cd4_count_at_art,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_art <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_art)

#Age of art start-> art start year, sex
cov1_age <- model.matrix(~
                      age_years + male + art_sd_year,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_age <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_age)

#ADE during the first 6 months of art start-> age of ART start, cd4 at art start, ART scheme
cov1_ade <- model.matrix(~
                       early_ade + age_years + haart_type +
                       cd4_count_at_art,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_ade <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_ade)

#Mode of transmission-> no adjustment
cov1_mode <- model.matrix(~
                      real_mode,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_mode <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_mode)

#Year of art start->No adjustment
cov1_year <- model.matrix(~
                      art_sd_year,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_year <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_year)

#CD4 count at art start-> age at ART start, year of ART start
cov1_cd4 <- model.matrix(~
                       cd4_count_at_art +age_years + art_sd_year,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_cd4 <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_cd4)


#ART switch prior to event-> initial ART scheme, CD4 count at ART start
cov1_switch <- model.matrix(~
                       art_susp_early_f + haart_type+ cd4_count_at_art,
                     data = filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)))[, -1] 

hr_sub_oi_switch <- crr(ftime=with(filter(basic_cox, !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                            event_oi_surv), 
            fstatus=with(filter(basic_cox,  !is.na(cd4_count_at_art)
                                   & !is.na(art_susp_early_f)),
                         event_oi), 
            cov1=cov1_switch)



#Predictions for CD4
crf_df_cd4 <- data.frame(cd4_count_at_art = seq(100, 500, 100),
                         age_years = rep(40,5),
                         art_sd_year = rep(2010, 5))
 
crf_pred_cd4 <- predict(hr_sub_oi_cd4, cov1 = as.matrix(crf_df_cd4))

crf_pred_cd4_ej <- predict(hr_sub_ej, cov1 = as.matrix(crf_df_cd4))
 
#Each object here takes as reference a cd4 count of 300
 crf_hr_cd4_100 <-round(crf_pred_cd4[625,2]/crf_pred_cd4[625,4],2)
 crf_hr_cd4_200 <-round(crf_pred_cd4[625,3]/crf_pred_cd4[625,4],2)
 crf_hr_cd4_400 <-round(crf_pred_cd4[625,5]/crf_pred_cd4[625,4],2)
 crf_hr_cd4_500 <-round(crf_pred_cd4[625,6]/crf_pred_cd4[625,4],2)

crf_hr_cd4_100_ej <-round(crf_pred_cd4_ej[625,2]/crf_pred_cd4_ej[625,4],2)
 crf_hr_cd4_200_ej <-round(crf_pred_cd4_ej[625,3]/crf_pred_cd4_ej[625,4],2)
 crf_hr_cd4_400_ej <-round(crf_pred_cd4[625,5]/crf_pred_cd4[625,4],2)
 crf_hr_cd4_500_ej <-round(crf_pred_cd4[625,6]/crf_pred_cd4[625,4],2)
 
#HR and conf interval for cd4 count
ci_risk_100 <- str_c(round(crf_hr_cd4_100 - (1.96*.0003),2),"-",
       round(crf_hr_cd4_100 + (1.96*.0003),2))
 
ci_risk_200 <- str_c(round(crf_hr_cd4_200 - (1.96*.0003),2),"-",
       round(crf_hr_cd4_200 + (1.96*.0003),2))
 
ci_risk_400 <- str_c(round(crf_hr_cd4_400 - (1.96*.0003),2),"-",
       round(crf_hr_cd4_400 + (1.96*.0003),2))
 
ci_risk_500 <- str_c(round(crf_hr_cd4_500 - (1.96*.0003),2),"-",
       round(crf_hr_cd4_500 + (1.96*.0003),2))
#
 
 str_c(crf_hr_cd4_500 - (1.96*.0003),"-",
       crf_hr_cd4_500 + (1.96*.0003))

#Predictions for age at art start
 crf_df_age <- data.frame(age_years = c(20,30,40,50,60),
                          male =  rep(T, 5), 
                          art_sd_year = rep(2015, 5))
 
crf_pred_age <- predict(hr_sub_oi_age, cov1 = as.matrix(crf_df_age))
#Age 40 is the reference
 crf_hr_age_20 <-round(crf_pred_age[625,2]/crf_pred_age[625,4], 2)
 crf_hr_age_30 <-round(crf_pred_age[625,3]/crf_pred_age[625,4], 2)
 crf_hr_age_50 <-round(crf_pred_age[625,5]/crf_pred_age[625,4], 2)
 crf_hr_age_60 <-round(crf_pred_age[625,6]/crf_pred_age[625,4], 2)

 
#HR and conf int for age at art start
ci_risk_20 <- str_c(round(crf_hr_age_20 - (1.96*.00368), 2),"-",
       round(crf_hr_age_20 + (1.96*.00368),2))

ci_risk_30 <- str_c(round(crf_hr_age_30 - (1.96*.00368), 2),"-",
       round(crf_hr_age_30 + (1.96*.00368),2))

ci_risk_50 <- str_c(round(crf_hr_age_50 - (1.96*.00368), 2),"-",
       round(crf_hr_age_50 + (1.96*.00368),2))

ci_risk_60 <- str_c(round(crf_hr_age_60 - (1.96*.00386), 2),"-",
       round(crf_hr_age_60 + (1.96*.00386),2))
#

#Predictions for year of art start
 crf_df_year <- data.frame(art_sd_year = c(2005,2010,2015))
 
crf_pred_year <- predict(hr_sub_oi_year, cov1 = as.matrix(crf_df_year))
#Year 2010 is the reference
 crf_hr_year_2005 <- round(crf_pred_year[377,2]/crf_pred_year[377,3], 2)
 crf_hr_year_2015 <- round(crf_pred_year[377,4]/crf_pred_year[377,3], 2)
 
#HRs and conf int for year of art start
ci_risk_2005 <- str_c(round(crf_hr_year_2005 - (1.96*.01013), 2),"-",
       round(crf_hr_year_2005 + (1.96*.01013),2))

ci_risk_2015 <- str_c(round(crf_hr_year_2015 - (1.96*.01013), 2),"-",
       round(crf_hr_year_2015 + (1.96*.01013),2))
#


#Subdistribution hazard ratios data frame
hr_sub_df <- data.frame(variables = c("male", "has_nnrti", "has_pi", "has_ii", "age_years20",
                                   "age_years30", "age_years40",  "age_years50",
                                   "age_years60","early_ade", "real_mode_homosexual", "real_modehetero", 
                                   "real_modeother", "real_modeunknown", "art_sd_year2005",
                                   "art_sd_year2010","art_sd_year2015",
                                   "cd4_100", "cd4_200", "cd4_300", "cd4_400", "cd4_500","art_susp_early_f"),
                        hr = round(c(summary(hr_sub_oi_sex)$conf.int["male",1],
                            1,
                            summary(hr_sub_oi_art)$conf.int["haart_typePI",1],
                            summary(hr_sub_oi_art)$conf.int["haart_typeII",1],
                            crf_hr_age_20, crf_hr_age_30,1, crf_hr_age_50, crf_hr_age_60,
                            summary(hr_sub_oi_ade)$conf.int["early_adeTRUE",1],
                            1,
                            summary(hr_sub_oi_mode)$conf.int["real_modehetero",1],
                            summary(hr_sub_oi_mode)$conf.int["real_modeother",1],
                            summary(hr_sub_oi_mode)$conf.int["real_modeunknown",1],
                            crf_hr_year_2005, 1, crf_hr_year_2015,
                            crf_hr_cd4_100, crf_hr_cd4_200, 1, crf_hr_cd4_400, crf_hr_cd4_500,
                            summary(hr_sub_oi_switch)$conf.int["art_susp_early_fTRUE",1]),2),
                     ci = c(str_c(round(summary(hr_sub_oi_sex)$conf.int["male",3],2),"-",
                                        round(summary(hr_sub_oi_sex)$conf.int["male",4],2)),
                            str_c("Reference"),
                            str_c(round(summary(hr_sub_oi_art)$conf.int["haart_typePI",3],2),"-",
                                  round(summary(hr_sub_oi_art)$conf.int["haart_typePI",4],2)),
                            str_c(round(summary(hr_sub_oi_art)$conf.int["haart_typeII",3],2),"-",
                                  round(summary(hr_sub_oi_art)$conf.int["haart_typeII",4],2)),
                            ci_risk_20, ci_risk_30, "-", ci_risk_50, ci_risk_60,
                            str_c(round(summary(hr_sub_oi_ade)$conf.int["early_adeTRUE",3],2),"-",
                                  round(summary(hr_sub_oi_ade)$conf.int["early_adeTRUE",4],2)),
                            "Reference",
                            str_c(round(summary(hr_sub_oi_mode)$conf.int["real_modehetero",3],2),"-",
                                  round(summary(hr_sub_oi_mode)$conf.int["real_modehetero",4],2)),
                            str_c(round(summary(hr_sub_oi_mode)$conf.int["real_modeother",3],2),"-",
                                  round(summary(hr_sub_oi_mode)$conf.int["real_modeother",4],2)),
                            str_c(round(summary(hr_sub_oi_mode)$conf.int["real_modeunknown",3],2),"-",
                                  round(summary(hr_sub_oi_mode)$conf.int["real_modeunknown",4],2)),
                            ci_risk_2005, "-",ci_risk_2015,
                            ci_risk_100, ci_risk_200, "-", ci_risk_400, ci_risk_500,
                            str_c(round(summary(hr_sub_oi_switch)$conf.int["art_susp_early_fTRUE",3],2),"-",
                                  round(summary(hr_sub_oi_switch)$conf.int["art_susp_early_fTRUE",4],2))))


hr_sub_df[] <- lapply(hr_sub_df, as.character)
```


```{r Table 2, echo = F}
#Names of the variables to be included in table 2
variables_2 <- c(NA, "male", "has_nnrti", "has_pi", "has_ii", "age_years20",
                                   "age_years30", "age_years40",  "age_years50",
                                   "age_years60","early_ade", "real_modeHomosexual", "real_modehetero", 
                                   "real_modeother", "real_modeunknown", "art_sd_year2005",
                                   "art_sd_year2010","art_sd_year2015",
                                   "cd4_100", "cd4_200", "cd4_300", "cd4_400", "cd4_500","art_susp_early_f")


table_2_df <- data.frame(Variable = variables_2,
                         "Unadjusted HR" = c(NA, 
                                             str_c(round(as.numeric(cox_u_df$hr), 2), " (", as.character(cox_u_df$ci), ")")),
                         'Adjusted HR' = c(NA, 
                                        str_c(round(as.numeric(cox_mv_df$hr),2), " (", cox_mv_df$ci, ")")),
                         'Subdistribution HR' = c(NA, 
                                       str_c(hr_sub_df$hr, " (", as.character(hr_sub_df$ci), ")")))


```


```{r Estimates cumulative incidence graphs, echo = F}
#Estimated cumulative incidence at various time points with OI as the event of interest
comp_risk_oi_2 <- cuminc(ftime = basic_cox$event_oi_surv/365, fstatus = factor(basic_cox$event_oi_mod,
                                                                               levels = c( "Death", "LOI","LTFU", "None")),
                         cencode = "None", group = "")



#Cumulative incidence graphs
cuminc_graph_2 <- ggcompetingrisks(comp_risk_oi_2, palette = "lancet",
                 legend = "top",
                 ggtheme = theme_survminer(),
                 multiple_panels = T,
                 xlab = "Time (Years)",
                 ylab = "Estimated cumulative incidence",
                 legend.title = "Event:",
                 xticks.by = 1,
                 xlim = c(0.5,15),
                 ylim = c(0,1),
                 conf.int = T)+
  ggtitle(label = "")

ggsave("plot_figure_1.jpg", cuminc_graph_2, width = 200, height = 200, units = "mm")
```

