# Clearing space and loading libraries rm(list=ls()) options(width=180) #------------------------------# # Loading libraries # #------------------------------# library(Hmisc) library(rms) library(ggplot2) library(cmprsk) #------------------------------# # Loading functions # #------------------------------# source('pediatric_mortality_functions.R') art <- read.csv('Data/art_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) basic <- read.csv('Data/basic_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) follow <- read.csv('Data/follow_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) cd4 <- read.csv('Data/lab_cd4_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) rna <- read.csv('Data/lab_rna_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) visit <- read.csv('Data/visit_p.csv',header=TRUE,na.strings=c(NA,''),stringsAsFactors=FALSE) #-------------------------------------------# # Creating unique patient identifiers # #-------------------------------------------# art$patient <- tolower(art$patient) basic$patient <- tolower(basic$patient) follow$patient <- tolower(follow$patient) cd4$patient <- tolower(cd4$patient) rna$patient <- tolower(rna$patient) visit$patient <- tolower(visit$patient) art <- upData(art, rename=c('patient'='unique_id')) basic <- upData(basic, rename=c('patient'='unique_id', 'baseline_d'='enrol_d')) # Changed baseline_d name because of changes made in Mark's ADE at entry code for first_failure3.nw follow <- upData(follow, rename=c('patient'='unique_id')) cd4 <- upData(cd4, rename=c('patient'='unique_id')) rna <- upData(rna, rename=c('patient'='unique_id')) visit <- upData(visit, rename=c('patient'='unique_id')) # Making regimen classes more manageable art$art_class <- tolower(art$art_class) #-----------------------------------# # Formatting dates # #-----------------------------------# art$art_sd <- ifelse(art$art_sd == '1990-01-01',NA,art$art_sd) art$art_sd <- as.Date(art$art_sd, format='%Y-%m-%d') art$art_ed <- as.Date(art$art_ed, format='%Y-%m-%d') follow$death_d <- as.Date(follow$death_d, format='%Y-%m-%d') follow$l_alive_d <- as.Date(follow$l_alive_d, format='%Y-%m-%d') basic$enrol_d <- as.Date(basic$enrol_d, format='%Y-%m-%d') basic$hivdiagnosis_d <- as.Date(basic$hivdiagnosis_d, format='%Y-%m-%d') basic$birth_d <- as.Date(basic$birth_d, format='%Y-%m-%d') cd4$cd4_d <- as.Date(cd4$cd4_d, format='%Y-%m-%d') rna$rna_d <- as.Date(rna$rna_d, format='%Y-%m-%d') #-----------------------------------# # ART in relation to birth # #-----------------------------------# # Removing records from ART with missing start dates art <- subset(art, !is.na(art_sd)) art$birth_d <- basic[match(art$unique_id, basic$unique_id, nomatch=NA),'birth_d'] art$days_from_birth <- as.numeric(art$art_sd - art$birth_d) art$days_on_art <- as.numeric(art$art_ed - art$art_sd) art$within_16_days <- ifelse(art$days_from_birth <= 16,1,0) #-----------------------------------# # First HAART # #-----------------------------------# art <- art[order(art$unique_id, art$art_sd),] which.first.haart <- which(art$art_class == 'haart') fhaart <- art[which.first.haart,] fhaart <- subset(fhaart, !duplicated(unique_id)) fhaart$year_fhaart <- apply(subset(fhaart, select='art_sd'),1,FUN=function(y){ tmp1 <- strsplit(as.character(y['art_sd']),split='-')[[1]] tmp_yr <- tmp1[1] return(tmp_yr) }) fhaart$year_fhaart <- as.numeric(fhaart$year_fhaart) fhaart$fhaart_class <- ifelse(fhaart$nnrti >= 1 & fhaart$pi == 0,1, ifelse(fhaart$nnrti == 0 & fhaart$pi >= 1,2,3)) fhaart$fhaart_class.factor <- factor(fhaart$fhaart_class, levels=c(1,2,3), labels=c('NNRTI-only','PI-only','Other')) fhaart$fhaart_class_dichot <- ifelse(fhaart$fhaart_class == 2,0, ifelse(!is.na(fhaart$fhaart_class),1,NA)) fhaart$fhaart_class_dichot.factor <- factor(fhaart$fhaart_class_dichot, levels=c(0,1), labels=c('PI-only', 'non-PI')) fhaart <- subset(fhaart, select=c(unique_id, art_sd, art_id, fhaart_class, fhaart_class.factor, year_fhaart,fhaart_class_dichot, fhaart_class_dichot.factor)) names(fhaart) <- c('unique_id','date_fhaart','fhaart','fhaart_class','fhaart_class.factor','year_fhaart','fhaart_class_dichot','fhaart_class_dichot.factor') # Creating collapsed mode variable from basic data # Defining age of HIV diagnosis basic$hivdiagnosis_d <- ifelse(basic$hivdiagnosis_d == '1900-01-01',NA,basic$hivdiagnosis_d) basic$hivdiagnosis_d <- as.Date(basic$hivdiagnosis_d, origin='1970-01-01') basic$age_hiv_dx <- as.numeric(basic$hivdiagnosis_d - basic$birth_d)/365.25 basic$inclusion <- ifelse(basic$mode == 'Perinatal' | ((basic$mode == 'Unknown' | is.na(basic$mode)) & (basic$age_hiv_dx <= 10.0 | is.na(basic$age_hiv_dx))),1,0) lgth.with.perinatal <- length(which(basic$mode == 'Perinatal')) lgth.with.unk.age.le.10 <- length(which(basic$mode == 'Unknown' & basic$age_hiv_dx <= 10.0)) basic$mode.factor <- factor(basic$mode) # Defining prior ART art$date_fhaart <- fhaart[match(art$unique_id, fhaart$unique_id, nomatch=NA),'date_fhaart'] art$death_y <- follow[match(art$unique_id, follow$unique_id, nomatch=NA),'death_y'] art$death_d <- follow[match(art$unique_id, follow$unique_id, nomatch=NA),'death_d'] art$drop_y <- follow[match(art$unique_id, follow$unique_id, nomatch=NA),'drop_y'] pre_fhaart <- subset(art, (!is.na(art_sd) & !is.na(date_fhaart) & art_sd < date_fhaart) | is.na(date_fhaart)) art$prior_art <- ifelse(art$unique_id %in% pre_fhaart$unique_id,1,0) ## Note: Originally had coded those who started ART w/in 16 days of birth and were on it longer than 60 days as having a pmtct value of 99 ## We have since decided to consider these cases as PMTCT so have changed it to 1 art$pmtct <- ifelse(art$prior_art == 1 & ((art$art_sd < art$date_fhaart) | is.na(art$date_fhaart)) & art$within_16_days == 1 & art$days_on_art <= 60,1, ifelse(art$prior_art == 1 & ((art$art_sd < art$date_fhaart) | is.na(art$date_fhaart)) & art$within_16_days == 1 & art$days_on_art > 60,1,0)) # Prior ART should combine my definition with recart_y (note that recart_y = 1 only if the ART was within a certain time frame of the enrollment date) art$recart_y <- basic[match(art$unique_id, basic$unique_id, nomatch=NA),'recart_y'] art$prior_art <- ifelse(art$prior_art == 1 | art$recart_y == 1,1,0) # In case need to distinguish those who started PMTCT but had no other prior ART any_pmtct <- lapply(split(subset(art,select=c('unique_id','pmtct')),art$unique_id),FUN=function(y){ if(all(is.na(y[,'pmtct']))){ df <- data.frame('unique_id'=y[1,'unique_id'], 'any_pmtct'=NA) } else { df=data.frame('unique_id'=y[1,'unique_id'], 'any_pmtct'=max(y[,'pmtct'],na.rm=TRUE)) } return(df) }) any_pmtct_df <- do.call('rbind',any_pmtct) art$any_pmtct <- any_pmtct_df[match(art$unique_id, any_pmtct_df$unique_id, nomatch=NA),'any_pmtct'] pre_fhaart$any_pmtct <- any_pmtct_df[match(pre_fhaart$unique_id, any_pmtct_df$unique_id, nomatch=NA),'any_pmtct'] any_recarty <- lapply(split(subset(art,select=c('unique_id','recart_y')),art$unique_id),FUN=function(y){ df=data.frame('unique_id'=y[1,'unique_id'], 'any_recarty'=max(y[,'recart_y'],na.rm=TRUE)) return(df) }) any_recarty_df <- do.call('rbind',any_recarty) #-----------------------------------# # Defining baseline CD4/CD4% # #-----------------------------------# # There are some subjects with duplicate records on a given date so will average the values on those days cd4$id_date <- paste(cd4$unique_id, cd4$cd4_d,sep='.') cd4$cd4_v_new <- unlist(lapply(split(cd4[,c('cd4_v')],cd4[,c('id_date')]),FUN=function(y){ if(length(y) > 1){ tmp <- (sum(y,na.rm=TRUE)/length(y[which(!is.na(y))])) return(rep(tmp, times=length(y[which(!is.na(y))]))) } else { return(y) } })) cd4$cd4_per_new <- unlist(lapply(split(cd4[,c('cd4_per')],cd4[,c('id_date')]),FUN=function(y){ if(length(y) > 1){ tmp <- (sum(y,na.rm=TRUE)/length(y[which(!is.na(y))])) return(rep(tmp, times=length(y))) } else { return(y) } })) # Now removing those records where the id_date is duplicated cd4 <- subset(cd4, !duplicated(id_date)) cd4 <- upData(cd4, drop=c('cd4_v','cd4_per')) cd4 <- upData(cd4, rename=c('cd4_v_new'='cd4_v', 'cd4_per_new'='cd4_per')) cd4$date_fhaart <- fhaart[match(cd4$unique_id, fhaart$unique_id, nomatch=NA),'date_fhaart'] cd4$date_diff <- as.numeric(cd4$cd4_d - cd4$date_fhaart) cd4$window_ind <- ifelse(cd4$date_diff >= -180 & cd4$date_diff <= 7,1,0) baseline_cd4 <- lapply(split(cd4[,c('unique_id','window_ind','date_diff','cd4_v','cd4_per')],cd4[,'unique_id']),FUN=function(y){ tmp <- subset(y,window_ind == 1) if(dim(tmp)[1] > 0){ min_diff <- min(abs(tmp[,'date_diff']),na.rm=TRUE) which.row <- which(abs(tmp[,'date_diff']) == min_diff) b_cd4 <- tmp[which.row,'cd4_v'] b_cd4_per <- tmp[which.row,'cd4_per'] } else { b_cd4 <- b_cd4_per <- NA } df <- data.frame('unique_id'=y[1,'unique_id'], 'b_cd4'=b_cd4, 'b_cd4_per'=b_cd4_per) return(df) }) baseline_cd4_df <- do.call('rbind',baseline_cd4) #-----------------------------------# # Defining baseline HIV-1 RNA # #-----------------------------------# # Setting negative ones to cut off. Will use 399 rna$rna_v <- ifelse(abs(rna$rna_v) <= 400, 399, rna$rna_v) # There are some subjects with duplicate records on a given date so will average the values on those days rna$id_date <- paste(rna$unique_id, rna$rna_d,sep='.') rna$rna_v_new <- unlist(lapply(split(rna[,c('rna_v')],rna[,c('id_date')]),FUN=function(y){ if(length(y) > 1){ tmp <- (sum(y,na.rm=TRUE)/length(y[which(!is.na(y))])) return(rep(tmp, times=length(y[which(!is.na(y))]))) } else { return(y) } })) # Now removing those records where the id_date is duplicated rna <- subset(rna, !duplicated(id_date)) rna <- upData(rna, drop=c('rna_v')) rna <- upData(rna, rename=c('rna_v_new'='rna_v')) rna$date_fhaart <- fhaart[match(rna$unique_id, fhaart$unique_id, nomatch=NA),'date_fhaart'] rna$date_diff <- as.numeric(rna$rna_d - rna$date_fhaart) rna$window_ind <- ifelse(rna$date_diff >= -180 & rna$date_diff <= 7,1,0) baseline_rna <- lapply(split(rna[,c('unique_id','window_ind','date_diff','rna_v')],rna[,'unique_id']),FUN=function(y){ tmp <- subset(y,window_ind == 1) if(dim(tmp)[1] > 0){ min_diff <- min(abs(tmp[,'date_diff']),na.rm=TRUE) which.row <- which(abs(tmp[,'date_diff']) == min_diff) b_rna <- tmp[which.row,'rna_v'] } else { b_rna <- NA } df <- data.frame('unique_id'=y[1,'unique_id'], 'b_rna'=b_rna) return(df) }) baseline_rna_df <- do.call('rbind',baseline_rna) #-----------------------------------# # Combining data # #-----------------------------------# # Removed aids_cl_y and aids_cl_d from list on 3/31/2015 since the updated data does not have that variable in table basic peds_data <- merge(subset(basic, select=c('unique_id','birth_d','enrol_d','hivdiagnosis_d','aids_d','site','center','recart_y','aids_y','male','mode.factor', 'birth_mode','inclusion')), subset(fhaart, unique_id %in% basic$unique_id), all=TRUE) peds_data <- merge(peds_data, subset(follow, select=c('unique_id','death_d','drop_d','l_alive_d','drop_y','death_y')),all=TRUE) peds_data$drop_y <- ifelse(is.na(peds_data$drop_y),0,peds_data$drop_y) peds_data$days_from_baseline_to_fhaart <- with(peds_data, as.numeric(date_fhaart - enrol_d)) peds_data$months_from_baseline_to_fhaart <- 12*peds_data$days_from_baseline_to_fhaart/365.25 peds_data$mode.factor <- ifelse(peds_data$mode.factor == 'Transfusion nonhemophilia related','Transfusion', ifelse(peds_data$mode.factor == 'Other (specify in mode_oth)','Other',as.character(peds_data$mode.factor))) peds_data$mode.factor <- factor(peds_data$mode.factor) peds_data$any_pmtct <- any_pmtct_df[match(peds_data$unique_id, any_pmtct_df$unique_id, nomatch=NA),'any_pmtct'] peds_data$any_recart_y <- any_recarty_df[match(peds_data$unique_id, any_recarty_df$unique_id, nomatch=NA),'any_recarty'] # Age fhaart peds_data$age_fhaart <- as.numeric(peds_data$date_fhaart - peds_data$birth_d)/365.25 # Adding in baseline CD4 and CD4% values from baseline_cd4_df peds_data$baseline_cd4 <- baseline_cd4_df[match(peds_data$unique_id, baseline_cd4_df$unique_id,nomatch=NA),'b_cd4'] peds_data$baseline_cd4_per <- baseline_cd4_df[match(peds_data$unique_id, baseline_cd4_df$unique_id,nomatch=NA),'b_cd4_per'] # Adding in baseline HIV1-RNA peds_data$baseline_rna <- baseline_rna_df[match(peds_data$unique_id, baseline_rna_df$unique_id,nomatch=NA),'b_rna'] which.rows <- which(peds_data$baseline_rna > 0) peds_data$log10_baseline_rna <- NA peds_data$log10_baseline_rna[which.rows] <- log10(peds_data$baseline_rna[which.rows]) # Adding in prior ART indicator peds_data$prior_art <- art[match(peds_data$unique_id, art$unique_id, nomatch=NA),'prior_art'] peds_data$prior_art.factor <- factor(peds_data$prior_art, levels=c(0,1), labels=c('No','Yes')) ids_w_prior_art <- unique(peds_data$unique_id[which(peds_data$prior_art == 1)]) ids_w_recart_y <- unique(peds_data$unique_id[which(peds_data$recart_y == 1)]) ids_w_pmtct <- unique(peds_data$unique_id[which(peds_data$any_pmtct == 1)]) ids_in_pre_fhaart <- unique(peds_data$unique_id[which(peds_data$unique_id %in% pre_fhaart$unique_id)]) ids_in_pre_fhaart_wo_pmtct <- unique(pre_fhaart$unique_id[which(pre_fhaart$unique_id %in% peds_data$unique_id & pre_fhaart$any_pmtct == 0)]) intersect(ids_w_pmtct,ids_w_recart_y) setdiff(ids_w_pmtct,ids_w_recart_y) # Creating unique site by center variable peds_data$site_center <- ifelse(peds_data$site %in% c('brazil','honduras'),paste(peds_data$site,peds_data$center,sep='.'),peds_data$site) peds_data$site_center.factor <- factor(peds_data$site_center) # Adding in close dates by site to compute lost to follow-up: # Brazil UFMG: 2013-12-18 # Brazil UNIFESP: 2014-05-23 # Haiti: 2013-12-30 # Argentina: 2013-07-30 # Honduras: 2014-08-26 peds_data$cohort_close_date <- ifelse(peds_data$site_center.factor == 'argentina','2013-07-30', ifelse(peds_data$site_center.factor == 'brazil.UFMG','2013-12-18', ifelse(peds_data$site_center.factor == 'brazil.UNIFESP','2014-05-23', ifelse(peds_data$site_center.factor == 'haiti','2013-12-30', ifelse(peds_data$site_center.factor %in% c('honduras.Escuela','honduras.IHSS'),'2014-08-26',NA))))) peds_data$cohort_close_date <- as.Date(peds_data$cohort_close_date, format='%Y-%m-%d') #-----------------------------------# # Creating factors # #-----------------------------------# peds_data$death_y.factor <- factor(peds_data$death_y, levels=c(0,1), labels=c('Alive','Died')) peds_data$drop_y.factor <- factor(peds_data$drop_y, levels=c(0,1), labels=c('No','Yes')) peds_data$male.factor <- factor(peds_data$male, levels=c(0,1), labels=c('Female','Male')) peds_data$aids_y.factor <- factor(peds_data$aids_y, levels=c(0,1), labels=c('No','Yes')) peds_data$recart_y.factor <- factor(peds_data$recart_y, levels=c(0,1), labels=c('No','Yes')) #-----------------------------------# # Removing subjects from Peru # #-----------------------------------# lgth.ge.18 <- length(which(peds_data$age_fhaart >= 18.0 & peds_data$site_center.factor != 'peru')) lgth.ge.18.by.site <- tapply(peds_data$age_fhaart, INDEX=peds_data$site_center.factor, FUN=function(y){length(which(y >= 18.0))}) lgth.non.perinatal <- length(which(peds_data$age_fhaart[which(peds_data$site_center.factor != 'peru')] < 18.0 & peds_data$inclusion[which(peds_data$site_center.factor != 'peru')] == 0)) lgth.non.perinatal.by.site <- lapply(split(peds_data[,c('age_fhaart','inclusion')], peds_data$site_center.factor), FUN=function(y){ length(which(y[,'age_fhaart'] < 18.0 & y[,'inclusion'] == 0)) }) lgth.no.fhaart <- length(which(is.na(peds_data$age_fhaart[which(peds_data$site_center.factor != 'peru')]))) lgth.no.fhaart.by.site <- tapply(peds_data$age_fhaart, INDEX=peds_data$site_center.factor, FUN=function(y){length(which(is.na(y)))}) peds_data <- subset(peds_data, site_center.factor != 'peru' & age_fhaart < 18 & inclusion == 1) peds_data$site_center.factor <- factor(peds_data$site_center) #-----------------------------------# # Follow-up time # #-----------------------------------# peds_data$followup_time <- with(peds_data, as.numeric(l_alive_d - date_fhaart)) peds_data$followup_time_yrs <- peds_data$followup_time/365.25 #-----------------------------------# # Lost to follow-up # #-----------------------------------# # This is defined as having >= 12 months between a person's l_alive_d and the maximum last_alive_d (changed to cohort close date 10/22/2014) for a given cohort. # First need to find the maximum l_alive_d for each cohort max_l_alive_d <- aggregate(peds_data$l_alive_d, by=list(Site=peds_data$site_center.factor), FUN=max) names(max_l_alive_d)[2] <- 'max_l_alive_d' peds_data$max_l_alive_d <- max_l_alive_d[match(peds_data$site_center.factor, max_l_alive_d$Site, nomatch=NA),'max_l_alive_d'] # peds_data$months_from_l_alive_to_close <- 12*as.numeric(peds_data$max_l_alive_d - peds_data$l_alive_d)/365.25 peds_data$months_from_l_alive_to_close <- ifelse(!is.na(peds_data$death_d), 12*as.numeric(peds_data$death_d - peds_data$l_alive_d)/365.25, 12*as.numeric(peds_data$cohort_close_date - peds_data$l_alive_d)/365.25) peds_data$l2fu <- ifelse(peds_data$months_from_l_alive_to_close > 12.0,1,0) peds_data$l2fu.factor <- factor(peds_data$l2fu, levels=c(0,1), labels=c('No','Yes')) # IDs lost to follow-up per site ids_l2fu_by_site <- peds_data[which(peds_data$l2fu == 1),c('unique_id','site_center.factor','l_alive_d','max_l_alive_d')] ids_l2fu_by_site <- ids_l2fu_by_site[order(ids_l2fu_by_site$unique_id, ids_l2fu_by_site$site_center.factor),] # Removing IDs not in peds_data ids_w_recart_y <- ids_w_recart_y[which(ids_w_recart_y %in% peds_data$unique_id)] ids_w_pmtct <- ids_w_pmtct[which(ids_w_pmtct %in% peds_data$unique_id)] ids_w_prior_art <- ids_w_prior_art[which(ids_w_prior_art %in% peds_data$unique_id)] ids_in_pre_fhaart <- ids_in_pre_fhaart[which(ids_in_pre_fhaart %in% peds_data$unique_id)] ids_w_recart_y_only <- setdiff(ids_w_recart_y, c(ids_in_pre_fhaart, ids_w_pmtct)) # ------------------------------------------ # # Hard coding per Carina's request: # # Recoding recart_y to be 0 for those in # # Argentina who had recart_y_only = 1 # # ------------------------------------------ # peds_data$recart_y <- ifelse(peds_data$unique_id %in% ids_w_recart_y_only & peds_data$site == 'argentina',0,peds_data$recart_y) which.argentina <- sapply(ids_w_recart_y_only, FUN=function(y){ site <- strsplit(y, split="[.]")[[1]][1] if(site == 'ar') { return(1) } else { return(0) } }) argentina_ids_2_remove <- ids_w_recart_y_only[which(which.argentina == 1)] ids_w_recart_y_only <- ids_w_recart_y_only[-which(which.argentina == 1)] recart_y_only_by_site <- table(peds_data$site_center.factor[which(peds_data$unique_id %in% ids_w_recart_y_only)]) ids_w_art_after_enrollment <- setdiff(ids_in_pre_fhaart, c(ids_w_pmtct, ids_w_recart_y)) art_after_enrollment_by_site <- table(peds_data$site_center.factor[which(peds_data$unique_id %in% ids_w_art_after_enrollment)]) ids_w_pmtct_only <- setdiff(ids_w_pmtct, c(ids_w_recart_y_only, ids_w_art_after_enrollment)) ids_w_both_recart_y_and_prefhaart <- peds_data$unique_id[which(peds_data$unique_id %nin% c(argentina_ids_2_remove,ids_w_recart_y_only,ids_w_art_after_enrollment,ids_w_pmtct_only) & peds_data$prior_art == 1)] both_recart_y_and_prefhaart_by_site <- table(peds_data$site_center.factor[which(peds_data$unique_id %in% ids_w_both_recart_y_and_prefhaart)]) # -------------------------------------- # # Removing those whose first HAART is # # before enrol_d # # -------------------------------------- # num_w_fhaart_before_enrold <- length(which(peds_data$date_fhaart < peds_data$enrol_d)) peds_data$fhaart_before_enrol_d <- ifelse(peds_data$date_fhaart < peds_data$enrol_d,1,0) num_w_fhaart_before_enrold_by_site <- tapply(peds_data$fhaart_before_enrol_d, INDEX=peds_data$site_center.factor, FUN=function(y){sum(y)}) peds_data$pre_art_aids.clinical <- AIDSstatus[match(peds_data$unique_id, AIDSstatus$unique_id,nomatch=NA),'preARTAIDS.clinical'] peds_data$pre_art_aids.clinical <- factor(peds_data$pre_art_aids.clinical, levels=c('not AIDS','AIDS','Unknown')) # Creating indicators for missing of CD4 and CD4% to include in Table 1 peds_data$missing_cd4 <- ifelse(is.na(peds_data$baseline_cd4),1,0) peds_data$missing_cd4.factor <- factor(peds_data$missing_cd4, levels=c(0,1), labels=c('No','Yes')) peds_data$missing_viral_load <- ifelse(is.na(peds_data$baseline_rna),1,0) peds_data$missing_viral_load.factor <- factor(peds_data$missing_viral_load, levels=c(0,1), labels=c('No','Yes')) peds_data$missing_cd4_per <- ifelse(is.na(peds_data$baseline_cd4_per),1,0) peds_data$missing_cd4_per.factor <- factor(peds_data$missing_cd4_per, levels=c(0,1), labels=c('No','Yes')) peds_data$mode.factor <- factor(peds_data$mode.factor, levels=c('Perinatal','Unknown')) peds_data$year_fhaart_cat.factor <- factor(peds_data$year_fhaart, levels=seq(1997,2013)) label(peds_data$death_y.factor) <- 'Status' label(peds_data$days_from_baseline_to_fhaart) <- 'Days from enrollment to first cART' label(peds_data$months_from_baseline_to_fhaart) <- 'Months from enrollment to first cART' label(peds_data$year_fhaart) <- 'Year of first cART' label(peds_data$fhaart_class.factor) <- 'First cART class' label(peds_data$fhaart_class_dichot.factor) <- 'First cART class (dichotomized)' label(peds_data$followup_time) <- 'First cART to time to death' label(peds_data$male.factor) <- 'Sex' label(peds_data$aids_y.factor) <- 'AIDS' label(peds_data$recart_y.factor) <- 'recart_y' label(peds_data$mode.factor) <- 'Route of infection' label(peds_data$baseline_cd4) <- 'Baseline CD4' label(peds_data$baseline_cd4_per) <- 'Baseline CD4 percent' label(peds_data$log10_baseline_rna) <- 'Baseline HIV-1 RNA (log10)' label(peds_data$pre_art_aids.clinical) <- 'Clinical AIDS (created)' label(peds_data$missing_cd4.factor) <- 'Missing CD4' label(peds_data$missing_cd4_per.factor) <- 'Missing CD4 percent' label(peds_data$missing_viral_load.factor) <- 'Missing viral load' label(peds_data$site_center.factor) <- 'Site' label(peds_data$l2fu.factor) <- 'Lost to follow-up' label(peds_data$followup_time_yrs) <- 'Follow-up time (yrs)' label(peds_data$prior_art.factor) <- 'Prior ART (including PMTCT)' label(peds_data$year_fhaart_cat.factor) <- 'Year of first cART categorized' #-----------------------------------# # Subsets based on age # #-----------------------------------# peds_data$log10_followup_time_yrs <- log10(peds_data$followup_time_yrs) peds_data$site.factor <- factor(peds_data$site) peds_data$alternate_site_center <- ifelse(peds_data$site_center %in% c('honduras.Escuela','honduras.IHSS'),'honduras',peds_data$site_center) peds_data$alternate_site_center.factor <- factor(peds_data$alternate_site_center) # Year of enrollment peds_data$yr_enrollment <- apply(subset(peds_data, select='enrol_d'),1,FUN=function(y){ as.numeric(unlist(strsplit(as.character(y),split='-')[[1]][1])) }) label(peds_data$yr_enrollment) <- 'Year of enrollment' # Adding in new PMTCT variable peds_data$pmtct <- basic[match(peds_data$unique_id, basic$unique_id, nomatch=NA),'pmtct'] peds_data$pmtct.factor <- factor(peds_data$pmtct, levels=c(0,1), labels=c('No','Yes')) label(peds_data$pmtct.factor) <- 'PMTCT' peds_data$any_pmtct.factor <- factor(peds_data$any_pmtct, levels=c(0,1), labels=c('No','Yes')) label(peds_data$any_pmtct.factor) <- 'PMTCT (estimated from ART data)' # Subsetting out those subjects with 0 follow-up time lgth.w.0.fu.time <- length(which(peds_data$followup_time_yrs == 0.0)) lgth.w.0.fu.time.by.site <- tapply(peds_data$followup_time_yrs, INDEX=peds_data$site_center.factor, FUN=function(y){length(which(y == 0.0))}) no_fu_time <- subset(peds_data, followup_time_yrs == 0.0) peds_data <- subset(peds_data, followup_time_yrs > 0.0) # Need to remove subjects who had prior ART or recart_y # Need to update ids_w_recart_y_only since one was apparently removed above ids_w_recart_y_only <- ids_w_recart_y_only[which(ids_w_recart_y_only %in% peds_data$unique_id)] peds_data <- subset(peds_data, unique_id %nin% c(ids_w_recart_y_only, ids_w_art_after_enrollment, ids_w_both_recart_y_and_prefhaart)) peds_data <- subset(peds_data, fhaart_before_enrol_d == 0) tmp_le_5 <- subset(peds_data, age_fhaart < 5) tmp_gt_5 <- subset(peds_data, age_fhaart >= 5) art$l2fu <- basic[match(art$unique_id, art$l2fu, nomatch=NA),'l2fu'] prior_art_df <- data.frame('unique_id'=c(ids_w_recart_y_only, ids_w_art_after_enrollment, ids_w_pmtct_only, ids_w_both_recart_y_and_prefhaart), 'type_of_prior_art'=c(rep('recart y only',times=length(ids_w_recart_y_only)), rep('ART after enrollment',times=length(ids_w_art_after_enrollment)), rep('PMTCT only',times=length(ids_w_pmtct_only)), rep('Both recart y and ART after enrollment',times=length(ids_w_both_recart_y_and_prefhaart)))) # A table of subjects who will be excluded due to recart_y = 1 or prior_art = 1 tbl0 <- data.frame('Reason'=c('Age >= 18','Non-perinatal mode of transmission','No first cART','Prior ART designated by recart y only', 'Prior ART with ART after enrollment only','Prior ART after enrollment and designated with recart y','No follow-up time','First cART prior to enrollment'), 'Frequency'=c(lgth.ge.18,lgth.non.perinatal,lgth.no.fhaart,length(ids_w_recart_y_only),length(ids_w_art_after_enrollment),length(ids_w_both_recart_y_and_prefhaart),lgth.w.0.fu.time, num_w_fhaart_before_enrold),check.names=FALSE) # By site per Carina's request: tbl0_by_site <- data.frame('Reason'= c('Age >= 18','Non-perinatal mode of transmission','No first cART','Prior ART designated by recart y only', 'Prior ART with ART after enrollment only','Prior ART after enrollment and designated with recart y','No follow-up time','First cART prior to enrollment'), 'Argentina'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'argentina')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'argentina')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'argentina')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'argentina')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'argentina')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'argentina')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'argentina')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'argentina')]), 'Brazil - UFMG'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'brazil.UFMG')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'brazil.UFMG')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'brazil.UFMG')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'brazil.UFMG')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'brazil.UFMG')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'brazil.UFMG')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'brazil.UFMG')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'brazil.UFMG')]), 'Brazil - UNIFESP'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'brazil.UNIFESP')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'brazil.UNIFESP')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'brazil.UNIFESP')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'brazil.UNIFESP')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'brazil.UNIFESP')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'brazil.UNIFESP')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'brazil.UNIFESP')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'brazil.UNIFESP')]), 'Haiti'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'haiti')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'haiti')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'haiti')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'haiti')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'haiti')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'haiti')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'haiti')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'haiti')]), 'Honduras - Escuela'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'honduras.Escuela')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'honduras.Escuela')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'honduras.Escuela')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'honduras.Escuela')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'honduras.Escuela')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'honduras.Escuela')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'honduras.Escuela')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'honduras.Escuela')]), 'Honduras - IHSS'=c(lgth.ge.18.by.site[which(names(lgth.ge.18.by.site) == 'honduras.IHSS')], unlist(lgth.non.perinatal.by.site[which(names(lgth.non.perinatal.by.site) == 'honduras.IHSS')]), lgth.no.fhaart.by.site[which(names(lgth.no.fhaart.by.site) == 'honduras.IHSS')], recart_y_only_by_site[which(names(recart_y_only_by_site) == 'honduras.IHSS')], art_after_enrollment_by_site[which(names(art_after_enrollment_by_site) == 'honduras.IHSS')], both_recart_y_and_prefhaart_by_site[which(names(both_recart_y_and_prefhaart_by_site) == 'honduras.IHSS')], lgth.w.0.fu.time.by.site[which(names(lgth.w.0.fu.time.by.site) == 'honduras.IHSS')], num_w_fhaart_before_enrold_by_site[which(names(num_w_fhaart_before_enrold_by_site) == 'honduras.IHSS')]), check.names=FALSE) # Table 1 by site label(peds_data$age_fhaart) <- 'Age at cART initiation' tbl1 <- summary(site_center.factor ~ age_fhaart + death_y.factor + l2fu.factor + months_from_baseline_to_fhaart + year_fhaart + year_fhaart_cat.factor + fhaart_class.factor + fhaart_class_dichot.factor + followup_time_yrs + male.factor + mode.factor + pmtct.factor + any_pmtct.factor + prior_art.factor + baseline_cd4 + missing_cd4.factor + baseline_cd4_per + missing_cd4_per.factor + log10_baseline_rna + missing_viral_load.factor + pre_art_aids.clinical + yr_enrollment, data=peds_data, method='reverse',overall=TRUE) tbl1_le5 <- summary(site_center.factor ~ age_fhaart + death_y.factor + l2fu.factor + months_from_baseline_to_fhaart + year_fhaart + year_fhaart_cat.factor + fhaart_class.factor + fhaart_class_dichot.factor + followup_time_yrs + male.factor + mode.factor + any_pmtct.factor + pmtct.factor + prior_art.factor + baseline_cd4 + missing_cd4.factor + baseline_cd4_per + missing_cd4_per.factor + log10_baseline_rna + missing_viral_load.factor + pre_art_aids.clinical + yr_enrollment, data=tmp_le_5,overall=TRUE, method='reverse') tbl1_gt5 <- summary(site_center.factor ~ age_fhaart + death_y.factor + l2fu.factor + months_from_baseline_to_fhaart + year_fhaart + year_fhaart_cat.factor + fhaart_class.factor + fhaart_class_dichot.factor + followup_time_yrs + male.factor + mode.factor + any_pmtct.factor + pmtct.factor + prior_art.factor + baseline_cd4 + missing_cd4.factor + baseline_cd4_per + missing_cd4_per.factor + log10_baseline_rna + missing_viral_load.factor + pre_art_aids.clinical + yr_enrollment, data=tmp_gt_5,overall=TRUE,method='reverse') # Reporting the first regimens for subjects and their subsequent regimens art <- art[order(art$unique_id, art$art_sd),] first_reg <- art[!duplicated(art$unique_id),c('unique_id','art_id','art_class')] names(first_reg) <- c('unique_id','first_regimen','first_regimen_class') art$first_regimen <- first_reg[match(art$unique_id, first_reg$unique_id,nomatch=NA),'first_regimen'] art$first_regimen_class <- first_reg[match(art$unique_id, first_reg$unique_id,nomatch=NA),'first_regimen_class'] sub_art <- subset(art, first_regimen_class == 'non-haart') unique_first_reg <- unique(sub_art$first_regimen) subsequent_regimens <- vector('list', length(unique_first_reg)) names(subsequent_regimens) <- unique_first_reg for(i in 1:length(unique_first_reg)){ tmp <- subset(sub_art, first_regimen == unique_first_reg[i],select=c('unique_id','art_id','art_class','first_regimen')) tmp$first_record <- 0 tmp$first_record[!duplicated(tmp$unique_id)] <- 1 tmp <- subset(tmp, first_record == 0) subsequent_regimens[[i]] <- unique(tmp$art_id) } regimen_df <- data.frame('first regimen'=unique_first_reg, 'subsequent regimens'=sapply(subsequent_regimens, FUN=function(y){paste(y,collapse='; ')},simplify=TRUE),check.names=FALSE) cum_inc_le_5 <- cuminc(ftime=tmp_le_5$followup_time_yrs, fstatus=tmp_le_5$death_y, cencode='0',rho=0) cum_inc_gt_5 <- cuminc(ftime=tmp_gt_5$followup_time_yrs, fstatus=tmp_gt_5$death_y, cencode='0',rho=0) plot(cum_inc_le_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)] ~ cum_inc_le_5[[1]]$time[which(cum_inc_le_5[[1]]$time <= 10.0)], xlim=c(0,10), type='l', axes=FALSE, xlab='Years from cART initiation',ylab='Cumulative incidence',lty=1,col='black', ylim=range(c(cum_inc_le_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)], cum_inc_gt_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)]),na.rm=TRUE),lwd=1.5) lines(cum_inc_gt_5[[1]]$est[which(cum_inc_gt_5[[1]]$time <= 10.0)] ~ cum_inc_gt_5[[1]]$time[which(cum_inc_gt_5[[1]]$time <= 10.0)], lty=1,col='grey',lwd=1.5) axis(1, at=c(0,2,4,6,8,10), labels=c(0,2,4,6,8,10)) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.95) # Need overall cuminc cum_inc_overall <- cuminc(ftime=peds_data$followup_time_yrs, fstatus=peds_data$death_y, cencode='0',rho=0) which.1 <- max(which(ceiling(cum_inc_overall[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_overall[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_overall[[1]]$time) == 5.0)) var_overall <- cum_inc_overall[[1]]$var inc1_overall <- cum_inc_overall[[1]]$est[which.1] lci1_overall <- inc1_overall - qnorm(0.975)*sqrt(var_overall[which.1]) uci1_overall <- inc1_overall + qnorm(0.975)*sqrt(var_overall[which.1]) inc3_overall <- cum_inc_overall[[1]]$est[which.3] lci3_overall <- inc3_overall - qnorm(0.975)*sqrt(var_overall[which.3]) uci3_overall <- inc3_overall + qnorm(0.975)*sqrt(var_overall[which.3]) inc5_overall <- cum_inc_overall[[1]]$est[which.5] lci5_overall <- inc5_overall - qnorm(0.975)*sqrt(var_overall[which.5]) uci5_overall <- inc5_overall + qnorm(0.975)*sqrt(var_overall[which.5]) # For those < 5 var_le5 <- cum_inc_le_5[[1]]$var which.1 <- max(which(ceiling(cum_inc_le_5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_le_5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_le_5[[1]]$time) == 5.0)) inc1_le5 <- cum_inc_le_5[[1]]$est[which.1] lci1_le5 <- inc1_le5 - qnorm(0.975)*sqrt(var_le5[which.1]) uci1_le5 <- inc1_le5 + qnorm(0.975)*sqrt(var_le5[which.1]) inc3_le5 <- cum_inc_le_5[[1]]$est[which.3] lci3_le5 <- inc3_le5 - qnorm(0.975)*sqrt(var_le5[which.3]) uci3_le5 <- inc3_le5 + qnorm(0.975)*sqrt(var_le5[which.3]) inc5_le5 <- cum_inc_le_5[[1]]$est[which.5] lci5_le5 <- inc5_le5 - qnorm(0.975)*sqrt(var_le5[which.5]) uci5_le5 <- inc5_le5 + qnorm(0.975)*sqrt(var_le5[which.5]) var5_le5 <- var_le5[which.5] # For those >= 5 var_gt5 <- cum_inc_gt_5[[1]]$var which.1 <- max(which(ceiling(cum_inc_gt_5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_gt_5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_gt_5[[1]]$time) == 5.0)) inc1_gt5 <- cum_inc_gt_5[[1]]$est[which.1] lci1_gt5 <- inc1_gt5 - qnorm(0.975)*sqrt(var_gt5[which.1]) uci1_gt5 <- inc1_gt5 + qnorm(0.975)*sqrt(var_gt5[which.1]) inc3_gt5 <- cum_inc_gt_5[[1]]$est[which.3] lci3_gt5 <- inc3_gt5 - qnorm(0.975)*sqrt(var_gt5[which.3]) uci3_gt5 <- inc3_gt5 + qnorm(0.975)*sqrt(var_gt5[which.3]) inc5_gt5 <- cum_inc_gt_5[[1]]$est[which.5] lci5_gt5 <- inc5_gt5 - qnorm(0.975)*sqrt(var_gt5[which.5]) uci5_gt5 <- inc5_gt5 + qnorm(0.975)*sqrt(var_gt5[which.5]) var5_gt5 <- var_gt5[which.5] # Computing p-value comparing 5-yr cumulative incidence of death between those < 5 and those >= 5 z5_dth <- (inc5_le5 - inc5_gt5)/sqrt(var5_le5 + var5_gt5) p5_dth <- 2*(1-pnorm(z5_dth)) # Overall by site_center cum_inc_brufmg <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'brazil.UFMG')], fstatus=peds_data$death_y[which(peds_data$site_center.factor == 'brazil.UFMG')], cencode='0',rho=0) var_brufmg <- cum_inc_brufmg[[1]]$var which.1 <- max(which(ceiling(cum_inc_brufmg[[1]]$time) == 1.0)) which.3 <- max(which(round(cum_inc_brufmg[[1]]$time) == 3.0)) which.5 <- max(which(round(cum_inc_brufmg[[1]]$time) == 5.0)) inc1_brufmg <- cum_inc_brufmg[[1]]$est[which.1] lci1_brufmg <- inc1_brufmg - qnorm(0.975)*sqrt(var_brufmg[which.1]) uci1_brufmg <- inc1_brufmg + qnorm(0.975)*sqrt(var_brufmg[which.1]) inc3_brufmg <- cum_inc_brufmg[[1]]$est[which.3] lci3_brufmg <- inc3_brufmg - qnorm(0.975)*sqrt(var_brufmg[which.3]) uci3_brufmg <- inc3_brufmg + qnorm(0.975)*sqrt(var_brufmg[which.3]) inc5_brufmg <- cum_inc_brufmg[[1]]$est[which.5] lci5_brufmg <- inc5_brufmg - qnorm(0.975)*sqrt(var_brufmg[which.5]) uci5_brufmg <- inc5_brufmg + qnorm(0.975)*sqrt(var_brufmg[which.5]) cum_inc_brunifesp <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'brazil.UNIFESP')], fstatus=peds_data$death_y[which(peds_data$site_center.factor == 'brazil.UNIFESP')], cencode='0',rho=0) var_brunifesp <- cum_inc_brunifesp[[1]]$var which.1 <- max(which(ceiling(cum_inc_brunifesp[[1]]$time) == 1.0)) which.3 <- min(which(cum_inc_brunifesp[[1]]$time > 3.0)) which.5 <- min(which(cum_inc_brunifesp[[1]]$time >= 5.0)) inc1_brunifesp <- cum_inc_brunifesp[[1]]$est[which.1] lci1_brunifesp <- inc1_brunifesp - qnorm(0.975)*sqrt(var_brunifesp[which.1]) uci1_brunifesp <- inc1_brunifesp + qnorm(0.975)*sqrt(var_brunifesp[which.1]) inc3_brunifesp <- cum_inc_brunifesp[[1]]$est[which.3] lci3_brunifesp <- inc3_brunifesp - qnorm(0.975)*sqrt(var_brunifesp[which.3]) uci3_brunifesp <- inc3_brunifesp + qnorm(0.975)*sqrt(var_brunifesp[which.3]) inc5_brunifesp <- cum_inc_brunifesp[[1]]$est[which.5] lci5_brunifesp <- inc5_brunifesp - qnorm(0.975)*sqrt(var_brunifesp[which.5]) uci5_brunifesp <- inc5_brunifesp + qnorm(0.975)*sqrt(var_brunifesp[which.5]) cum_inc_haiti <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'haiti')], fstatus=peds_data$death_y[which(peds_data$site_center.factor == 'haiti')], cencode='0',rho=0) var_haiti <- cum_inc_haiti[[1]]$var which.1 <- max(which(ceiling(cum_inc_haiti[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_haiti[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_haiti[[1]]$time) == 5.0)) inc1_haiti <- cum_inc_haiti[[1]]$est[which.1] lci1_haiti <- inc1_haiti - qnorm(0.975)*sqrt(var_haiti[which.1]) uci1_haiti <- inc1_haiti + qnorm(0.975)*sqrt(var_haiti[which.1]) inc3_haiti <- cum_inc_haiti[[1]]$est[which.3] lci3_haiti <- inc3_haiti - qnorm(0.975)*sqrt(var_haiti[which.3]) uci3_haiti <- inc3_haiti + qnorm(0.975)*sqrt(var_haiti[which.3]) inc5_haiti <- cum_inc_haiti[[1]]$est[which.5] lci5_haiti <- inc5_haiti - qnorm(0.975)*sqrt(var_haiti[which.5]) uci5_haiti <- inc5_haiti + qnorm(0.975)*sqrt(var_haiti[which.5]) cum_inc_honescuela <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'honduras.Escuela')], fstatus=peds_data$death_y[which(peds_data$site_center.factor == 'honduras.Escuela')], cencode='0',rho=0) var_honescuela <- cum_inc_honescuela[[1]]$var which.1 <- max(which(ceiling(cum_inc_honescuela[[1]]$time) == 1.0)) which.3 <- min(which(cum_inc_honescuela[[1]]$time >= 3.0)) which.5 <- max(which(ceiling(cum_inc_honescuela[[1]]$time) == 5.0)) inc1_honescuela <- cum_inc_honescuela[[1]]$est[which.1] lci1_honescuela <- inc1_honescuela - qnorm(0.975)*sqrt(var_honescuela[which.1]) uci1_honescuela <- inc1_honescuela + qnorm(0.975)*sqrt(var_honescuela[which.1]) inc3_honescuela <- cum_inc_honescuela[[1]]$est[which.3] lci3_honescuela <- inc3_honescuela - qnorm(0.975)*sqrt(var_honescuela[which.3]) uci3_honescuela <- inc3_honescuela + qnorm(0.975)*sqrt(var_honescuela[which.3]) inc5_honescuela <- cum_inc_honescuela[[1]]$est[which.5] lci5_honescuela <- inc5_honescuela - qnorm(0.975)*sqrt(var_honescuela[which.5]) uci5_honescuela <- inc5_honescuela + qnorm(0.975)*sqrt(var_honescuela[which.5]) # Need to create a 3-level censoring variable that is 0=No event, 1=LTFU, 2=Death # For those < 5 tmp_le_5$l2fu_censor <- ifelse(tmp_le_5$l2fu == 1,1, ifelse(tmp_le_5$death_y == 1,2, ifelse(!is.na(tmp_le_5$l2fu) & !is.na(tmp_le_5$death_y),0,NA))) cum_inc_l2fu_le5 <- cuminc(ftime=tmp_le_5$followup_time_yrs, fstatus=tmp_le_5$l2fu_censor, cencode='0',rho=0) # For those >= 5 years old: tmp_gt_5$l2fu_censor <- ifelse(tmp_gt_5$l2fu == 1,1, ifelse(tmp_gt_5$death_y == 1,2, ifelse(!is.na(tmp_gt_5$l2fu) & !is.na(tmp_gt_5$death_y),0,NA))) cum_inc_l2fu_gt5 <- cuminc(ftime=tmp_gt_5$followup_time_yrs, fstatus=tmp_gt_5$l2fu_censor, cencode='0',rho=0) cum_inc_l2fu_overall <- cuminc(ftime=peds_data$followup_time_yrs, fstatus=peds_data$l2fu, cencode='0',rho=0) which.1 <- max(which(ceiling(cum_inc_l2fu_overall[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_l2fu_overall[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_l2fu_overall[[1]]$time) == 5.0)) which.7 <- max(which(ceiling(cum_inc_l2fu_overall[[1]]$time) == 7.0)) var_l2fu_overall <- cum_inc_l2fu_overall[[1]]$var inc1_l2fu_overall <- cum_inc_l2fu_overall[[1]]$est[which.1] lci1_l2fu_overall <- inc1_l2fu_overall - qnorm(0.975)*sqrt(var_l2fu_overall[which.1]) uci1_l2fu_overall <- inc1_l2fu_overall + qnorm(0.975)*sqrt(var_l2fu_overall[which.1]) inc3_l2fu_overall <- cum_inc_l2fu_overall[[1]]$est[which.3] lci3_l2fu_overall <- inc3_l2fu_overall - qnorm(0.975)*sqrt(var_l2fu_overall[which.3]) uci3_l2fu_overall <- inc3_l2fu_overall + qnorm(0.975)*sqrt(var_l2fu_overall[which.3]) inc5_l2fu_overall <- cum_inc_l2fu_overall[[1]]$est[which.5] lci5_l2fu_overall <- inc5_l2fu_overall - qnorm(0.975)*sqrt(var_l2fu_overall[which.5]) uci5_l2fu_overall <- inc5_l2fu_overall + qnorm(0.975)*sqrt(var_l2fu_overall[which.5]) inc7_l2fu_overall <- cum_inc_l2fu_overall[[1]]$est[which.7] lci7_l2fu_overall <- inc7_l2fu_overall - qnorm(0.975)*sqrt(var_l2fu_overall[which.7]) uci7_l2fu_overall <- inc7_l2fu_overall + qnorm(0.975)*sqrt(var_l2fu_overall[which.7]) # For those < 5 var_l2fu_le5 <- cum_inc_l2fu_le5[[1]]$var which.1 <- max(which(ceiling(cum_inc_l2fu_le5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_l2fu_le5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_l2fu_le5[[1]]$time) == 5.0)) which.7 <- max(which(ceiling(cum_inc_l2fu_le5[[1]]$time) == 7.0)) inc1_l2fu_le5 <- cum_inc_l2fu_le5[[1]]$est[which.1] lci1_l2fu_le5 <- inc1_l2fu_le5 - qnorm(0.975)*sqrt(var_l2fu_le5[which.1]) uci1_l2fu_le5 <- inc1_l2fu_le5 + qnorm(0.975)*sqrt(var_l2fu_le5[which.1]) inc3_l2fu_le5 <- cum_inc_l2fu_le5[[1]]$est[which.3] lci3_l2fu_le5 <- inc3_l2fu_le5 - qnorm(0.975)*sqrt(var_l2fu_le5[which.3]) uci3_l2fu_le5 <- inc3_l2fu_le5 + qnorm(0.975)*sqrt(var_l2fu_le5[which.3]) inc5_l2fu_le5 <- cum_inc_l2fu_le5[[1]]$est[which.5] lci5_l2fu_le5 <- inc5_l2fu_le5 - qnorm(0.975)*sqrt(var_l2fu_le5[which.5]) uci5_l2fu_le5 <- inc5_l2fu_le5 + qnorm(0.975)*sqrt(var_l2fu_le5[which.5]) inc7_l2fu_le5 <- cum_inc_l2fu_le5[[1]]$est[which.7] lci7_l2fu_le5 <- inc7_l2fu_le5 - qnorm(0.975)*sqrt(var_l2fu_le5[which.7]) uci7_l2fu_le5 <- inc7_l2fu_le5 + qnorm(0.975)*sqrt(var_l2fu_le5[which.7]) var5_l2fu_le5 <- var_l2fu_le5[which.5] # For those >= 5 var_l2fu_gt5 <- cum_inc_l2fu_gt5[[1]]$var which.1 <- max(which(ceiling(cum_inc_l2fu_gt5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_l2fu_gt5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_l2fu_gt5[[1]]$time) == 5.0)) which.7 <- max(which(ceiling(cum_inc_l2fu_gt5[[1]]$time) == 7.0)) inc1_l2fu_gt5 <- cum_inc_l2fu_gt5[[1]]$est[which.1] lci1_l2fu_gt5 <- inc1_l2fu_gt5 - qnorm(0.975)*sqrt(var_l2fu_gt5[which.1]) uci1_l2fu_gt5 <- inc1_l2fu_gt5 + qnorm(0.975)*sqrt(var_l2fu_gt5[which.1]) inc3_l2fu_gt5 <- cum_inc_l2fu_gt5[[1]]$est[which.3] lci3_l2fu_gt5 <- inc3_l2fu_gt5 - qnorm(0.975)*sqrt(var_l2fu_gt5[which.3]) uci3_l2fu_gt5 <- inc3_l2fu_gt5 + qnorm(0.975)*sqrt(var_l2fu_gt5[which.3]) inc5_l2fu_gt5 <- cum_inc_l2fu_gt5[[1]]$est[which.5] lci5_l2fu_gt5 <- inc5_l2fu_gt5 - qnorm(0.975)*sqrt(var_l2fu_gt5[which.5]) uci5_l2fu_gt5 <- inc5_l2fu_gt5 + qnorm(0.975)*sqrt(var_l2fu_gt5[which.5]) var5_l2fu_gt5 <- var_l2fu_gt5[which.5] inc7_l2fu_gt5 <- cum_inc_l2fu_gt5[[1]]$est[which.7] lci7_l2fu_gt5 <- inc7_l2fu_gt5 - qnorm(0.975)*sqrt(var_l2fu_gt5[which.7]) uci7_l2fu_gt5 <- inc7_l2fu_gt5 + qnorm(0.975)*sqrt(var_l2fu_gt5[which.7]) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) rge <- range(c(cum_inc_l2fu_le5[[1]]$est, cum_inc_l2fu_gt5[[1]]$est)) plot(cum_inc_l2fu_le5[[1]]$est ~ cum_inc_l2fu_le5[[1]]$time, type='l',axes=FALSE, xlab='Years from cART initiation',ylab='Cumulative incidence of loss to follow-up',ylim=rge,col='black',lwd=1.5) lines(cum_inc_l2fu_gt5[[1]]$est ~ cum_inc_l2fu_gt5[[1]]$time,col='grey',lwd=1.5) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_le_5 <- Surv(tmp_le_5$followup_time_yrs,tmp_le_5$death_y) Srv_gt_5 <- Surv(tmp_gt_5$followup_time_yrs,tmp_gt_5$death_y) peds_data <- upData(peds_data, drop=c('birth_mode','prior_art','prior_art.factor'))## Added 1/9/15 dd <- datadist(peds_data) options(datadist='dd') tmp_le_5$baseline_cd4_sqrt <- sqrt(tmp_le_5$baseline_cd4) var_list_le_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, pre_art_aids.clinical, strat(site_center.factor)) var_list_gt_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, pre_art_aids.clinical, strat(site_center.factor)) cph.le5.fmla <- as.formula(paste('Srv_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.gt5.fmla <- as.formula(paste('Srv_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# tmp_le_5 <- upData(tmp_le_5, drop=c('birth_mode')) ## Added 1/9/15 dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_impute <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) mod_le_5 <- fit.mult.impute(cph.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------------------------------# # Another imputation method for comparison # #----------------------------------------------------# set.seed(2) impute_mod_cd4 <- glm(baseline_cd4_sqrt ~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, family='gaussian') newdata.cd4 <- data.frame('log10_followup_time_yrs'=tmp_le_5$log10_followup_time_yrs, 'death_y'=tmp_le_5$death_y, 'age_fhaart'=tmp_le_5$age_fhaart, 'male.factor'=tmp_le_5$male.factor, 'year_fhaart'=tmp_le_5$year_fhaart, 'fhaart_class_dichot.factor'=tmp_le_5$fhaart_class_dichot.factor, 'pre_art_aids.clinical'=tmp_le_5$pre_art_aids.clinical, 'site_center.factor'=tmp_le_5$site_center.factor) p.cd4 <- predict(impute_mod_cd4, newdata=newdata.cd4, na.action=na.omit) cd4.na <- length(which(is.na(tmp_le_5$baseline_cd4_sqrt))) r.cd4 <- residuals(impute_mod_cd4) mod_le_5_personal_alt <- impute.fun(p.cd4,r.cd4,tmp_le_5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5 <- upData(tmp_gt_5, drop=c('birth_mode')) ## Added 1/9/15 dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_impute <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) mod_gt_5 <- fit.mult.impute(cph.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) le_5_df <- output_mgmt_no_hiv(mod_le_5,le5_ind=1) gt_5_df <- output_mgmt_no_hiv(mod_gt_5,le5_ind=0) # Combining < 5 and >= 5 mortality Cox results dth_df_combined <- cbind(le_5_df, gt_5_df[,-1]) # Splitting up art_id art$date_fhaart <- fhaart[match(art$unique_id, fhaart$unique_id,nomatch=NA),'date_fhaart'] art_after_fhaart <- subset(art, art_sd >= date_fhaart) # Numbering regimens to see when switch occurs art_after_fhaart$regimen_number <- NA split(art_after_fhaart$regimen_number, art_after_fhaart$unique_id) <- lapply(split(art_after_fhaart$unique_id, art_after_fhaart$unique_id),FUN=function(y){ seq(1,length(y))}) druglist <- strsplit(art_after_fhaart$art_id,split=',') drugnames <- unique(unlist(druglist)) n.drugs <- length(unique(unlist(druglist))) x <- matrix(0,nrow=nrow(art_after_fhaart),ncol=n.drugs) colnames(x) <- paste('ind_',unique(unlist(druglist)),sep='') for(i in seq_along(druglist)){ x[i,] <- as.numeric(drugnames %in% druglist[[i]]) } # A switch can be the addition of 2 or more ARVs to the initial cART regimen art_after_fhaart$numdrugs_initial_cART <- ifelse(art_after_fhaart$art_sd == art_after_fhaart$date_fhaart,art_after_fhaart$numdrugs,NA) tmp <- subset(art_after_fhaart, !duplicated(unique_id) & !is.na(date_fhaart), select=c(unique_id, numdrugs_initial_cART)) art_after_fhaart$numdrugs_initial_cART <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id,nomatch=NA),'numdrugs_initial_cART'] tmp <- subset(art_after_fhaart, !duplicated(unique_id) & !is.na(date_fhaart), select=c(unique_id, art_id)) art_after_fhaart$fhaart <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'art_id'] art_after_fhaart$potential_switch <- ifelse(art_after_fhaart$art_id == art_after_fhaart$fhaart,0, ifelse(!is.na(art_after_fhaart$fhaart) & art_after_fhaart$art_id != art_after_fhaart$fhaart,1,NA)) art_after_fhaart$diff_num_drugs <- art_after_fhaart$numdrugs - art_after_fhaart$numdrugs_initial_cART art_after_fhaart$length_drug_substitutions <- apply(art_after_fhaart[,c('diff_num_drugs','art_id','fhaart')],1,FUN=function(y){ length(setdiff(unlist(strsplit(y['art_id'],split=',')),unlist(strsplit(y['fhaart'],split=',')))) }) art_after_fhaart$number_arvs_added <- apply(art_after_fhaart[,c('art_id','fhaart')],1,FUN=function(y){ length(setdiff(unlist(strsplit(y['art_id'],split=',')),unlist(strsplit(y['fhaart'],split=',')))) }) art_after_fhaart$number_arvs_dropped <- apply(art_after_fhaart[,c('art_id','fhaart')],1,FUN=function(y){ length(setdiff(unlist(strsplit(y['fhaart'],split=',')),unlist(strsplit(y['art_id'],split=',')))) }) art_after_fhaart$net_arv_change <- with(art_after_fhaart, number_arvs_added - number_arvs_dropped) #---------------------------# # Definition 1 # #---------------------------# art_after_fhaart$switch_def_1a <- ifelse(art_after_fhaart$diff_num_drugs == 0 & art_after_fhaart$length_drug_substitutions >= 2,1,0) art_after_fhaart$switch_def_1 <- ifelse(art_after_fhaart$net_arv_change == 0 & art_after_fhaart$number_arvs_added >= 2,1,0) art_after_fhaart$switch_def_1_date <- ifelse(art_after_fhaart$switch_def_1 == 1, art_after_fhaart$art_sd,NA) tmp <- subset(art_after_fhaart, !is.na(switch_def_1_date), select=c(unique_id, switch_def_1_date)) tmp <- tmp[!duplicated(tmp$unique_id),] art_after_fhaart$switch_def_1_date <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'switch_def_1_date'] art_after_fhaart$switch_def_1_date <- as.Date(art_after_fhaart$switch_def_1_date, origin='1970-01-01') #---------------------------# # Definition 2 # #---------------------------# art_after_fhaart$switch_def_2a <- ifelse(art_after_fhaart$diff_num_drugs >= 2,1,0) art_after_fhaart$switch_def_2 <- ifelse(art_after_fhaart$number_arvs_added >= 2,1,0) art_after_fhaart$switch_def_2_date <- ifelse(art_after_fhaart$switch_def_2 == 1, art_after_fhaart$art_sd,NA) tmp <- subset(art_after_fhaart, !is.na(switch_def_2_date), select=c(unique_id, switch_def_2_date)) tmp <- tmp[!duplicated(tmp$unique_id),] art_after_fhaart$switch_def_2_date <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'switch_def_2_date'] art_after_fhaart$switch_def_2_date <- as.Date(art_after_fhaart$switch_def_2_date, origin='1970-01-01') #---------------------------# # Definition 3 # #---------------------------# art_after_fhaart$switch_def_3a <- ifelse(art_after_fhaart$net_arv_change <= -1 & art_after_fhaart$art_class == 'non-haart',1,0) art_after_fhaart$switch_def_3 <- ifelse(art_after_fhaart$number_arvs_dropped >= 1 & art_after_fhaart$art_class == 'non-haart',1,0) art_after_fhaart$switch_def_3_date <- ifelse(art_after_fhaart$switch_def_3 == 1, art_after_fhaart$art_sd,NA) tmp <- subset(art_after_fhaart, !is.na(switch_def_3_date), select=c(unique_id, switch_def_3_date)) tmp <- tmp[!duplicated(tmp$unique_id),] art_after_fhaart$switch_def_3_date <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'switch_def_3_date'] art_after_fhaart$switch_def_3_date <- as.Date(art_after_fhaart$switch_def_3_date, origin='1970-01-01') #---------------------------# # Definition 4 # #---------------------------# art_after_fhaart$single_drug_suba <- ifelse(art_after_fhaart$diff_num_drugs == 0 & art_after_fhaart$length_drug_substitutions == 1,1,0) art_after_fhaart$single_drug_sub <- ifelse(art_after_fhaart$net_arv_change == 0 & art_after_fhaart$number_arvs_added == 1,1,0) art_after_fhaart$single_drug_sub_date <- ifelse(art_after_fhaart$single_drug_sub == 1, art_after_fhaart$art_sd,NA) tmp <- subset(art_after_fhaart, !is.na(single_drug_sub_date), select=c(unique_id, single_drug_sub_date)) tmp <- tmp[!duplicated(tmp$unique_id),] art_after_fhaart$single_drug_sub_date <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'single_drug_sub_date'] art_after_fhaart$single_drug_sub_date <- as.Date(art_after_fhaart$single_drug_sub_date, origin='1970-01-01') # Need to determine if the single drug substitution is in our out of class sds <- lapply(split(subset(art_after_fhaart, select=c(unique_id,single_drug_sub_date, art_sd, art_id, fhaart)),art_after_fhaart$unique_id),FUN=function(y){ which.row <- which(!is.na(y[,'single_drug_sub_date']) & y[,'single_drug_sub_date'] == y[,'art_sd']) if(length(which.row) > 0){ old_drug <- setdiff(unlist(strsplit(y[which.row,'fhaart'],split=',')),unlist(strsplit(y[which.row,'art_id'],split=','))) new_drug <- setdiff(unlist(strsplit(y[which.row,'art_id'],split=',')), unlist(strsplit(y[which.row,'fhaart'],split=','))) old_drug_class <- ifelse(old_drug %in% Cs(RTV,TPV,LPV,DRV,SQV,NFV,ATV,IDV,APV,FPV),'PI', ifelse(old_drug %in% Cs(AZT,D4T,ABC,FTC,TDF,DDI),'NRTI', ifelse(old_drug %in% c('3TC','ETR','NVP','EFV'),'NNRTI', ifelse(old_drug %in% Cs(RAL),'Integrase inhibitor', ifelse(old_drug %in% Cs(MVC,ENF),'Fusion inhibitor',NA))))) new_drug_class <- ifelse(new_drug %in% Cs(RTV,TPV,LPV,DRV,SQV,NFV,ATV,IDV,APV,FPV),'PI', ifelse(new_drug %in% Cs(AZT,D4T,ABC,FTC,TDF,DDI),'NRTI', ifelse(new_drug %in% c('3TC','ETR','NVP','EFV'),'NNRTI', ifelse(new_drug %in% Cs(RAL),'Integrase inhibitor', ifelse(new_drug %in% Cs(MVC,ENF),'Fusion inhibitor',NA))))) switch_def_4 <- 0 switch_def_4_date <- NA if(!is.na(old_drug_class) & !is.na(new_drug_class) & old_drug_class != new_drug_class){ switch_def_4 <- 1 switch_def_4_date <- y[which.row,'single_drug_sub_date'] } } else { old_drug <- new_drug <- old_drug_class <- new_drug_class <- switch_def_4 <- switch_def_4_date <- NA } df <- data.frame('unique_id'=y[1,'unique_id'], 'old_drug'=old_drug, 'new_drug'=new_drug, 'old_drug_class'=old_drug_class, 'new_drug_class'=new_drug_class, 'switch_def_4'=switch_def_4, 'switch_def_4_date'=switch_def_4_date) return(df) }) sds.df <- do.call('rbind',sds) # Adding definition 4 information to art_after_fhaart art_after_fhaart$old_drug <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'old_drug'] art_after_fhaart$new_drug <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'new_drug'] art_after_fhaart$old_drug_class <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'old_drug_class'] art_after_fhaart$new_drug_class <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'new_drug_class'] art_after_fhaart$switch_def_4 <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'switch_def_4'] art_after_fhaart$switch_def_4_date <- sds.df[match(art_after_fhaart$unique_id, sds.df$unique_id, nomatch=NA),'switch_def_4_date'] art_after_fhaart$switch_def_4_date <- as.Date(art_after_fhaart$switch_def_4_date, origin='1970-01-01') art_after_fhaart$any_switch <- apply(art_after_fhaart[,c('switch_def_1','switch_def_2','switch_def_3','switch_def_4')],1,FUN=function(y){as.numeric(any(y == 1))}) art_after_fhaart$any_switch <- ifelse(is.na(art_after_fhaart$any_switch),0,art_after_fhaart$any_switch) art_after_fhaart$first_switch_date <- apply(subset(art_after_fhaart, select=c('switch_def_1_date','switch_def_2_date','switch_def_3_date','switch_def_4_date')),1,FUN=function(y){ tmp <- min(y,na.rm=TRUE) if(tmp %in% c(Inf,NA)) tmp <- NA return(tmp) }) art_after_fhaart$first_switch_date <- as.Date(art_after_fhaart$first_switch_date, format='%Y-%m-%d') art_after_fhaart$regimen_switched_to <- NA split(art_after_fhaart$regimen_switched_to, art_after_fhaart$unique_id) <- lapply(split(art_after_fhaart[,c('unique_id','art_id','first_switch_date','art_sd')],art_after_fhaart$unique_id),FUN=function(y){ print(y[1,'unique_id']) if(!all(is.na(y[,'first_switch_date']))){ regimen_switched_to <- rep(y[which(y[,'art_sd'] == y[,'first_switch_date']),'art_id'], times=nrow(y)) } else { regimen_switched_to <- rep(NA,nrow(y)) } }) art_after_fhaart$regimen_switched_to[which(art_after_fhaart$art_sd == art_after_fhaart$first_switch_date & !is.na(art_after_fhaart$first_switch_date))] <- art_after_fhaart$art_id[which(art_after_fhaart$art_sd == art_after_fhaart$first_switch_date & !is.na(art_after_fhaart$first_switch_date))] art_after_fhaart$next_start <- NA split(art_after_fhaart$next_start, art_after_fhaart$unique_id) <- lapply(split(subset(art_after_fhaart, select=c(art_sd)),art_after_fhaart$unique_id),FUN=function(y){ if(nrow(y) > 1){ c(y[2:nrow(y),'art_sd'],NA) } else { NA } }) art_after_fhaart$next_start <- as.Date(art_after_fhaart$next_start, origin='1970-01-01') art_after_fhaart$gap_in_trt <- as.numeric(art_after_fhaart$next_start - art_after_fhaart$art_ed) art_after_fhaart$gap_after_fhaart <- NA art_after_fhaart$gap_after_fhaart[which(art_after_fhaart$date_fhaart == art_after_fhaart$art_sd)] <- as.numeric(art_after_fhaart$next_start[which(art_after_fhaart$date_fhaart == art_after_fhaart$art_sd)] - art_after_fhaart$art_ed[which(art_after_fhaart$date_fhaart == art_after_fhaart$art_sd)]) tmp <- subset(art_after_fhaart, !is.na(gap_after_fhaart),select=c(unique_id, gap_after_fhaart)) art_after_fhaart$gap_after_fhaart <- tmp[match(art_after_fhaart$unique_id, tmp$unique_id, nomatch=NA),'gap_after_fhaart'] peds_data$first_switch_date <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'first_switch_date'] peds_data$regimen_switched_to <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'regimen_switched_to'] peds_data$fhaart <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'fhaart'] peds_data$switch_def_1_date <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'switch_def_1_date'] peds_data$switch_def_2_date <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'switch_def_2_date'] peds_data$switch_def_3_date <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'switch_def_3_date'] peds_data$switch_def_4_date <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'switch_def_4_date'] peds_data$switch_level <- ifelse(!is.na(peds_data$first_switch_date) & !is.na(peds_data$switch_def_1_date) & peds_data$first_switch_date == peds_data$switch_def_1_date,1, ifelse(!is.na(peds_data$first_switch_date) & !is.na(peds_data$switch_def_2_date) & peds_data$first_switch_date == peds_data$switch_def_2_date,2, ifelse(!is.na(peds_data$first_switch_date) & !is.na(peds_data$switch_def_3_date) & peds_data$first_switch_date == peds_data$switch_def_3_date,3, ifelse(!is.na(peds_data$first_switch_date) & !is.na(peds_data$switch_def_4_date) & peds_data$first_switch_date == peds_data$switch_def_4_date,4,0)))) peds_data$switch_level.factor <- factor(peds_data$switch_level, levels=c(0,1,2,3,4), labels=c('No switch','Switch >= 2 ARVs','Add >= 2 ARVs','Remove >= 1 ARVs leading to non-HAART class', 'Single out-of-class drug substitution')) peds_data$switch_outcome <- ifelse(peds_data$switch_level > 0,1,0) peds_data$switch_outcome.factor <- factor(peds_data$switch_outcome, levels=c(0,1), labels=c('No','Yes')) peds_data$switch_or_death_date <- apply(subset(peds_data, select=c(death_d,first_switch_date)),1,FUN=function(y){ min(y,na.rm=TRUE) }) peds_data$switch_or_death_date <- as.Date(peds_data$switch_or_death_date, format='%Y-%m-%d') peds_data$switch_or_death_outcome <- ifelse(!is.na(peds_data$switch_or_death_date),1,0) peds_data$time_from_fhaart_to_switch <- as.numeric(peds_data$first_switch_date - peds_data$date_fhaart) peds_data$time_from_fhaart_to_switch_or_death <- as.numeric(peds_data$switch_or_death_date - peds_data$date_fhaart) peds_data$followup_time_switch <- ifelse(is.na(peds_data$time_from_fhaart_to_switch), peds_data$followup_time, peds_data$time_from_fhaart_to_switch) peds_data$followup_time_switch_yrs <- peds_data$followup_time_switch/365.25 peds_data$log10_followup_time_switch_yrs <- log10(peds_data$followup_time_switch_yrs) peds_data$followup_time_combined <- ifelse(is.na(peds_data$time_from_fhaart_to_switch_or_death),peds_data$followup_time,peds_data$time_from_fhaart_to_switch_or_death) peds_data$followup_time_combined_yrs <- peds_data$followup_time_combined/365.25 peds_data$gap_after_fhaart <- art_after_fhaart[match(peds_data$unique_id, art_after_fhaart$unique_id, nomatch=NA),'gap_after_fhaart'] peds_data$crr_indicator <- ifelse(peds_data$switch_outcome == 0 & peds_data$death_y == 0,0, ifelse(peds_data$switch_outcome == 1,1, ifelse(peds_data$death_y == 1,2,NA))) # Removing those with a missing date of first HAART peds_data <- subset(peds_data, !is.na(date_fhaart)) label(peds_data$switch_level.factor) <- 'Switch definition' label(peds_data$time_from_fhaart_to_switch) <- 'Days from first cART to first switch' label(peds_data$time_from_fhaart_to_switch_or_death) <- 'Days from first cART to switch or death' label(peds_data$gap_after_fhaart) <- 'Days between end of first cART and start of next regimen' label(peds_data$followup_time_switch_yrs) <- 'Follow-up time switch outcome (yrs)' tbl1_switch <- summary(site_center.factor ~ switch_level.factor + time_from_fhaart_to_switch + time_from_fhaart_to_switch_or_death + gap_after_fhaart, data=peds_data, overall=TRUE, method='reverse') tmp_le_5 <- subset(peds_data, age_fhaart < 5) tmp_gt_5 <- subset(peds_data, age_fhaart >= 5) tmp_gt_5$alternate_site_center <- ifelse(tmp_gt_5$site_center %in% c('honduras.Escuela','honduras.IHSS'),'honduras',tmp_gt_5$site_center) tmp_gt_5$alternate_site_center.factor <- factor(tmp_gt_5$alternate_site_center) cum_inc_switch_le_5 <- cuminc(ftime=tmp_le_5$followup_time_switch_yrs, fstatus=tmp_le_5$crr_indicator, cencode='0',rho=0) cum_inc_switch_gt_5 <- cuminc(ftime=tmp_gt_5$followup_time_switch_yrs, fstatus=tmp_gt_5$crr_indicator, cencode='0',rho=0) plot(cum_inc_switch_le_5[[1]]$est[which(cum_inc_switch_le_5[[1]]$time <= 10)] ~ cum_inc_switch_le_5[[1]]$time[which(cum_inc_switch_le_5[[1]]$time <= 10)], type='l',axes=FALSE,xlim=c(0,10), xlab='Years from cART initiation',ylab='Cumulative incidence of switching regimen',lty=1,col='black', ylim=range(c(cum_inc_switch_le_5[[1]]$est[which(cum_inc_switch_le_5[[1]]$time <= 10)],cum_inc_switch_gt_5[[1]]$est[which(cum_inc_switch_gt_5[[1]]$time <= 10)])),lwd=1.5) lines(cum_inc_switch_gt_5[[1]]$est[which(cum_inc_switch_gt_5[[1]]$time <= 10)] ~ cum_inc_switch_gt_5[[1]]$time[which(cum_inc_switch_gt_5[[1]]$time <= 10)], lty=1,col='grey',lwd=1.5) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.95) # Overall cum_inc_switch_overall <- cuminc(ftime=peds_data$followup_time_switch_yrs, fstatus=peds_data$crr_indicator, cencode='0',rho=0) var_switch_overall <- cum_inc_switch_overall[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_overall[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_overall[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_overall[[1]]$time) == 5.0)) inc1_switch_overall <- cum_inc_switch_overall[[1]]$est[which.1] lci1_switch_overall <- inc1_switch_overall - qnorm(0.975)*sqrt(var_switch_overall[which.1]) uci1_switch_overall <- inc1_switch_overall + qnorm(0.975)*sqrt(var_switch_overall[which.1]) inc3_switch_overall <- cum_inc_switch_overall[[1]]$est[which.3] lci3_switch_overall <- inc3_switch_overall - qnorm(0.975)*sqrt(var_switch_overall[which.3]) uci3_switch_overall <- inc3_switch_overall + qnorm(0.975)*sqrt(var_switch_overall[which.3]) inc5_switch_overall <- cum_inc_switch_overall[[1]]$est[which.5] lci5_switch_overall <- inc5_switch_overall - qnorm(0.975)*sqrt(var_switch_overall[which.5]) uci5_switch_overall <- inc5_switch_overall + qnorm(0.975)*sqrt(var_switch_overall[which.5]) # For those < 5 var_switch_le5 <- cum_inc_switch_le_5[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_le_5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_le_5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_le_5[[1]]$time) == 5.0)) inc1_switch_le5 <- cum_inc_switch_le_5[[1]]$est[which.1] lci1_switch_le5 <- inc1_switch_le5 - qnorm(0.975)*sqrt(var_switch_le5[which.1]) uci1_switch_le5 <- inc1_switch_le5 + qnorm(0.975)*sqrt(var_switch_le5[which.1]) inc3_switch_le5 <- cum_inc_switch_le_5[[1]]$est[which.3] lci3_switch_le5 <- inc3_switch_le5 - qnorm(0.975)*sqrt(var_switch_le5[which.3]) uci3_switch_le5 <- inc3_switch_le5 + qnorm(0.975)*sqrt(var_switch_le5[which.3]) inc5_switch_le5 <- cum_inc_switch_le_5[[1]]$est[which.5] lci5_switch_le5 <- inc5_switch_le5 - qnorm(0.975)*sqrt(var_switch_le5[which.5]) uci5_switch_le5 <- inc5_switch_le5 + qnorm(0.975)*sqrt(var_switch_le5[which.5]) var5_switch_le5 <- var_switch_le5[which.5] # For those >= 5 var_switch_gt5 <- cum_inc_switch_gt_5[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_gt_5[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_gt_5[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_gt_5[[1]]$time) == 5.0)) inc1_switch_gt5 <- cum_inc_switch_gt_5[[1]]$est[which.1] lci1_switch_gt5 <- inc1_switch_gt5 - qnorm(0.975)*sqrt(var_switch_gt5[which.1]) uci1_switch_gt5 <- inc1_switch_gt5 + qnorm(0.975)*sqrt(var_switch_gt5[which.1]) inc3_switch_gt5 <- cum_inc_switch_gt_5[[1]]$est[which.3] lci3_switch_gt5 <- inc3_switch_gt5 - qnorm(0.975)*sqrt(var_switch_gt5[which.3]) uci3_switch_gt5 <- inc3_switch_gt5 + qnorm(0.975)*sqrt(var_switch_gt5[which.3]) inc5_switch_gt5 <- cum_inc_switch_gt_5[[1]]$est[which.5] lci5_switch_gt5 <- inc5_switch_gt5 - qnorm(0.975)*sqrt(var_switch_gt5[which.5]) uci5_switch_gt5 <- inc5_switch_gt5 + qnorm(0.975)*sqrt(var_switch_gt5[which.5]) var5_switch_gt5 <- var_switch_gt5[which.5] # Computing p-value comparing 5-yr cumulative incidence of death between those < 5 and those >= 5 z5_switch <- (inc5_switch_le5 - inc5_switch_gt5)/sqrt(var5_switch_le5 + var5_switch_gt5) p5_switch <- 2*(1-pnorm(z5_switch)) par(mfrow=c(1,2),mar=c(4,4,1,1),oma=c(1,1,1,1),pty='s') plot(cum_inc_le_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)] ~ cum_inc_le_5[[1]]$time[which(cum_inc_le_5[[1]]$time <= 10.0)], xlim=c(0,10), type='l', axes=FALSE, xlab='Years from cART initiation',ylab='Cumulative incidence',lty=1,col='black',cex.lab=0.75,cex.main=0.85, ylim=range(c(cum_inc_le_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)], cum_inc_gt_5[[1]]$est[which(cum_inc_le_5[[1]]$time <= 10.0)]),na.rm=TRUE),main='A. Mortality',lwd=1.5) lines(cum_inc_gt_5[[1]]$est[which(cum_inc_gt_5[[1]]$time <= 10.0)] ~ cum_inc_gt_5[[1]]$time[which(cum_inc_gt_5[[1]]$time <= 10.0)], lty=1,col='grey',lwd=1.5) axis(1, at=c(0,2,4,6,8,10), labels=c(0,2,4,6,8,10)) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.75) rge <- range(c(cum_inc_l2fu_le5[[1]]$est[which(cum_inc_l2fu_le5[[1]]$time <= 10.0)], cum_inc_l2fu_gt5[[1]]$est[which(cum_inc_l2fu_gt5[[1]]$time <= 10.0)])) plot(cum_inc_l2fu_le5[[1]]$est ~ cum_inc_l2fu_le5[[1]]$time, type='l',axes=FALSE,xlim=c(0,10),cex.lab=0.85,cex.main=0.85, xlab='Years from cART initiation',ylab='Cumulative incidence',ylim=rge,col='black',main='B. Lost to Follow-up',lwd=1.5) lines(cum_inc_l2fu_gt5[[1]]$est ~ cum_inc_l2fu_gt5[[1]]$time,col='grey',lwd=1.5) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.75) # Overall by site_center # Argentina cum_inc_switch_argentina <- cuminc(ftime=peds_data$followup_time_switch_yrs[which(peds_data$site_center.factor == 'argentina')], fstatus=peds_data$crr_indicator[which(peds_data$site_center.factor == 'argentina')], cencode='0',rho=0) var_switch_argentina <- cum_inc_switch_argentina[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_argentina[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_argentina[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_argentina[[1]]$time) == 5.0)) inc1_switch_argentina <- cum_inc_switch_argentina[[1]]$est[which.1] lci1_switch_argentina <- inc1_switch_argentina - qnorm(0.975)*sqrt(var_switch_argentina[which.1]) uci1_switch_argentina <- inc1_switch_argentina + qnorm(0.975)*sqrt(var_switch_argentina[which.1]) inc3_switch_argentina <- cum_inc_switch_argentina[[1]]$est[which.3] lci3_switch_argentina <- inc3_switch_argentina - qnorm(0.975)*sqrt(var_switch_argentina[which.3]) uci3_switch_argentina <- inc3_switch_argentina + qnorm(0.975)*sqrt(var_switch_argentina[which.3]) inc5_switch_argentina <- cum_inc_switch_argentina[[1]]$est[which.5] lci5_switch_argentina <- inc5_switch_argentina - qnorm(0.975)*sqrt(var_switch_argentina[which.5]) uci5_switch_argentina <- inc5_switch_argentina + qnorm(0.975)*sqrt(var_switch_argentina[which.5]) cum_inc_switch_brufmg <- cuminc(ftime=peds_data$followup_time_switch_yrs[which(peds_data$site_center.factor == 'brazil.UFMG')], fstatus=peds_data$crr_indicator[which(peds_data$site_center.factor == 'brazil.UFMG')], cencode='0',rho=0) var_switch_brufmg <- cum_inc_switch_brufmg[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_brufmg[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_brufmg[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_brufmg[[1]]$time) == 5.0)) inc1_switch_brufmg <- cum_inc_switch_brufmg[[1]]$est[which.1] lci1_switch_brufmg <- inc1_switch_brufmg - qnorm(0.975)*sqrt(var_switch_brufmg[which.1]) uci1_switch_brufmg <- inc1_switch_brufmg + qnorm(0.975)*sqrt(var_switch_brufmg[which.1]) inc3_switch_brufmg <- cum_inc_switch_brufmg[[1]]$est[which.3] lci3_switch_brufmg <- inc3_switch_brufmg - qnorm(0.975)*sqrt(var_switch_brufmg[which.3]) uci3_switch_brufmg <- inc3_switch_brufmg + qnorm(0.975)*sqrt(var_switch_brufmg[which.3]) inc5_switch_brufmg <- cum_inc_switch_brufmg[[1]]$est[which.5] lci5_switch_brufmg <- inc5_switch_brufmg - qnorm(0.975)*sqrt(var_switch_brufmg[which.5]) uci5_switch_brufmg <- inc5_switch_brufmg + qnorm(0.975)*sqrt(var_switch_brufmg[which.5]) cum_inc_switch_brunifesp <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'brazil.UNIFESP')], fstatus=peds_data$crr_indicator[which(peds_data$site_center.factor == 'brazil.UNIFESP')], cencode='0',rho=0) var_switch_brunifesp <- cum_inc_switch_brunifesp[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_brunifesp[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_brunifesp[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_brunifesp[[1]]$time) == 5.0)) inc1_switch_brunifesp <- cum_inc_switch_brunifesp[[1]]$est[which.1] lci1_switch_brunifesp <- inc1_switch_brunifesp - qnorm(0.975)*sqrt(var_switch_brunifesp[which.1]) uci1_switch_brunifesp <- inc1_switch_brunifesp + qnorm(0.975)*sqrt(var_switch_brunifesp[which.1]) inc3_switch_brunifesp <- cum_inc_switch_brunifesp[[1]]$est[which.3] lci3_switch_brunifesp <- inc3_switch_brunifesp - qnorm(0.975)*sqrt(var_switch_brunifesp[which.3]) uci3_switch_brunifesp <- inc3_switch_brunifesp + qnorm(0.975)*sqrt(var_switch_brunifesp[which.3]) inc5_switch_brunifesp <- cum_inc_switch_brunifesp[[1]]$est[which.5] lci5_switch_brunifesp <- inc5_switch_brunifesp - qnorm(0.975)*sqrt(var_switch_brunifesp[which.5]) uci5_switch_brunifesp <- inc5_switch_brunifesp + qnorm(0.975)*sqrt(var_switch_brunifesp[which.5]) cum_inc_switch_haiti <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'haiti')], fstatus=peds_data$crr_indicator[which(peds_data$site_center.factor == 'haiti')], cencode='0',rho=0) var_switch_haiti <- cum_inc_switch_haiti[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_haiti[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_switch_haiti[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_haiti[[1]]$time) == 5.0)) inc1_switch_haiti <- cum_inc_switch_haiti[[1]]$est[which.1] lci1_switch_haiti <- inc1_switch_haiti - qnorm(0.975)*sqrt(var_switch_haiti[which.1]) uci1_switch_haiti <- inc1_switch_haiti + qnorm(0.975)*sqrt(var_switch_haiti[which.1]) inc3_switch_haiti <- cum_inc_switch_haiti[[1]]$est[which.3] lci3_switch_haiti <- inc3_switch_haiti - qnorm(0.975)*sqrt(var_switch_haiti[which.3]) uci3_switch_haiti <- inc3_switch_haiti + qnorm(0.975)*sqrt(var_switch_haiti[which.3]) inc5_switch_haiti <- cum_inc_switch_haiti[[1]]$est[which.5] lci5_switch_haiti <- inc5_switch_haiti - qnorm(0.975)*sqrt(var_switch_haiti[which.5]) uci5_switch_haiti <- inc5_switch_haiti + qnorm(0.975)*sqrt(var_switch_haiti[which.5]) cum_inc_switch_honescuela <- cuminc(ftime=peds_data$followup_time_yrs[which(peds_data$site_center.factor == 'honduras.Escuela')], fstatus=peds_data$crr_indicator[which(peds_data$site_center.factor == 'honduras.Escuela')], cencode='0',rho=0) var_switch_honescuela <- cum_inc_switch_honescuela[[1]]$var which.1 <- max(which(ceiling(cum_inc_switch_honescuela[[1]]$time) == 1.0)) which.3 <- min(which(cum_inc_switch_honescuela[[1]]$time >= 3.0)) which.5 <- max(which(ceiling(cum_inc_switch_honescuela[[1]]$time) == 5.0)) inc1_switch_honescuela <- '0' inc3_switch_honescuela <- cum_inc_switch_honescuela[[1]]$est[which.3] lci3_switch_honescuela <- inc3_switch_honescuela - qnorm(0.975)*sqrt(var_switch_honescuela[which.3]) uci3_switch_honescuela <- inc3_switch_honescuela + qnorm(0.975)*sqrt(var_switch_honescuela[which.3]) inc5_switch_honescuela <- cum_inc_switch_honescuela[[1]]$est[which.5] lci5_switch_honescuela <- inc5_switch_honescuela - qnorm(0.975)*sqrt(var_switch_honescuela[which.5]) uci5_switch_honescuela <- inc5_switch_honescuela + qnorm(0.975)*sqrt(var_switch_honescuela[which.5]) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_switch_le_5 <- Surv(tmp_le_5$followup_time_switch_yrs,tmp_le_5$switch_outcome) Srv_switch_gt_5 <- Surv(tmp_gt_5$followup_time_switch_yrs,tmp_gt_5$switch_outcome) dd <- datadist(peds_data) options(datadist='dd') cph.switch.le5.fmla <- as.formula(paste('Srv_switch_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.switch.gt5.fmla <- as.formula(paste('Srv_switch_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE) mod_le_5_switch <- fit.mult.impute(cph.switch.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5 <- upData(tmp_gt_5, drop='switch_def_3_date') dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) mod_gt_5_switch <- fit.mult.impute(cph.switch.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) le_5_switch_df <- output_mgmt_no_hiv(mod_le_5_switch,le5_ind=1) gt_5_switch_df <- output_mgmt_no_hiv(mod_gt_5_switch,le5_ind=0) switch_df_combined <- cbind(le_5_switch_df, gt_5_switch_df[,-1]) # Table of cumulative incidences cum_inc_df <- data.frame('Outcome'=c('Death',rep('',times=2), 'Regimen change',rep('',times=2), 'Loss to follow-up',rep('',times=3)), 'Time'=c(rep(c('1 yr','3 yr','5 yr'),times=3),'7 yr'), '< 5'=c(paste(round(inc1_le5,2),' (',round(lci1_le5,2),', ',round(uci1_le5,2),')',sep=''), paste(round(inc3_le5,2),' (',round(lci3_le5,2),', ',round(uci3_le5,2),')',sep=''), paste(round(inc5_le5,2),' (',round(lci5_le5,2),', ',round(uci5_le5,2),')',sep=''), paste(round(inc1_switch_le5,2),' (',round(lci1_switch_le5,2),', ',round(uci1_switch_le5,2),')',sep=''), paste(round(inc3_switch_le5,2),' (',round(lci3_switch_le5,2),', ',round(uci3_switch_le5,2),')',sep=''), paste(round(inc5_switch_le5,2),' (',round(lci5_switch_le5,2),', ',round(uci5_switch_le5,2),')',sep=''), paste(round(inc1_l2fu_le5,2),' (',round(lci1_l2fu_le5,2),', ',round(uci1_l2fu_le5,2),')',sep=''), paste(round(inc3_l2fu_le5,2),' (',round(lci3_l2fu_le5,2),', ',round(uci3_l2fu_le5,2),')',sep=''), paste(round(inc5_l2fu_le5,2),' (',round(lci5_l2fu_le5,2),', ',round(uci5_l2fu_le5,2),')',sep=''), paste(round(inc7_l2fu_le5,2),' (',round(lci7_l2fu_le5,2),', ',round(uci7_l2fu_le5,2),')',sep='')), '>= 5'=c(paste(round(inc1_gt5,2),' (',round(lci1_gt5,2),', ',round(uci1_gt5,2),')',sep=''), paste(round(inc3_gt5,2),' (',round(lci3_gt5,2),', ',round(uci3_gt5,2),')',sep=''), paste(round(inc5_gt5,2),' (',round(lci5_gt5,2),', ',round(uci5_gt5,2),')',sep=''), paste(round(inc1_switch_gt5,2),' (',round(lci1_switch_gt5,2),', ',round(uci1_switch_gt5,2),')',sep=''), paste(round(inc3_switch_gt5,2),' (',round(lci3_switch_gt5,2),', ',round(uci3_switch_gt5,2),')',sep=''), paste(round(inc5_switch_gt5,2),' (',round(lci5_switch_gt5,2),', ',round(uci5_switch_gt5,2),')',sep=''), paste(round(inc1_l2fu_gt5,2),' (',round(lci1_l2fu_gt5,2),', ',round(uci1_l2fu_gt5,2),')',sep=''), paste(round(inc3_l2fu_gt5,2),' (',round(lci3_l2fu_gt5,2),', ',round(uci3_l2fu_gt5,2),')',sep=''), paste(round(inc5_l2fu_gt5,2),' (',round(lci5_l2fu_gt5,2),', ',round(uci5_l2fu_gt5,2),')',sep=''), paste(round(inc7_l2fu_gt5,2),' (',round(lci7_l2fu_gt5,2),', ',round(uci7_l2fu_gt5,2),')',sep='')), 'Overall'=c(paste(round(inc1_overall,2),' (',round(lci1_overall,2),', ',round(uci1_overall,2),')',sep=''), paste(round(inc3_overall,2),' (',round(lci3_overall,2),', ',round(uci3_overall,2),')',sep=''), paste(round(inc5_overall,2),' (',round(lci5_overall,2),', ',round(uci5_overall,2),')',sep=''), paste(round(inc1_switch_overall,2),' (',round(lci1_switch_overall,2),', ',round(uci1_switch_overall,2),')',sep=''), paste(round(inc3_switch_overall,2),' (',round(lci3_switch_overall,2),', ',round(uci3_switch_overall,2),')',sep=''), paste(round(inc5_switch_overall,2),' (',round(lci5_switch_overall,2),', ',round(uci5_switch_overall,2),')',sep=''), paste(round(inc1_l2fu_overall,2),' (',round(lci1_l2fu_overall,2),', ',round(uci1_l2fu_overall,2),')',sep=''), paste(round(inc3_l2fu_overall,2),' (',round(lci3_l2fu_overall,2),', ',round(uci3_l2fu_overall,2),')',sep=''), paste(round(inc5_l2fu_overall,2),' (',round(lci5_l2fu_overall,2),', ',round(uci5_l2fu_overall,2),')',sep=''), paste(round(inc7_l2fu_overall,2),' (',round(lci7_l2fu_overall,2),', ',round(uci7_l2fu_overall,2),')',sep='')),check.names=FALSE) # Table of cumulative incidences by site cum_inc_site_df <- data.frame('Site'=c('Argentina','','Brazil UFMG','','Brazil UNIFESP','','Haiti','','Honduras Escuela','','Honduras IHSS',''), 'Outcome'=rep(c('Death','Regimen change'),times=6), '1 yr'=c('0',paste(round(inc1_switch_argentina,2),' (',round(lci1_switch_argentina,2),', ',round(uci1_switch_argentina,2),')',sep=''), paste(round(inc1_brufmg,2),' (',round(lci1_brufmg,2),', ',round(uci1_brufmg,2),')',sep=''), paste(round(inc1_switch_brufmg,2),' (',round(lci1_switch_brufmg,2),', ',round(uci1_switch_brufmg,2),')',sep=''), paste(round(inc1_brunifesp,2),' (',round(lci1_brunifesp,2),', ',round(uci1_brunifesp,2),')',sep=''), paste(round(inc1_switch_brunifesp,2),' (',round(lci1_switch_brunifesp,2),', ',round(uci1_switch_brunifesp,2),')',sep=''), paste(round(inc1_haiti,2),' (',round(lci1_haiti,2),', ',round(uci1_haiti,2),')',sep=''), paste(round(inc1_switch_haiti,2),' (',round(lci1_switch_haiti,2),', ',round(uci1_switch_haiti,2),')',sep=''), paste(round(inc1_honescuela,2),' (',round(lci1_honescuela,2),', ',round(uci1_honescuela,2),')',sep=''), '0', '0','0'), '3 yr'=c('0',paste(round(inc3_switch_argentina,2),' (',round(lci3_switch_argentina,2),', ',round(uci3_switch_argentina,2),')',sep=''), paste(round(inc3_brufmg,2),' (',round(lci3_brufmg,2),', ',round(uci3_brufmg,2),')',sep=''), paste(round(inc3_switch_brufmg,2),' (',round(lci3_switch_brufmg,2),', ',round(uci3_switch_brufmg,2),')',sep=''), paste(round(inc3_brunifesp,2),' (',round(lci3_brunifesp,2),', ',round(uci3_brunifesp,2),')',sep=''), paste(round(inc3_switch_brunifesp,2),' (',round(lci3_switch_brunifesp,2),', ',round(uci3_switch_brunifesp,2),')',sep=''), paste(round(inc3_haiti,2),' (',round(lci3_haiti,2),', ',round(uci3_haiti,2),')',sep=''), paste(round(inc3_switch_haiti,2),' (',round(lci3_switch_haiti,2),', ',round(uci3_switch_haiti,2),')',sep=''), paste(round(inc3_honescuela,2),' (',round(lci3_honescuela,2),', ',round(uci3_honescuela,2),')',sep=''),'0', '0','0'), '5 yr'=c('0',paste(round(inc5_switch_argentina,2),' (',round(lci5_switch_argentina,2),', ',round(uci5_switch_argentina,2),')',sep=''), paste(round(inc5_brufmg,2),' (',round(lci5_brufmg,2),', ',round(uci5_brufmg,2),')',sep=''), paste(round(inc5_switch_brufmg,2),' (',round(lci5_switch_brufmg,2),', ',round(uci5_switch_brufmg,2),')',sep=''), paste(round(inc5_brunifesp,2),' (',round(lci5_brunifesp,2),', ',round(uci5_brunifesp,2),')',sep=''), paste(round(inc5_switch_brunifesp,2),' (',round(lci5_switch_brunifesp,2),', ',round(uci5_switch_brunifesp,2),')',sep=''), paste(round(inc5_haiti,2),' (',round(lci5_haiti,2),', ',round(uci5_haiti,2),')',sep=''), paste(round(inc5_switch_haiti,2),' (',round(lci5_switch_haiti,2),', ',round(uci5_switch_haiti,2),')',sep=''), paste(round(inc5_honescuela,2),' (',round(lci5_honescuela,2),', ',round(uci5_honescuela,2),')',sep=''), paste(round(inc5_switch_honescuela,2),' (',round(lci5_switch_honescuela,2),', ',round(uci5_switch_honescuela,2),')',sep=''),'0','0'),check.names=FALSE) sub_peds_data <- subset(peds_data, site_center.factor %in% c('brazil.UFMG','brazil.UNIFESP','haiti')) sub_peds_data$site_center.factor <- factor(sub_peds_data$site_center.factor) sub_peds_data$alternate_site_center.factor <- factor(sub_peds_data$alternate_site_center.factor) tmp_le_5 <- subset(sub_peds_data, age_fhaart < 5) tmp_gt_5 <- subset(sub_peds_data, age_fhaart >= 5) cum_inc_le_5_sub <- cuminc(ftime=tmp_le_5$followup_time_yrs, fstatus=tmp_le_5$death_y, cencode='0',rho=0) cum_inc_gt_5_sub <- cuminc(ftime=tmp_gt_5$followup_time_yr, fstatus=tmp_gt_5$death_y, cencode='0',rho=0) plot(cum_inc_le_5_sub[[1]]$est[which(cum_inc_le_5_sub[[1]]$time <= 10)] ~ cum_inc_le_5_sub[[1]]$time[which(cum_inc_le_5_sub[[1]]$time <= 10)], type='l',axes=FALSE,xlab='Years from cART initiation', xlim=c(0,10),ylab='Cumulative incidence',col='black', ylim=range(c(cum_inc_le_5_sub[[1]]$est[which(cum_inc_le_5_sub[[1]]$time <= 10)],cum_inc_gt_5_sub[[1]]$est[which(cum_inc_gt_5_sub[[1]]$time <= 10)])),lwd=1.5) lines(cum_inc_gt_5_sub[[1]]$est[which(cum_inc_gt_5_sub[[1]]$time <= 10)] ~ cum_inc_gt_5_sub[[1]]$time[which(cum_inc_gt_5_sub[[1]]$time <= 10)], col='grey',lwd=1.5) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.75) # Need overall cuminc cum_inc_overall_sub <- cuminc(ftime=sub_peds_data$followup_time_yrs, fstatus=sub_peds_data$death_y, cencode='0',rho=0) which.1 <- max(which(ceiling(cum_inc_overall_sub[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_overall_sub[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_overall_sub[[1]]$time) == 5.0)) var_overall_sub <- cum_inc_overall_sub[[1]]$var inc1_overall_sub <- cum_inc_overall_sub[[1]]$est[which.1] lci1_overall_sub <- inc1_overall_sub - qnorm(0.975)*sqrt(var_overall_sub[which.1]) uci1_overall_sub <- inc1_overall_sub + qnorm(0.975)*sqrt(var_overall_sub[which.1]) inc3_overall_sub <- cum_inc_overall_sub[[1]]$est[which.3] lci3_overall_sub <- inc3_overall_sub - qnorm(0.975)*sqrt(var_overall_sub[which.3]) uci3_overall_sub <- inc3_overall_sub + qnorm(0.975)*sqrt(var_overall_sub[which.3]) inc5_overall_sub <- cum_inc_overall_sub[[1]]$est[which.5] lci5_overall_sub <- inc5_overall_sub - qnorm(0.975)*sqrt(var_overall_sub[which.5]) uci5_overall_sub <- inc5_overall_sub + qnorm(0.975)*sqrt(var_overall_sub[which.5]) # For those < 5 var_le5_sub <- cum_inc_le_5_sub[[1]]$var which.1 <- max(which(ceiling(cum_inc_le_5_sub[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_le_5_sub[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_le_5_sub[[1]]$time) == 5.0)) inc1_le5_sub <- cum_inc_le_5_sub[[1]]$est[which.1] lci1_le5_sub <- inc1_le5_sub - qnorm(0.975)*sqrt(var_le5_sub[which.1]) uci1_le5_sub <- inc1_le5_sub + qnorm(0.975)*sqrt(var_le5_sub[which.1]) inc3_le5_sub <- cum_inc_le_5_sub[[1]]$est[which.3] lci3_le5_sub <- inc3_le5_sub - qnorm(0.975)*sqrt(var_le5_sub[which.3]) uci3_le5_sub <- inc3_le5_sub + qnorm(0.975)*sqrt(var_le5_sub[which.3]) inc5_le5_sub <- cum_inc_le_5_sub[[1]]$est[which.5] lci5_le5_sub <- inc5_le5_sub - qnorm(0.975)*sqrt(var_le5_sub[which.5]) uci5_le5_sub <- inc5_le5_sub + qnorm(0.975)*sqrt(var_le5_sub[which.5]) var5_le5_sub <- var_le5_sub[which.5] # For those >= 5 var_gt5_sub <- cum_inc_gt_5_sub[[1]]$var which.1 <- max(which(ceiling(cum_inc_gt_5_sub[[1]]$time) == 1.0)) which.3 <- max(which(ceiling(cum_inc_gt_5_sub[[1]]$time) == 3.0)) which.5 <- max(which(ceiling(cum_inc_gt_5_sub[[1]]$time) == 5.0)) inc1_gt5_sub <- cum_inc_gt_5_sub[[1]]$est[which.1] lci1_gt5_sub <- inc1_gt5_sub - qnorm(0.975)*sqrt(var_gt5_sub[which.1]) uci1_gt5_sub <- inc1_gt5_sub + qnorm(0.975)*sqrt(var_gt5_sub[which.1]) inc3_gt5_sub <- cum_inc_gt_5_sub[[1]]$est[which.3] lci3_gt5_sub <- inc3_gt5_sub - qnorm(0.975)*sqrt(var_gt5_sub[which.3]) uci3_gt5_sub <- inc3_gt5_sub + qnorm(0.975)*sqrt(var_gt5_sub[which.3]) inc5_gt5_sub <- cum_inc_gt_5_sub[[1]]$est[which.5] lci5_gt5_sub <- inc5_gt5_sub - qnorm(0.975)*sqrt(var_gt5_sub[which.5]) uci5_gt5_sub <- inc5_gt5_sub + qnorm(0.975)*sqrt(var_gt5_sub[which.5]) var5_gt5_sub <- var_gt5_sub[which.5] # Table of cumulative incidences cum_inc_df_sub <- data.frame('Outcome'=c('Death',rep('',times=2)), 'Time'=c(rep(c('1 yr','3 yr','5 yr'),times=1)), '< 5'=c(paste(round(inc1_le5_sub,2),' (',round(lci1_le5_sub,2),', ',round(uci1_le5_sub,2),')',sep=''), paste(round(inc3_le5_sub,2),' (',round(lci3_le5_sub,2),', ',round(uci3_le5_sub,2),')',sep=''), paste(round(inc5_le5_sub,2),' (',round(lci5_le5_sub,2),', ',round(uci5_le5_sub,2),')',sep='')), '>= 5'=c(paste(round(inc1_gt5_sub,2),' (',round(lci1_gt5_sub,2),', ',round(uci1_gt5_sub,2),')',sep=''), paste(round(inc3_gt5_sub,2),' (',round(lci3_gt5_sub,2),', ',round(uci3_gt5_sub,2),')',sep=''), paste(round(inc5_gt5_sub,2),' (',round(lci5_gt5_sub,2),', ',round(uci5_gt5_sub,2),')',sep='')), 'Overall'=c(paste(round(inc1_overall_sub,2),' (',round(lci1_overall_sub,2),', ',round(uci1_overall_sub,2),')',sep=''), paste(round(inc3_overall_sub,2),' (',round(lci3_overall_sub,2),', ',round(uci3_overall_sub,2),')',sep=''), paste(round(inc5_overall_sub,2),' (',round(lci5_overall_sub,2),', ',round(uci5_overall_sub,2),')',sep='')), check.names=FALSE) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_le_5 <- Surv(tmp_le_5$followup_time_yrs,tmp_le_5$death_y) Srv_gt_5 <- Surv(tmp_gt_5$followup_time_yrs,tmp_gt_5$death_y) var_list_le_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4_per),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, pre_art_aids.clinical, strat(site_center.factor)) var_list_gt_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, pre_art_aids.clinical, strat(site_center.factor)) cph.le5.fmla <- as.formula(paste('Srv_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.gt5.fmla <- as.formula(paste('Srv_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_impute_cd4 <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4_per) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5 <- fit.mult.impute(cph.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_impute_cd4, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_impute_cd4 <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5 <- fit.mult.impute(cph.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_impute_cd4, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_df <- output_mgmt_cd4_per(sub_mod_le_5,le5_ind=1) sub_gt_5_df <- output_mgmt_no_hiv(sub_mod_gt_5,le5_ind=0) cum_inc_switch_le_5_sub <- cuminc(ftime=tmp_le_5$followup_time_switch_yrs, fstatus=tmp_le_5$crr_indicator, cencode='0',rho=0) cum_inc_switch_gt_5_sub <- cuminc(ftime=tmp_gt_5$followup_time_switch_yrs, fstatus=tmp_gt_5$crr_indicator, cencode='0',rho=0) plot(cum_inc_switch_le_5_sub[[1]]$est[which(cum_inc_switch_le_5_sub[[1]]$time <= 10)] ~ cum_inc_switch_le_5_sub[[1]]$time[which(cum_inc_switch_le_5_sub[[1]]$time <= 10)], type='l',axes=FALSE, xlim=c(0,10), xlab='Years from cART initiation',ylab='Cumulative incidence of switching regimen',col='black', ylim=range(c(cum_inc_switch_le_5_sub[[1]]$est[which(cum_inc_switch_le_5_sub[[1]]$time <= 10)],cum_inc_switch_gt_5_sub[[1]]$est[which(cum_inc_switch_gt_5_sub[[1]]$time <= 10)])),lwd=1.5) lines(cum_inc_switch_gt_5_sub[[1]]$est[which(cum_inc_switch_gt_5_sub[[1]]$time <= 10)] ~ cum_inc_switch_gt_5_sub[[1]]$time[which(cum_inc_switch_gt_5_sub[[1]]$time <= 10)], col='grey',lwd=1.5) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),fill=c('black','grey'),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_switch_le_5 <- Surv(tmp_le_5$followup_time_switch_yrs,tmp_le_5$switch_outcome) Srv_switch_gt_5 <- Surv(tmp_gt_5$followup_time_switch_yrs,tmp_gt_5$switch_outcome) cph.switch.le5.fmla <- as.formula(paste('Srv_switch_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.switch.gt5.fmla <- as.formula(paste('Srv_switch_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4_per) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5_switch <- fit.mult.impute(cph.switch.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5$site.factor <- factor(tmp_gt_5$site) tmp_gt_5 <- upData(tmp_gt_5, drop='switch_def_3_date') dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5_switch <- fit.mult.impute(cph.switch.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_switch_df <- output_mgmt_cd4_per(sub_mod_le_5_switch,le5_ind=1) sub_gt_5_switch_df <- output_mgmt_no_hiv(sub_mod_gt_5_switch,le5_ind=0) sub_peds_data2 <- subset(peds_data, site_center.factor %in% c('argentina','brazil.UFMG','brazil.UNIFESP')) sub_peds_data2$site_center.factor <- factor(sub_peds_data2$site_center.factor) sub_peds_data2$alternate_site_center.factor <- factor(sub_peds_data2$alternate_site_center.factor) # Because of issues with the first cART class by death, combining 'Other' with 'PI-only' cART sub_peds_data2$fhaart_class <- ifelse(sub_peds_data2$fhaart_class == 3, 2, sub_peds_data2$fhaart_class) sub_peds_data2$fhaart_class.factor <- factor(sub_peds_data2$fhaart_class, levels=c(1,2), labels=c('NNRTI-only','PI-only or Other')) tmp_le_5 <- subset(sub_peds_data2, age_fhaart < 5) tmp_gt_5 <- subset(sub_peds_data2, age_fhaart >= 5) cum_inc_le_5_subvl <- cuminc(ftime=tmp_le_5$followup_time_yrs, fstatus=tmp_le_5$death_y, cencode='0',rho=0) cum_inc_gt_5_subvl <- cuminc(ftime=tmp_gt_5$followup_time_yr, fstatus=tmp_gt_5$death_y, cencode='0',rho=0) plot(cum_inc_le_5_subvl[[1]]$est[which(cum_inc_le_5_subvl[[1]]$time <= 10)] ~ cum_inc_le_5_subvl[[1]]$time[which(cum_inc_le_5_subvl[[1]]$time <= 10)], type='l',axes=FALSE, xlim=c(0,10),xlab='Years from cART initiation',ylab='Cumulative incidence',col='black', lty=2, ylim=range(c(cum_inc_le_5_subvl[[1]]$est[which(cum_inc_le_5_subvl[[1]]$time <= 10)] ,cum_inc_gt_5_subvl[[1]]$est[which(cum_inc_gt_5_subvl[[1]]$time <= 10)] )),lwd=1.5) lines(cum_inc_gt_5_subvl[[1]]$est[which(cum_inc_gt_5_subvl[[1]]$time <= 10)] ~ cum_inc_gt_5_subvl[[1]]$time[which(cum_inc_gt_5_subvl[[1]]$time <= 10)], col='grey',lwd=1.5,lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),lty=c(2,3),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_le_5 <- Surv(tmp_le_5$followup_time_yrs,tmp_le_5$death_y) Srv_gt_5 <- Surv(tmp_gt_5$followup_time_yrs,tmp_gt_5$death_y) var_list_le_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, rcs(log10_baseline_rna,3), pre_art_aids.clinical, strat(site_center.factor)) var_list_gt_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, rcs(log10_baseline_rna,3), pre_art_aids.clinical, strat(site_center.factor)) cph.le5.fmla <- as.formula(paste('Srv_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.gt5.fmla <- as.formula(paste('Srv_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_impute_vl <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5_vl <- fit.mult.impute(cph.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_impute_vl, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# dd <- datadist(tmp_gt_5) options(datadist='dd') set.seed(10) gt_5_impute_vl <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5_vl <- fit.mult.impute(cph.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_impute_vl, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_vl_df <- output_mgmt(sub_mod_le_5_vl,le5_ind=1,vl_restrict=1) sub_gt_5_vl_df <- output_mgmt(sub_mod_gt_5_vl,le5_ind=0,vl_restrict=1) cum_inc_switch_le_5_subvl <- cuminc(ftime=tmp_le_5$followup_time_switch_yrs, fstatus=tmp_le_5$crr_indicator, cencode='0',rho=0) cum_inc_switch_gt_5_subvl <- cuminc(ftime=tmp_gt_5$followup_time_switch_yrs, fstatus=tmp_gt_5$crr_indicator, cencode='0',rho=0) plot(cum_inc_switch_le_5_subvl[[1]]$est[which(cum_inc_switch_le_5_subvl[[1]]$time <= 10)] ~ cum_inc_switch_le_5_subvl[[1]]$time[which(cum_inc_switch_le_5_subvl[[1]]$time <= 10)], type='l',axes=FALSE, xlim=c(0,10),xlab='Years from cART initiation',ylab='Cumulative incidence of switching regimen',col='black',lty=2, ylim=range(c(cum_inc_switch_le_5_subvl[[1]]$est[which(cum_inc_switch_le_5_subvl[[1]]$time <= 10)], cum_inc_switch_gt_5_subvl[[1]]$est[which(cum_inc_switch_gt_5_subvl[[1]]$time <= 10)])),lwd=1.5) lines(cum_inc_switch_gt_5_subvl[[1]]$est[which(cum_inc_switch_gt_5_subvl[[1]]$time <= 10)] ~ cum_inc_switch_gt_5_subvl[[1]]$time[which(cum_inc_switch_gt_5_subvl[[1]]$time <= 10)], col='grey',lwd=1.5,lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),lty=c(2,3),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_switch_le_5 <- Surv(tmp_le_5$followup_time_switch_yrs,tmp_le_5$switch_outcome) Srv_switch_gt_5 <- Surv(tmp_gt_5$followup_time_switch_yrs,tmp_gt_5$switch_outcome) cph.switch.le5.fmla <- as.formula(paste('Srv_switch_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.switch.gt5.fmla <- as.formula(paste('Srv_switch_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5_switch_vl <- fit.mult.impute(cph.switch.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5$site.factor <- factor(tmp_gt_5$site) tmp_gt_5 <- upData(tmp_gt_5, drop='switch_def_3_date') dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5_switch_vl <- fit.mult.impute(cph.switch.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_switch_vl_df <- output_mgmt(sub_mod_le_5_switch_vl,le5_ind=1,vl_restrict=1) sub_gt_5_switch_vl_df <- output_mgmt(sub_mod_gt_5_switch_vl,le5_ind=0,vl_restrict=1) sub_peds_data3 <- subset(peds_data, site_center.factor %in% c('brazil.UFMG','brazil.UNIFESP')) sub_peds_data3$site_center.factor <- factor(sub_peds_data3$site_center.factor) sub_peds_data3$alternate_site_center.factor <- factor(sub_peds_data3$alternate_site_center.factor) # Because of issues with the first cART class by death, combining 'Other' with 'PI-only' cART sub_peds_data3$fhaart_class <- ifelse(sub_peds_data3$fhaart_class == 3, 2, sub_peds_data3$fhaart_class) sub_peds_data3$fhaart_class.factor <- factor(sub_peds_data3$fhaart_class, levels=c(1,2), labels=c('NNRTI-only','PI-only or Other')) tmp_le_5 <- subset(sub_peds_data3, age_fhaart < 5) tmp_gt_5 <- subset(sub_peds_data3, age_fhaart >= 5) cum_inc_le_5_subvlcd4 <- cuminc(ftime=tmp_le_5$followup_time_yrs, fstatus=tmp_le_5$death_y, cencode='0',rho=0) cum_inc_gt_5_subvlcd4 <- cuminc(ftime=tmp_gt_5$followup_time_yr, fstatus=tmp_gt_5$death_y, cencode='0',rho=0) plot(cum_inc_le_5_subvlcd4[[1]]$est[which(cum_inc_le_5_subvlcd4[[1]]$time <= 10)] ~ cum_inc_le_5_subvlcd4[[1]]$time[which(cum_inc_le_5_subvlcd4[[1]]$time <= 10)], type='l',axes=FALSE, xlim=c(0,10),xlab='Years from cART initiation',ylab='Cumulative incidence',col='black',lty=2, ylim=range(c(cum_inc_le_5_subvlcd4[[1]]$est[which(cum_inc_le_5_subvlcd4[[1]]$time <= 10)],cum_inc_gt_5_subvlcd4[[1]]$est[which(cum_inc_gt_5_subvlcd4[[1]]$time <= 10)])),lwd=1.5) lines(cum_inc_gt_5_subvlcd4[[1]]$est[which(cum_inc_gt_5_subvlcd4[[1]]$time <= 10)] ~ cum_inc_gt_5_subvlcd4[[1]]$time[which(cum_inc_gt_5_subvlcd4[[1]]$time <= 10)], col='black',lwd=1.5,lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),lty=c(2,3),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_le_5 <- Surv(tmp_le_5$followup_time_yrs,tmp_le_5$death_y) Srv_gt_5 <- Surv(tmp_gt_5$followup_time_yrs,tmp_gt_5$death_y) var_list_le_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4_per),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, rcs(log10_baseline_rna,3), pre_art_aids.clinical, strat(site_center.factor)) var_list_gt_5 <- Cs(rcs(age_fhaart,3), male.factor, rcs(sqrt(baseline_cd4),3), rcs(year_fhaart,3), fhaart_class_dichot.factor, rcs(log10_baseline_rna,3), pre_art_aids.clinical, strat(site_center.factor)) cph.le5.fmla <- as.formula(paste('Srv_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.gt5.fmla <- as.formula(paste('Srv_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_impute_cd4_vl <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4_per) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5_vl_cd4 <- fit.mult.impute(cph.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_impute_cd4_vl, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5 <- upData(tmp_gt_5, drop='switch_def_3_date') dd <- datadist(tmp_gt_5) options(datadist='dd') set.seed(10) gt_5_impute_cd4_vl <- aregImpute(~ log10_followup_time_yrs*death_y + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5_vl_cd4 <- fit.mult.impute(cph.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_impute_cd4_vl, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_vl_cd4_df <- output_mgmt_vl_cd4_per(sub_mod_le_5_vl_cd4,le5_ind=1,vl_restrict=1) sub_gt_5_vl_cd4_df <- output_mgmt(sub_mod_gt_5_vl_cd4,le5_ind=0,vl_restrict=1) cum_inc_switch_le_5_subvlcd4 <- cuminc(ftime=tmp_le_5$followup_time_switch_yrs, fstatus=tmp_le_5$crr_indicator, cencode='0',rho=0) cum_inc_switch_gt_5_subvlcd4 <- cuminc(ftime=tmp_gt_5$followup_time_switch_yrs, fstatus=tmp_gt_5$crr_indicator, cencode='0',rho=0) plot(cum_inc_switch_le_5_subvlcd4[[1]]$est[which(cum_inc_switch_le_5_subvlcd4[[1]]$time <= 10)] ~ cum_inc_switch_le_5_subvlcd4[[1]]$time[which(cum_inc_switch_le_5_subvlcd4[[1]]$time <= 10)], type='l', xlim=c(0,10),axes=FALSE,xlab='Years from cART initiation',ylab='Cumulative incidence of switching regimen',col='black', ylim=range(c(cum_inc_switch_le_5_subvlcd4[[1]]$est[which(cum_inc_switch_le_5_subvlcd4[[1]]$time <= 10)],cum_inc_switch_gt_5_subvlcd4[[1]]$est[which(cum_inc_switch_gt_5_subvlcd4[[1]]$time <= 10)])),lwd=1.5, lty=2) lines(cum_inc_switch_gt_5_subvlcd4[[1]]$est[which(cum_inc_switch_gt_5_subvlcd4[[1]]$time <= 10)] ~ cum_inc_switch_gt_5_subvlcd4[[1]]$time[which(cum_inc_switch_gt_5_subvlcd4[[1]]$time <= 10)], col='red',lwd=1.5, lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=expression(phantom() < '5 years old', phantom() >= '5 years old'),lty=c(2,3),bty='n',cex=0.75) #------------------------------------# # Time-to-event analysis # #------------------------------------# Srv_switch_le_5 <- Surv(tmp_le_5$followup_time_switch_yrs,tmp_le_5$switch_outcome) Srv_switch_gt_5 <- Surv(tmp_gt_5$followup_time_switch_yrs,tmp_gt_5$switch_outcome) cph.switch.le5.fmla <- as.formula(paste('Srv_switch_le_5 ~ ',paste(var_list_le_5, collapse=' + '),sep='')) cph.switch.gt5.fmla <- as.formula(paste('Srv_switch_gt_5 ~ ',paste(var_list_gt_5, collapse=' + '),sep='')) #---------------------------# # < 5 years old # #---------------------------# dd <- datadist(tmp_le_5) options(datadist='dd') set.seed(2) le_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4_per) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + site_center.factor, data=tmp_le_5, x=TRUE,n.impute=5) sub_mod_le_5_switch_vl_cd4 <- fit.mult.impute(cph.switch.le5.fmla, data=tmp_le_5, fitter=cph, xtrans=le_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) #----------------------------# # >= 5 years old # #----------------------------# tmp_gt_5$site.factor <- factor(tmp_gt_5$site) tmp_gt_5 <- upData(tmp_gt_5, drop='switch_def_3_date') dd <- datadist(tmp_gt_5) options(datadist='dd') gt_5_switch_impute <- aregImpute(~ log10_followup_time_switch_yrs*switch_outcome + age_fhaart + male.factor + sqrt(baseline_cd4) + year_fhaart + fhaart_class_dichot.factor + log10_baseline_rna + pre_art_aids.clinical + alternate_site_center.factor, data=tmp_gt_5, x=TRUE,n.impute=5) sub_mod_gt_5_switch_vl_cd4 <- fit.mult.impute(cph.switch.gt5.fmla, data=tmp_gt_5, fitter=cph, xtrans=gt_5_switch_impute, x=TRUE, y=TRUE, surv=TRUE,n.impute=5) sub_le_5_switch_vl_cd4_df <- output_mgmt_vl_cd4_per(sub_mod_le_5_switch_vl_cd4,le5_ind=1,vl_restrict=1) sub_gt_5_switch_vl_cd4_df <- output_mgmt(sub_mod_gt_5_switch_vl_cd4,le5_ind=0,vl_restrict=1) # For abstract/manuscript cum_inc_switch_overall[[1]]$time[which(round(cum_inc_switch_overall[[1]]$est,3) == 0.50)] cum_inc_switch_le_5[[1]]$time[which(round(cum_inc_switch_le_5[[1]]$est,2) == 0.50)] cum_inc_switch_gt_5[[1]]$time[which(round(cum_inc_switch_gt_5[[1]]$est,2) == 0.50)] stuff<-with(peds_data, boxplot(age_fhaart~year_fhaart)) stuff.a<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="argentina")) stuff.b.ufmg<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="brazil.UFMG")) stuff.b.unifesp<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="brazil.UNIFESP")) stuff.h<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="haiti")) stuff.ho.escuela<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="honduras.Escuela")) stuff.ho.IHSS<-with(peds_data, boxplot(age_fhaart~year_fhaart,subset=site_center=="honduras.IHSS")) plot(c(1997,2013.2),c(0,18),type="n",xlab="Year of cART initiation",ylab="Median age",axes=FALSE) points(stuff.a$names,stuff.a$stats[3,],cex=sqrt(stuff.a$n)/(1.5*pi), pch=19, col=2) lines(stuff.a$names,stuff.a$stats[3,],col=2) points(stuff.b.ufmg$names,stuff.b.ufmg$stats[3,],cex=sqrt(stuff.b.ufmg$n)/(1.5*pi), pch=19,col=3) lines(stuff.b.ufmg$names,stuff.b.ufmg$stats[3,],col=3) points(stuff.b.unifesp$names,stuff.b.unifesp$stats[3,],cex=sqrt(stuff.b.unifesp$n)/(1.5*pi), pch=19,col=4) lines(stuff.b.unifesp$names,stuff.b.unifesp$stats[3,],col=4) points(stuff.h$names,stuff.h$stats[3,],cex=sqrt(stuff.h$n)/(1.5*pi), pch=19,col=5) lines(stuff.h$names,stuff.h$stats[3,],col=5) points(stuff.ho.escuela$names,stuff.ho.escuela$stats[3,],cex=sqrt(stuff.ho.escuela$n)/(1.5*pi), pch=19,col=6) lines(stuff.ho.escuela$names,stuff.ho.escuela$stats[3,],col=6) points(stuff.ho.IHSS$names,stuff.ho.IHSS$stats[3,],cex=sqrt(stuff.ho.IHSS$n)/(1.5*pi), pch=19,col=7) lines(stuff.ho.IHSS$names,stuff.ho.IHSS$stats[3,],col=7) lines(stuff$names,stuff$stats[3,],lwd=3) axis(1,at=seq(1997,2013,by=2)) axis(2) box() legend("topleft",c("HF-Argentina","UFMG-Brazil","UNIFESP-Brazil","GHESKIO-Haiti","HE-Honduras","IHSS-Honduras","Combined"), col=c(2:7,1),lty=rep(1,7),pch=rep(19,7))