--- 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") ```