#-----------------------------# # Some default settings # #-----------------------------# rm(list=ls()) options(width=180) #-----------------------------# # Loading libraries # #-----------------------------# library(Hmisc) library(cmprsk) library(rms) library(mi) library(plyr) #-----------------------------# # Functions # #-----------------------------# source('first_failure_virologic_failure_function.R') source('time_to_evt_function.R') source('/home/jenkinca/InfectiousDiseases/CathyMcGowan/3rdLineHaart/time_to_evt_nnrti_function.R') source('/home/jenkinca/InfectiousDiseases/CathyMcGowan/3rdLineHaart/rct_descriptives.R') source('/home/jenkinca/InfectiousDiseases/CathyMcGowan/3rdLineHaart/find_max_date.R') source('/home/jenkinca/InfectiousDiseases/CathyMcGowan/3rdLineHaart/aids_at_fhaart_function.R') source('CCASAnet_ADE_definition_R_code.R') source('unadj_time_2_evt.R') source('create_unadj_fmla.R') summary.fun <- function(include_haiti,mod,cd4.comp,age.comp,year.comp,vl.comp,reduced){ if(include_haiti | reduced){ summ. <- summary(mod, male.factor='Female', baseline_cd4=cd4.comp,aids_at_fhaart.factor='Not AIDS', fhaart_regimen_class.factor='NNRTI',nrti_type.factor = 'AZT', age_fhaart=age.comp, year_fhaart=year.comp) } else { summ. <- summary(mod, male.factor='Female', baseline_cd4=cd4.comp,aids_at_fhaart.factor='Not AIDS', fhaart_regimen_class.factor='NNRTI',nrti_type.factor = 'AZT', age_fhaart=age.comp, year_fhaart=year.comp, transmission_risk_larger_cat.factor='Heterosexual', baseline_log_rna=vl.comp) } return(summ.) } #-----------------------------# # Loading data # #-----------------------------# art <- read.csv('Data/art.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) art <- subset(art, !is.na(patient)) basic <- read.csv('Data/basic.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'','.')) basic$enrol_d <- basic$baseline_d basic <- upData(basic, drop='baseline_d') follow <- read.csv('Data/follow.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) cd4 <- read.csv('Data/lab_cd4.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) rna <- read.csv('Data/lab_rna.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) visit <- read.csv('Data/visit.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) ltfu <- read.csv('~/InfectiousDiseases/CathyMcGowan/3rdLineHaart/Data/third-line-ltfu.csv',header=TRUE,sep=',',stringsAsFactors=FALSE,na.strings=c(NA,'')) #------------------------------# # HAART regimen codes # #------------------------------# regimen_codes <- data.frame('Drug_type'=c(rep('Entry inhibitors',2), rep('NNRTI',4), rep('NRTI',7), 'Integrase inhibitors', rep('PI',9)), 'Drug_name'=c('Enfuvirtide','Maraviroc','Delavirdine', 'Efavirenz','Etravirine','Nevirapine', 'Abacavir','Didanosine','Emtricitabine', 'Lamivudine','Stavudine','Tenofovir', 'Zidovudine','Raltegravir','Atazanavir', 'Darunavir','Fosamprenavir','Indinavir', 'Lopinavir/Ritonavir','Nelfinavir', 'Ritonavir','Saquinavir','Tipranavir'), 'Drug_code'=c('ENV','MVC','DLV','EFV','ETR','NVP', 'ABC','ddl','FTC','3TC','d4T','TDF','ZDV/AZT', 'RAL','ATV','DRV','FPV', 'IDV','LPV/r','NFV','RTV','SQV','TPV')) #-------------------------------------------# # Creating unique patient identifiers # #-------------------------------------------# art$unique_id <- paste(art$site, art$patient, sep='.') basic$unique_id <- paste(basic$site, basic$patient, sep='.') follow$unique_id <- paste(follow$site, follow$patient, sep='.') cd4$unique_id <- paste(cd4$site, cd4$patient, sep='.') rna$unique_id <- paste(rna$site, rna$patient, sep='.') visit$unique_id <- paste(visit$site, visit$patient, sep='.') #----------------------------------------# # Setting 1900-01-01 dates to missing # #----------------------------------------# basic$birth_d <- ifelse(basic$birth_d == '1900-01-01',NA,basic$birth_d) basic$enrol_d <- ifelse(basic$enrol_d == '1900-01-01',NA,basic$enrol_d) follow$death_d <- ifelse(follow$death_d == '1900-01-01',NA,follow$death_d) visit$visit_d <- as.Date(visit$visit_d, format='%Y-%m-%d') visit$visit_d <- ifelse(visit$visit_d == '1900-01-01',NA,visit$visit_d) visit$visit_d <- as.Date(visit$visit_d, origin='1970-01-01') #---------------------------------------------------# # Removing subjects enrolled in a clinical trial # #---------------------------------------------------# ids.bf.clinical_trial <- length(unique(basic$unique_id)) clinical_trial_ids <- unique(art$unique_id[which(art$art_id %in% c('3TC,ABC,EFV,FTC,MRV/PBO,TDF', '3TC,AZT,EFV/PBO,MRV/PBO', 'DRV,ELV/PBO,RAL/PBO,RTV,TDF', 'EFV/PBO,FTC,RAL/PBO,TDF', 'ELV/PBO,FTC,LPV,RAL/PBO,RTV,TDF', 'EFV/PBO,FTC,RPV/PBO,TDF', 'FTC,MRV/PBO,RAL/PBO,TDF', 'LPV,RAL/PBO,RTV,TDF', 'DRV,ELV/PBO,ETR,RAL/PBO,RTV', 'EFV,FTC,MRV/PBO,TDF'))]) length.clinical.trial.ids <- length(clinical_trial_ids) ## Changed 6/18/13 -- changed back to include 6/9/14 art <- subset(art, unique_id %nin% clinical_trial_ids) # Basic also has an indicator 'clinicaltrial_y' that is incorrect for 2 IDs # Note: Last 3 from Mexico added 6/18/13 basic$clinicaltrial_y <- ifelse(basic$unique_id %in% c('argentina.1604','argentina.1808','brazil.25985','peru.1666','peru.901','mexico.157','mexico.431','mexico.682'),1, basic$clinicaltrial_y) # Hard coding this fix to honduras.194 regimen_class until I get new data (6/26/13) art$regimen_class[which(art$unique_id == 'honduras.194')] <- 'PI' # Saving information on those enrolled in a clinical trial determined either by clinicaltrial_y variable or as above based on regimen clinical_trial_ids <- unique(c(clinical_trial_ids,basic$unique_id[which(basic$clinicaltrial_y == 1)])) rct_basic.df <- subset(basic, unique_id %in% clinical_trial_ids) rct_art.df <- subset(art, unique_id %in% clinical_trial_ids) rct_follow.df <- subset(follow, unique_id %in% clinical_trial_ids) rct_visit.df <- subset(visit, unique_id %in% clinical_trial_ids) rct_cd4.df <- subset(cd4, unique_id %in% clinical_trial_ids) rct_rna.df <- subset(rna, unique_id %in% clinical_trial_ids) # Removing all who had clinicaltrial_y = 1 or 9 ids.w.unknown.clinical_trial_status <- length(which(basic$clinicaltrial_y == 9)) basic <- subset(basic, clinicaltrial_y == 0) art <- subset(art, unique_id %in% basic$unique_id) follow <- subset(follow, unique_id %in% basic$unique_id) visit <- subset(visit, unique_id %in% basic$unique_id) cd4 <- subset(cd4, unique_id %in% basic$unique_id) rna <- subset(rna, unique_id %in% basic$unique_id) ids.in.analysis <- length(unique(basic$unique_id)) # Calculating date of first haart (date_fhaart) art$art_sd <- as.Date(art$art_sd) art <- art[order(art$unique_id, art$art_sd),] art$art_ed <- as.Date(art$art_ed) art <- art[order(art$unique_id, art$art_sd),] date_first_regimen <- art[!duplicated(art$unique_id),c('unique_id','art_sd')] art$date_first_regimen <- date_first_regimen[match(art$unique_id, date_first_regimen$unique_id,nomatch=NA),'art_sd'] art$enrol_d <- basic[match(art$unique_id, basic$unique_id,nomatch=NA),'enrol_d'] art$enrol_d <- as.Date(art$enrol_d, format='%Y-%m-%d') # Not all regimens in 'art' are HAART regimens. Need to make sure grab the first HAART regimen. For those who # had ART before HAART, need to exclude them regimen_classes <- read.csv('/home/jenkinca/InfectiousDiseases/CathyMcGowan/3rdLineHaart/regimen-classes-06dec2012-ccm.csv',sep=',',header=TRUE,stringsAsFactors=FALSE,na.strings=c(NA,'')) art$regimen_class <- regimen_classes[match(art$art_id, regimen_classes$regimen,nomatch=NA),'class'] art <- art[order(art$unique_id,art$art_sd),] not_ART_naive <- art[!duplicated(art$unique_id),c('unique_id','regimen_class')] not_ART_naive$prior_ART <- ifelse(not_ART_naive$regimen_class == 'nonHAART',1,0) art$prior_art <- not_ART_naive[match(art$unique_id, not_ART_naive$unique_id,nomatch=NA),'prior_ART'] art$recart_y <- basic[match(art$unique_id,basic$unique_id,nomatch=NA),'recart_y'] art$birth_d <- basic[match(art$unique_id,basic$unique_id,nomatch=NA),'birth_d'] art$birth_d <- as.Date(art$birth_d,format=c('%Y-%m-%d')) #-----------------------------# # Subsetting based on # # incluson criteria # #-----------------------------# basic$prior_art <- art[match(basic$unique_id, art$unique_id,nomatch=NA),'prior_art'] basic$date_first_regimen <- art[match(basic$unique_id, art$unique_id, nomatch=NA),'date_first_regimen'] ids.w.recart_y.1 <- length(which(basic$recart_y == 1)) ids.w.recart_y.9 <- length(which(basic$recart_y == 9)) basic <- subset(basic, recart_y == 0) ids.after.recart_y <- dim(basic)[1] ids.w.prior_art.1 <- length(which(basic$prior_art == 1)) ids.w.no.art.data <- length(which(basic$unique_id %nin% art$unique_id)) ids.w.missing.regimen.not_ART_naive <- length(which(is.na(not_ART_naive$prior_ART) & not_ART_naive$unique_id %in% basic$unique_id)) basic <- subset(basic, prior_art == 0) ids.after.prior_art <- dim(basic)[1] ids.bf.jan12000 <- length(which(basic$date_first_regimen < '2000-01-01')) basic <- subset(basic, date_first_regimen >= '2000-01-01') ids.after.jan12000 <- dim(basic)[1] art <- subset(art, prior_art == 0 & date_first_regimen >= '2000-01-01' & recart_y == 0) date_fhaart <- art[!duplicated(art$unique_id),c('unique_id','art_sd')] # Setting date_fhaart art$date_fhaart <- date_fhaart[match(art$unique_id,date_fhaart$unique_id,nomatch=NA),'art_sd'] art$age_fhaart <- as.numeric(difftime(art$date_fhaart,art$birth_d,units='days'))/365.25 basic$age_fhaart <- art[match(basic$unique_id, art$unique_id, nomatch=NA),'age_fhaart'] ids.lt18 <- length(which(basic$age_fhaart < 18 | is.na(basic$age_fhaart))) basic <- subset(basic, age_fhaart >= 18) ids.after.lt18 <- dim(basic)[1] ids.w.missing.regimen <- unique(art$unique_id[which(is.na(art$regimen_class))]) basic <- subset(basic, unique_id %nin% ids.w.missing.regimen) art <- subset(art, unique_id %nin% ids.w.missing.regimen) ids.after.missing.regimen <- dim(basic)[1] ids.final <- ids.after.missing.regimen # Subsetting out pediatric patients; Also removing those with missing regimen classes (e.g., argentina.1121) art <- subset(art, age_fhaart >= 18 & !is.na(age_fhaart) & !is.na(regimen_class)) # Subsetting all other data sets to only include subjects in the art data set. follow <- subset(follow, unique_id %in% basic$unique_id) cd4 <- subset(cd4, unique_id %in% basic$unique_id) rna <- subset(rna, unique_id %in% basic$unique_id) visit <- subset(visit, unique_id %in% basic$unique_id) # Specifically, wants a categorical variable for NRTI type in model (AZT, d4T, TDF, Other) added to Table 1 cAZT<-function(a){ stuff<-ifelse(a=="AZT",1,0) sum(stuff) } cd4T<-function(a){ stuff<-ifelse(a=="D4T",1,0) sum(stuff) } cTDF<-function(a){ stuff<-ifelse(a=="TDF",1,0) sum(stuff) } # Only care about the art_id that was first HAART fhaart <- subset(art, art_sd == date_fhaart) fhaart$vazt<-sapply(lapply(strsplit(as.character(fhaart$art_id),','),cAZT),sum) fhaart$vd4t <- sapply(lapply(strsplit(as.character(fhaart$art_id),','),cd4T),sum) fhaart$vtdf <- sapply(lapply(strsplit(as.character(fhaart$art_id),','),cTDF),sum) # Defining those whose first HAART was an NNRTI nnrti_d <- subset(art, art_sd == date_first_regimen) nnrti_d <- subset(nnrti_d, regimen_class == 'NNRTI') nnrti_ids <- unique(nnrti_d$unique_id) #--------------------------------------------------# # Looking for year gap in viral load measurement # #--------------------------------------------------# rna$date_fhaart <- art[match(rna$unique_id, art$unique_id, nomatch=NA),'date_fhaart'] rna$rna_d <- as.Date(rna$rna_d, format='%Y-%m-%d') rna <- rna[order(rna$unique_id,rna$rna_d),] rna_gap <- lapply(split(subset(rna, select=c(unique_id,rna_d,rna_v,date_fhaart)),rna$unique_id),FUN=function(y){ y <- subset(y, rna_d >= date_fhaart & !is.na(rna_v)) if(dim(y)[1] > 1){ rna_d_lag <- c(y[1,'rna_d'],y[1:(dim(y)[1]-1),'rna_d']) date_diff <- as.numeric(y[,'rna_d'] - rna_d_lag) # For 1 year gap if(any(!is.na(date_diff) & date_diff > 365.25)){ which.row <- which(date_diff > 365.25)-1 year_gap_date <- min(y[which.row,'rna_d'],na.rm=TRUE) } else { year_gap_date <- NA } } else { year_gap_date <- NA } df <- data.frame('unique_id'=y[1,'unique_id'], 'year_gap_date'=year_gap_date) return(df) } ) rna_gap.df <- do.call('rbind',rna_gap) rna$year_gap_date <- rna_gap.df[match(rna$unique_id,rna_gap.df$unique_id,nomatch=NA),'year_gap_date'] rna$year_gap_date <- as.Date(rna$year_gap_date, origin='1970-01-01') no_w_year_rna_gap <- length(which(!is.na(rna$year_gap_date[!duplicated(rna$unique_id)]))) rna_new <- subset(rna, is.na(year_gap_date) | rna_d <= year_gap_date) # Adding to all other data sets basic$year_gap_date <- rna_new[match(basic$unique_id,rna_new$unique_id,nomatch=NA),'year_gap_date'] art$year_gap_date <- rna_new[match(art$unique_id,rna_new$unique_id,nomatch=NA),'year_gap_date'] visit$year_gap_date <- rna_new[match(visit$unique_id,rna_new$unique_id,nomatch=NA),'year_gap_date'] cd4$year_gap_date <- rna_new[match(cd4$unique_id,rna_new$unique_id,nomatch=NA),'year_gap_date'] follow$year_gap_date <- rna_new[match(follow$unique_id,rna_new$unique_id,nomatch=NA),'year_gap_date'] # Creating long data frame to find max date of follow-up and min date of event cd4$date <- cd4$cd4_d cd4$id_date <- paste(cd4$unique_id, cd4$date, sep='.') rna$date <- rna$rna_d rna$id_date <- paste(rna$unique_id, rna$date, sep='.') all_data <- find_max_date(art,cd4,rna,follow,basic) all_visits <- all_data[['all_visits']] art <- all_data[['art']] cd4 <- all_data[['cd4']] rna <- all_data[['rna']] follow <- all_data[['follow']] # # New way of defining AIDS at first HAART rather than AIDS at enrollment aids_at_fhaart <- aids_at_treatment_function(basic, art, cd4, visit) basic <- aids_at_fhaart[['basic']] art <- aids_at_fhaart[['art']] cd4 <- aids_at_fhaart[['cd4']] visit <- aids_at_fhaart[['visit']] aidsstatus <- aids_at_fhaart[['AIDSstatus']] aidsstatus$unique_id <- paste(aidsstatus$site, aidsstatus$patient,sep='.') all_visits$pre_artaids_clinical <- aidsstatus[match(all_visits$unique_id, aidsstatus$unique_id, nomatch=NA),'preARTAIDS.clinical'] all_visits$pre_artaids_all <- aidsstatus[match(all_visits$unique_id, aidsstatus$unique_id, nomatch=NA),'preARTAIDS.all'] # Using pre_artaids_clinical as aids_at_fhaart all_visits$aids_at_fhaart <- all_visits$pre_artaids_clinical all_visits$aids_at_fhaart.factor <- factor(all_visits$aids_at_fhaart, levels=c('not AIDS','AIDS'), labels=c('Not AIDS','AIDS')) # Adding in vazt, vd4t, and vtdf to all_visits all_visits$vazt <- fhaart[match(all_visits$unique_id, fhaart$unique_id, nomatch=NA),'vazt'] all_visits$vazt.factor <- factor(all_visits$vazt, levels=c(0,1), labels=c('No','Yes')) all_visits$vd4t <- fhaart[match(all_visits$unique_id, fhaart$unique_id, nomatch=NA),'vd4t'] all_visits$vd4t.factor <- factor(all_visits$vd4t, levels=c(0,1), labels=c('No','Yes')) all_visits$vtdf <- fhaart[match(all_visits$unique_id, fhaart$unique_id, nomatch=NA),'vtdf'] all_visits$vtdf.factor <- factor(all_visits$vtdf, levels=c(0,1), labels=c('No','Yes')) all_visits$gt_1_nrti <- apply(subset(all_visits, select=c(vazt,vd4t,vtdf)),1,FUN=function(y){if(sum(y) > 1){ 1 } else { 0 } }) all_visits$nrti_type <- ifelse(all_visits$vazt == 1 & all_visits$gt_1_nrti == 0,0, ifelse(all_visits$vd4t == 1 & all_visits$gt_1_nrti == 0,1, ifelse(all_visits$vtdf == 1 & all_visits$gt_1_nrti == 0,2, ifelse(all_visits$gt_1_nrti == 1 | (all_visits$vazt == 0 & all_visits$vd4t == 0 & all_visits$vtdf == 0),3,NA)))) all_visits$nrti_type.factor <- factor(all_visits$nrti_type, levels=c(0,1,2,3), labels=c('AZT','d4T','TDF','Other')) # Bryan wants initial HAART regimen class in Table 1 fhaart.df <- subset(art, date_fhaart == art_sd,select=c(unique_id,regimen_class)) all_visits$fhaart <- fhaart.df[match(all_visits$unique_id, fhaart.df$unique_id,nomatch=NA),'regimen_class'] all_visits$fhaart.factor <- factor(all_visits$fhaart) all_visits$site.factor <- factor(all_visits$site) all_visits$fu_ignoring_outcomes <- as.numeric(difftime(all_visits$max_date,all_visits$date_fhaart, units='days')/365.25) all_visits$fhaart_nnrti <- ifelse(all_visits$unique_id %in% nnrti_ids,1,0) all_visits$fu_nnrti <- NA all_visits$fu_nnrti[which(all_visits$fhaart_nnrti == 1)] <- all_visits$fu_ignoring_outcomes[which(all_visits$fhaart_nnrti == 1)] all_visits$transmission_risk <- as.factor(all_visits$transmission_risk) # Creating new transmission risk categorized variable (overwriting one created in find_max_date function) # Old one will be saved for models as transmission_risk_dichot all_visits$transmission_risk_dichot <- all_visits$transmission_risk_cat all_visits$transmission_risk_dichot.factor <- factor(all_visits$transmission_risk_cat, levels=c(0,1),labels=c('Heterosexual','Other')) all_visits$transmission_risk_cat <- ifelse(!is.na(all_visits$transmission_risk) & all_visits$transmission_risk == 'Heterosexual contact',0, ifelse(!is.na(all_visits$transmission_risk) & all_visits$transmission_risk %in% c('Bisexual','Homosexual contact'),1, ifelse(!is.na(all_visits$transmission_risk) & all_visits$transmission_risk %in% c('Heterosexual contact and Injecting drug user', 'Homo/Bisexual and Injecting drug user', 'Injecting drug user'),2, ifelse(!is.na(all_visits$transmission_risk) & all_visits$transmission_risk %in% c('Other (specify in mode_oth)', 'Perinatal', 'Transfusion, nonhemophilia related','Unknown'),3,NA)))) all_visits$transmission_risk_cat.factor <- factor(all_visits$transmission_risk_cat, levels=c(0,1,2,3), labels=c('Heterosexual','Homo/Bisexual','IDU','Other')) #----------------------------------------------# # Defining second line regimens # #----------------------------------------------# # First need to indicate which boosted PI switch indicates a switch to a second line regimen (requires switching the PI) art$regimen_class_new <- apply(subset(art, select=c(art_id,regimen_class)),1,FUN=function(y){ split_regimen <- unlist(strsplit(y['art_id'],',')) jk <- ifelse(y['regimen_class'] != 'OTHER',y['regimen_class'], ifelse(y['regimen_class'] == 'OTHER' & any(c('EFV','NVP','DLV','RPV','ETR') %in% split_regimen) & any(c('SQV','NFV','LPV','IDV','FPV','APV') %in% split_regimen) & 'RTV' %in% split_regimen,'NNRTI AND BOOSTED PI',y['regimen_class'])) }) # For those regimens listed as first HAART and are 'OTHER', need to be marked as incorrect art$incorrect_regimen <- ifelse(art$art_sd == art$date_fhaart & art$regimen_class == 'OTHER',1,0) previous_art_ed <- tapply(art$art_ed, INDEX=art$unique_id, FUN=function(y){ c(y[1],y[-length(y)]) }) previous_art_ed <- unlist(previous_art_ed) art$previous_art_ed <- previous_art_ed # art$previous_art_ed <- c(NA,art$art_ed[-nrow(art)]) art$previous_art_ed <- as.Date(art$previous_art_ed, origin='1970-01-01') art$month_gap <- ifelse(!is.na(art$previous_art_ed) & as.numeric(difftime(art$art_sd,art$previous_art_ed,units='days')) > 30,1, ifelse(is.na(art$previous_art_ed),NA,0)) art$month_gap_start <- ifelse(!is.na(art$month_gap) & art$month_gap == 1,art$previous_art_ed,NA) art$month_gap_end <- ifelse(!is.na(art$month_gap) & art$month_gap == 1,art$art_sd,NA) art$month_gap <- ifelse(is.na(art$month_gap),0,art$month_gap) # PI/r to PI/r in which backbone PI has been changed. A qualified OTHER regimen could be the 2nd line regimen. It qualifies for this if it contains an NNRTI and PI/r art <- art[order(art$unique_id,art$art_sd),] switch_pi <- lapply(split(subset(art, select=c(art_id,regimen_class_new,art_sd,art_ed,previous_art_ed,incorrect_regimen,unique_id)),art$unique_id),FUN=function(y){ if(!any(y[,'regimen_class_new'] %in% c('BOOSTED PI','NNRTI AND BOOSTED PI'))){ bpi_after_bpi <- rep(0, times=nrow(y)) second_line_start_date <- second_line_end_date <- NA } else { which.rows <- which(y[,'regimen_class_new'] %in% c('BOOSTED PI','NNRTI AND BOOSTED PI') & y[,'incorrect_regimen'] == 0) tmp.df <- y[which.rows,] if(nrow(tmp.df) == 1){ bpi_after_bpi <- 0 second_line_start_date <- second_line_end_date <- NA } else { bpi_after_bpi <- second_line_start_date <- second_line_end_date <- rep(NA,times=nrow(y)) for(j in 1:max(c(1,(nrow(tmp.df)-1)))){ tmp1 <- unlist(strsplit(tmp.df[j,'art_id'],',')) tmp2 <- unlist(strsplit(tmp.df[(j+1),'art_id'],',')) if(any(c('ATV','SQV','NFV','LPV','IDV','FPV','APV') %in% c(setdiff(tmp1,tmp2),setdiff(tmp2,tmp1)))){ bpi_after_bpi[which.rows[(j+1)]] <- 1 second_line_start_date[which.rows[(j+1)]] <- tmp.df[(j+1),'art_sd'] second_line_end_date[which.rows[(j+1)]] <- tmp.df[(j+1),'art_ed'] } else { bpi_after_bpi[which.rows[(j+1)]] <- 0 } } if(any(!is.na(bpi_after_bpi) & bpi_after_bpi == 1)){ min.bpi <- min(which(bpi_after_bpi == 1 & !is.na(bpi_after_bpi))) bpi_after_bpi[min.bpi:nrow(y)] <- 1 } if(nrow(y) > 1 & is.na(bpi_after_bpi[nrow(y)])){ bpi_after_bpi[nrow(y)] <- 0 second_line_start_date[nrow(y)] <- second_line_end_date[nrow(y)] <- NA } } } if(!all(is.na(second_line_start_date))){ min_second_line_start_date <- min(second_line_start_date,na.rm=TRUE) } else { min_second_line_start_date <- NA } df <- data.frame('unique_id'=y[,'unique_id'], 'bpi_after_bpi'=bpi_after_bpi, 'second_line_start_date'=second_line_start_date, 'second_line_end_date'=second_line_end_date, 'min_second_line_start_date'=min_second_line_start_date) return(df) }) switch_pi.df <- do.call('rbind',switch_pi) art <- art[order(art$unique_id,art$art_sd),] art <- cbind(art,'bpi_after_bpi'=switch_pi.df$bpi_after_bpi) art <- cbind(art, 'second_line_start_date_bpi'=switch_pi.df$second_line_start_date) art$second_line_start_date_bpi <- as.Date(art$second_line_start_date_bpi, origin='1970-01-01') art <- cbind(art, 'second_line_end_date_bpi'=switch_pi.df$second_line_end_date) art$second_line_end_date_bpi <- as.Date(art$second_line_end_date_bpi, origin='1970-01-01') art <- cbind(art, 'min_second_line_start_date_bpi'=switch_pi.df$min_second_line_start_date) art$bpi_after_bpi2 <- ifelse(is.na(art$bpi_after_bpi),0,art$bpi_after_bpi) art$bpi_after_bpi <- art$bpi_after_bpi2 # Now finding those with switch from NNRTI to boosted PI (in secondary analysis, need to subset on those whose first HAART was NNRTI and then check for change to second line regimen of PI/r) art <- art[order(art$unique_id,art$art_sd),] switch_nnrti <- lapply(split(subset(art, select=c(art_id,regimen_class_new,art_sd,art_ed,unique_id)),art$unique_id),FUN=function(y){ nnrti_switch <- second_line_start_date <- second_line_end_date <- rep(NA,nrow(y)) if(!any(y[,'regimen_class_new'] == 'NNRTI')){ bpi_after_nnrti <- rep(0, times=nrow(y)) } else { which.rows <- which(y[,'regimen_class_new'] == 'NNRTI') tmp.df <- y[min(which.rows):nrow(y),] if(!any(tmp.df[,'regimen_class_new'] %in% c('BOOSTED PI','NNRTI AND BOOSTED PI'))){ bpi_after_nnrti <- 0 } else { if(nrow(tmp.df) == 1){ bpi_after_nnrti <- 0 } else { which.bpi <- min(which(tmp.df[,'regimen_class_new'] %in% c('BOOSTED PI','NNRTI AND BOOSTED PI'))) bpi_after_nnrti <- second_line_start_date <- second_line_end_date <- rep(NA,times=nrow(tmp.df)) bpi_after_nnrti[which.bpi:nrow(tmp.df)] <- 1 second_line_start_date[which.bpi] <- tmp.df[which.bpi,'art_sd'] second_line_end_date[which.bpi] <- tmp.df[which.bpi,'art_ed'] if(min(which.rows) > 1){ bpi_after_nnrti <- c(rep(NA,times=length(seq(1,(min(which.rows)-1)))),bpi_after_nnrti) second_line_start_date <- c(rep(NA,times=length(seq(1,(min(which.rows)-1)))),second_line_start_date) second_line_end_date <- c(rep(NA,times=length(seq(1,(min(which.rows)-1)))),second_line_end_date) } } } } if(!all(is.na(second_line_start_date))){ min_second_line_start_date <- min(second_line_start_date,na.rm=TRUE) } else { min_second_line_start_date <- NA } df <- data.frame('unique_id'=y[,'unique_id'], 'bpi_after_nnrti'=bpi_after_nnrti, 'second_line_start_date'=second_line_start_date, 'second_line_end_date'=second_line_end_date, 'min_second_line_start_date'=min_second_line_start_date) return(df) }) switch_nnrti.df <- do.call('rbind',switch_nnrti) art <- art[order(art$unique_id,art$art_sd),] art <- cbind(art, 'bpi_after_nnrti'=switch_nnrti.df$bpi_after_nnrti) art <- cbind(art, 'second_line_start_date_nnrti'=switch_nnrti.df$second_line_start_date) art$second_line_start_date_nnrti <- as.Date(art$second_line_start_date_nnrti, origin='1970-01-01') art <- cbind(art, 'second_line_end_date_nnrti'=switch_nnrti.df$second_line_end_date) art$second_line_end_date_nnrti <- as.Date(art$second_line_end_date_nnrti, origin='1970-01-01') art <- cbind(art, 'min_second_line_start_date_nnrti'=switch_nnrti.df$min_second_line_start_date) art$bpi_after_nnrti2 <- ifelse(is.na(art$bpi_after_nnrti),0,art$bpi_after_nnrti) art$bpi_after_nnrti2 <- ifelse(is.na(art$bpi_after_nnrti),0,art$bpi_after_nnrti) art$bpi_after_nnrti <- art$bpi_after_nnrti2 art$bpi_after_bpi_or_nnrti <- ifelse(art$bpi_after_bpi == 1 | art$bpi_after_nnrti == 1,1,0) art$min_second_line_start_date_bpi_or_nnrti <- apply(subset(art, select=c(min_second_line_start_date_bpi,min_second_line_start_date_nnrti)),1,FUN=min,na.rm=TRUE) art$min_second_line_start_date_bpi_or_nnrti <- ifelse(art$min_second_line_start_date_bpi_or_nnrti == Inf,NA,art$min_second_line_start_date_bpi_or_nnrti) art$min_second_line_start_date_bpi_or_nnrti <- as.Date(art$min_second_line_start_date_bpi_or_nnrti,origin='1970-01-01') # Adding to all_visits all_visits$min_second_line_start_date_bpi_or_nnrti <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'min_second_line_start_date_bpi_or_nnrti'] # Adding in date second line regimen date information to rna data to check virologic failure on all viral loads from the first second line regimen forward rna$min_second_line_start_date_nnrti <- art[match(rna$unique_id,art$unique_id,nomatch=NA),'min_second_line_start_date_nnrti'] rna$min_second_line_start_date_nnrti <- as.Date(rna$min_second_line_start_date_nnrti, origin='1970-01-01') rna$min_second_line_start_date_bpi_or_nnrti <- art[match(rna$unique_id,art$unique_id,nomatch=NA),'min_second_line_start_date_bpi_or_nnrti'] rna$min_second_line_start_date_bpi_or_nnrti <- as.Date(rna$min_second_line_start_date_bpi_or_nnrti, origin='1970-01-01') # Need to be able to map in rna data when someone was on or not on a BPI or approved 'OTHER regimen after the switch rna$bpi_after_bpi <- rna$bpi_after_nnrti <- rna$bpi_after_bpi_or_nnrti <- 0 art.tmp <- subset(art, select=c(unique_id,art_sd,art_ed,bpi_after_bpi,bpi_after_nnrti,bpi_after_bpi_or_nnrti)) # Just for now entering in today for end dates that are missing art.tmp$art_ed[which(is.na(art.tmp$art_ed))] <- '2013-06-14' for(i in 1:nrow(rna)){ # For bpi_after_bpi indicator loop.tmp <- art.tmp[which(art.tmp$unique_id == rna$unique_id[i]),] loop.tmp$art_sd[which(loop.tmp$bpi_after_bpi == 0)] <- NA loop.tmp$art_ed[which(loop.tmp$bpi_after_bpi == 0)] <- NA inside <- !is.na(loop.tmp$art_sd) & loop.tmp$art_sd <= rna$rna_d[i] & loop.tmp$art_ed > rna$rna_d[i] if(any(inside)){ rna$bpi_after_bpi[i] <- loop.tmp$bpi_after_bpi[inside] } # For bpi_after_nnrti indicator loop.tmp <- art.tmp[which(art.tmp$unique_id == rna$unique_id[i]),] loop.tmp$art_sd[which(loop.tmp$bpi_after_nnrti == 0)] <- NA loop.tmp$art_ed[which(loop.tmp$bpi_after_nnrti == 0)] <- NA inside <- !is.na(loop.tmp$art_sd) & loop.tmp$art_sd <= rna$rna_d[i] & loop.tmp$art_ed > rna$rna_d[i] if(any(inside)){ rna$bpi_after_nnrti[i] <- loop.tmp$bpi_after_nnrti[inside] } } rna$bpi_after_bpi_or_nnrti <- ifelse(rna$bpi_after_bpi == 1 | rna$bpi_after_nnrti == 1,1,0) # Creating overall indicators for those who switched to second line with the PI/r to PI/r definition art$overall_second_line_bpi <- NA split(art$overall_second_line_bpi,art$unique_id) <- lapply(split(art$bpi_after_bpi,art$unique_id),FUN=function(y){ max(y) }) # Creating overall indicators for those who switched to second line with the NNRTI to PI/r definition art$overall_second_line_nnrti <- NA split(art$overall_second_line_nnrti,art$unique_id) <- lapply(split(art$bpi_after_nnrti,art$unique_id),FUN=function(y){ max(y) }) # Creating overall indicators for those who switched to second line with the PI/r to PI/r or NNRTI to PI/r definition art$overall_second_line <- NA split(art$overall_second_line,art$unique_id) <- lapply(split(art$bpi_after_bpi_or_nnrti,art$unique_id),FUN=function(y){ max(y) }) #--------------------------------------# # Assessing virologic failure # #--------------------------------------# # Either criteria for switch (BPI to BPI or NNRTI to BPI) vl_failure_all <- first_failure_virologic_failure_function(switch='all',rna,detectable_level=400) # Bryan said fail.1000 is the overall one now, not fail. art$overall_first_failure.1000 <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'vl_failure.1000'] art$overall_first_failure.1000 <- ifelse(is.na(art$overall_first_failure.1000),0,art$overall_first_failure.1000) art$overall_first_failure.1000.factor <- factor(art$overall_first_failure.1000, levels=c(0,1),labels=c('No','Yes')) art$days_from_fhaart_to_vl_failure.1000 <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure.1000'] art$date_first_failure.1000 <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'date_vl_failure.1000'] all_visits$date_first_failure.1000 <- vl_failure_all[match(all_visits$unique_id,vl_failure_all$unique_id,nomatch=NA),'date_vl_failure.1000'] all_visits$days_from_fhaart_to_vl_failure.1000 <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure.1000'] art$overall_first_failure <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'vl_failure'] art$overall_first_failure <- ifelse(is.na(art$overall_first_failure),0,art$overall_first_failure) art$overall_first_failure.factor <- factor(art$overall_first_failure, levels=c(0,1),labels=c('No','Yes')) art$days_from_fhaart_to_vl_failure <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure'] art$date_first_failure <- vl_failure_all[match(art$unique_id,vl_failure_all$unique_id,nomatch=NA),'date_vl_failure'] all_visits$date_first_failure <- vl_failure_all[match(all_visits$unique_id,vl_failure_all$unique_id,nomatch=NA),'date_vl_failure'] all_visits$days_from_fhaart_to_vl_failure <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure'] art$overall_first_failure.1000_new <- ifelse(!is.na(art$year_gap_date) & !is.na(art$date_first_failure.1000) & art$date_first_failure.1000 > art$year_gap_date,0, art$overall_first_failure.1000) art$overall_first_failure.1000_new.factor <- factor(art$overall_first_failure.1000_new, levels=c(0,1), labels=c('No','Yes')) art$date_first_failure.1000_new <- ifelse(!is.na(art$date_first_failure.1000) & !is.na(art$year_gap_date) & art$date_first_failure.1000 > art$year_gap_date,NA,art$date_first_failure.1000) art$date_first_failure.1000_new <- as.Date(art$date_first_failure.1000_new,origin='1970-01-01') art$days_from_fhaart_to_vl_failure.1000_new <- ifelse(is.na(art$date_first_failure.1000_new),NA,art$days_from_fhaart_to_vl_failure.1000) all_visits$date_first_failure.1000_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'date_first_failure.1000_new'] all_visits$days_from_fhaart_to_vl_failure.1000_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure.1000_new'] art$overall_first_failure_new <- ifelse(!is.na(art$year_gap_date) & !is.na(art$date_first_failure) & art$date_first_failure > art$year_gap_date,0, art$overall_first_failure) art$overall_first_failure_new.factor <- factor(art$overall_first_failure_new, levels=c(0,1), labels=c('No','Yes')) art$date_first_failure_new <- ifelse(!is.na(art$date_first_failure) & !is.na(art$year_gap_date) & art$date_first_failure > art$year_gap_date,NA,art$date_first_failure) art$date_first_failure_new <- as.Date(art$date_first_failure_new,origin='1970-01-01') art$days_from_fhaart_to_vl_failure_new <- ifelse(is.na(art$date_first_failure_new),NA,art$days_from_fhaart_to_vl_failure) all_visits$date_first_failure_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'date_first_failure_new'] all_visits$days_from_fhaart_to_vl_failure_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_vl_failure_new'] # Checking to see how many failures during month gap art$failure_during_gap <- apply(subset(art, select=c(month_gap,month_gap_start,month_gap_end,date_first_failure.1000,date_first_failure)),1,FUN=function(y){ if(y['month_gap'] == 1){ if(!is.na(y['date_first_failure.1000']) & y['date_first_failure.1000'] >= y['month_gap_start'] & y['date_first_failure.1000'] <= y['month_gap_end']){ failure_during_gap <- 1 } else { failure_during_gap <- 0 } } else { failure_during_gap <- NA } }) art$failure_on_art.1000 <- ifelse(art$date_first_failure.1000 >= art$art_sd & art$date_first_failure.1000 <= art$art_ed,1,0) stuff <- table(art$failure_on_art.1000,art$month_gap) #----------------------------------# # Adding in info to ART data # #----------------------------------# art$clinicaltrial_y <- basic[match(art$unique_id,basic$unique_id,nomatch=NA),'clinicaltrial_y'] art$male <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'male'] art$male.factor <- factor(art$male, levels=c(0,1), labels=c('Female','Male')) art$baseline_clinical_stage_new <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'baseline_clinical_stage_new'] art$aids_at_fhaart <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'aids_at_fhaart'] art$aids_at_fhaart.factor <- factor(art$aids_at_fhaart, levels=c(0,1),labels=c('Not AIDS','AIDS')) art$baseline_cd4 <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'baseline_cd4'] art$baseline_log_rna <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'baseline_log_rna'] art$site.factor <- factor(art$site) art$overall_death_indicator.factor <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'overall_death_indicator.factor'] all_visits$overall_death_indicator_new <- ifelse(!is.na(all_visits$year_gap_date) & !is.na(all_visits$death_d) & all_visits$death_d > all_visits$year_gap_date,0, all_visits$overall_death_indicator) all_visits$overall_death_indicator_new.factor <- factor(all_visits$overall_death_indicator_new, levels=c(0,1), labels=c('No','Yes')) art$fhaart <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'fhaart'] art$fhaart.factor <- factor(art$fhaart) art$transmission_risk_cat.factor <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'transmission_risk_cat.factor'] art$transmission_risk <- all_visits[match(art$unique_id,all_visits$unique_id,nomatch=NA),'transmission_risk'] art$transmission_risk <- factor(art$transmission_risk) art$year_gap_date <- rna[match(art$unique_id, rna$unique_id,nomatch=NA),'year_gap_date'] art$composite_outcome <- ifelse(art$overall_second_line == 1 | art$overall_first_failure == 1,1,0) art$composite_outcome.factor <- factor(art$composite_outcome, levels=c(0,1), labels=c('No','Yes')) art$composite_outcome.1000 <- ifelse(art$overall_second_line == 1 | art$overall_first_failure.1000 == 1,1,0) art$composite_outcome.1000.factor <- factor(art$composite_outcome.1000, levels=c(0,1), labels=c('No','Yes')) art$days_from_fhaart_to_second_line_start <- as.numeric(difftime(art$min_second_line_start_date_bpi_or_nnrti,art$date_fhaart,units='days')) all_visits$days_from_fhaart_to_second_line_start <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_second_line_start'] all_visits$overall_second_line <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_second_line'] all_visits$overall_first_failure <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_first_failure'] all_visits$overall_first_failure.1000 <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_first_failure.1000'] art$min_second_line_start_date_bpi_or_nnrti_new <- ifelse(!is.na(art$year_gap_date) & !is.na(art$min_second_line_start_date_bpi_or_nnrti) & art$min_second_line_start_date_bpi_or_nnrti > art$year_gap_date,NA,art$min_second_line_start_date_bpi_or_nnrti) art$min_second_line_start_date_bpi_or_nnrti_new <- as.Date(art$min_second_line_start_date_bpi_or_nnrti_new, origin='1970-01-01') all_visits$min_second_line_start_date_bpi_or_nnrti_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'min_second_line_start_date_bpi_or_nnrti_new'] art$overall_second_line_new <- ifelse(!is.na(art$year_gap_date) & !is.na(art$min_second_line_start_date_bpi_or_nnrti) & art$min_second_line_start_date_bpi_or_nnrti > art$year_gap_date,0,art$overall_second_line) art$overall_second_line_new.factor <- factor(art$overall_second_line_new, levels=c(0,1), labels=c('No','Yes')) art$composite_outcome_new <- ifelse(art$overall_second_line_new == 1 | art$overall_first_failure_new == 1,1,0) art$composite_outcome_new.factor <- factor(art$composite_outcome_new, levels=c(0,1), labels=c('No','Yes')) art$composite_outcome.1000_new <- ifelse(art$overall_second_line_new == 1 | art$overall_first_failure.1000_new == 1,1,0) art$composite_outcome.1000_new.factor <- factor(art$composite_outcome.1000_new, levels=c(0,1), labels=c('No','Yes')) art$days_from_fhaart_to_second_line_start_new <- as.numeric(difftime(art$min_second_line_start_date_bpi_or_nnrti_new,art$date_fhaart,units='days')) all_visits$days_from_fhaart_to_second_line_start_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'days_from_fhaart_to_second_line_start_new'] all_visits$overall_second_line_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_second_line_new'] all_visits$overall_first_failure_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_first_failure_new'] all_visits$overall_first_failure.1000_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'overall_first_failure.1000_new'] all_visits$composite_outcome_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome_new'] all_visits$composite_outcome_new.factor <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome_new.factor'] all_visits$composite_outcome.1000_new <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome.1000_new'] all_visits$composite_outcome.1000_new.factor <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome.1000_new.factor'] # tmp_2nd_line <- subset(art, clinicaltrial_y == 0 & !duplicated(unique_id)) art$pre_artaids_all <- all_visits[match(art$unique_id, all_visits$unique_id, nomatch=NA),'pre_artaids_all'] art$aids_at_fhaart <- all_visits[match(art$unique_id, all_visits$unique_id, nomatch=NA),'aids_at_fhaart'] art$aids_at_fhaart.factor <- all_visits[match(art$unique_id, all_visits$unique_id, nomatch=NA),'aids_at_fhaart.factor'] tmp_2nd_line <- subset(art, !duplicated(unique_id)) tmp_2nd_line$overall_second_line.factor <- factor(tmp_2nd_line$overall_second_line, levels=c(0,1), labels=c('No 2nd line regimen','2nd line regimen')) # tmp_2nd_line_b <- subset(tmp_2nd_line, overall_failure_on_second_line.1000 == 1) # Based on tmp_2nd_line, there are 9501 subjects not on a clinical trial. # Out of the 9501, table(tmp_2nd_line$overall_second_line[!duplicated(tmp_2nd_line$unique_id)]) # gives 1160 subjects who started a 2nd line regimen (12.2%); # of the 1160, 219 failed a 2nd regimen (ignoring if they started a 3rd line regimen) but 200 of the 1160 had no # HIV1-RNA measurements in the rna data. 9466 who were not in a clinical trial also did not start a 3rd line regimen. # 1135 of those who were not in a clinical trial and did not start a 3rd line regimen started a 2nd line regimen. # 208 of those who did not start a 3rd line regimen failed the 2nd line regimen (208 out of 9466 or 208 out of 1135 # or 18.3%. tbl_by_tbl <- with(tmp_2nd_line, table(overall_first_failure.1000.factor,overall_second_line.factor)) tmp_2nd_line_no_haiti <- subset(tmp_2nd_line, site.factor != 'haiti') tmp_2nd_line_no_haiti$site.factor <- factor(tmp_2nd_line_no_haiti$site) #------------------------------------------------------------------------------------------------# # Time-to-event analysis # #------------------------------------------------------------------------------------------------# # Need to add appropriate variables to all_visits in order to calculate the first endpoint all_visits$composite_outcome.1000 <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome.1000'] all_visits$composite_outcome <- art[match(all_visits$unique_id,art$unique_id,nomatch=NA),'composite_outcome'] all_visits$max_date_outcome <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti, date_first_failure, max_date)),1,FUN=function(y){ min(y,na.rm=TRUE) }) all_visits$max_date_outcome.1000 <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti, date_first_failure.1000, max_date)),1,FUN=function(y){ min(y,na.rm=TRUE) }) all_visits$max_date_outcome <- as.Date(all_visits$max_date_outcome, format='%Y-%m-%d') all_visits$max_date_outcome.1000 <- as.Date(all_visits$max_date_outcome.1000, format='%Y-%m-%d') all_visits$year_gap_date <- rna[match(all_visits$unique_id,rna$unique_id,nomatch=NA),'year_gap_date'] all_visits$max_date_outcome_new <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti_new, date_first_failure_new, max_date, year_gap_date)),1,FUN=function(y){ if(is.na(y['year_gap_date'])){ min(y[1:3],na.rm=TRUE) } else if(!is.na(y['year_gap_date'])){ min(y,na.rm=TRUE) } }) all_visits$max_date_outcome_new <- as.Date(all_visits$max_date_outcome_new, format='%Y-%m-%d') all_visits$max_date_outcome.1000_new <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti_new, date_first_failure.1000_new, max_date, year_gap_date)),1,FUN=function(y){ if(is.na(y['year_gap_date'])){ min(y[1:3],na.rm=TRUE) } else if(!is.na(y['year_gap_date'])){ min(y,na.rm=TRUE) } }) all_visits$max_date_outcome.1000_new <- as.Date(all_visits$max_date_outcome.1000_new, format='%Y-%m-%d') all_visits$min_date_of_evt <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti,date_first_failure,year_gap_date)),1,FUN=function(y){ if(all(is.na(y[1:2]))){ NA } else { min(y[1:2],na.rm=TRUE) } }) all_visits$min_date_of_evt.1000 <- apply(subset(all_visits, select=c(min_second_line_start_date_bpi_or_nnrti,date_first_failure.1000)),1,FUN=function(y){ if(all(is.na(y))){ NA } else { min(y,na.rm=TRUE) } }) all_visits$min_date_of_evt <- as.Date(all_visits$min_date_of_evt, format=c('%Y-%m-%d')) all_visits$min_date_of_evt.1000 <- as.Date(all_visits$min_date_of_evt.1000, format=c('%Y-%m-%d')) # Calculating time to event from date_fhaart to ... should this be min_date_of_evt or max_date_outcome??? all_visits$time_to_evt <- as.numeric(difftime(all_visits$min_date_of_evt,all_visits$date_fhaart,units='days')/365.25) all_visits$time_to_evt.1000 <- as.numeric(difftime(all_visits$min_date_of_evt.1000,all_visits$date_fhaart,units='days')/365.25) # all_visits_original$time_to_evt.1000 <- as.numeric(difftime(all_visits_original$min_date_of_evt.1000,all_visits_original$date_fhaart,units='days')/365.25) # Adding in an indicator for whether someone's initial HAART was an NNRTI or not all_visits$initial_haart_nnrti <- ifelse(all_visits$unique_id %in% nnrti_ids,1,0) all_visits$initial_haart_nnrti.factor <- factor(all_visits$initial_haart_nnrti, levels=c(1,0),labels=c('NNRTI','Other')) # Adding in ltfu information all_visits$lost1 <- ltfu[match(all_visits$unique_id,ltfu$unique_id,nomatch=NA),'lost1'] # Adding in regimen class information to all_visits to then categorize for the Cox models sub_art <- subset(art, date_fhaart == art_sd) sub_art$fhaart_regimen_class <- ifelse(sub_art$regimen_class == 'NNRTI',0, ifelse(sub_art$regimen_class %in% c('BOOSTED PI'),1, ifelse(!is.na(sub_art$regimen_class),2,NA))) sub_art$fhaart_regimen_class.factor <- factor(sub_art$fhaart_regimen_class, levels=c(0,1,2), labels=c('NNRTI','PI','OTHER')) all_visits$fhaart_regimen_class.factor <- sub_art[match(all_visits$unique_id,sub_art$unique_id,nomatch=NA),'fhaart_regimen_class.factor'] # Subsetting on records <= max_date_outcome d <- subset(all_visits, !duplicated(unique_id) & date <= max_date_outcome) d.1000 <- subset(all_visits, !duplicated(unique_id) & date <= max_date_outcome.1000) # d.1000_original <- subset(all_visits_original, !duplicated(unique_id) & date <= max_date_outcome.1000) d$upper_limit_of_fu <- apply(subset(d, select=c(min_date_of_evt, max_date)),1,FUN=function(y){ min(y,na.rm=TRUE) }) d.1000$upper_limit_of_fu.1000 <- apply(subset(d.1000, select=c(min_date_of_evt.1000, max_date)),1,FUN=function(y){ min(y,na.rm=TRUE) }) d$total_followup <- as.numeric(difftime(d$max_date_outcome,d$date_fhaart,units='days')/365.25) d.1000$total_followup.1000 <- as.numeric(difftime(d.1000$max_date_outcome.1000,d.1000$date_fhaart,units='days')/365.25) d$fu_ignoring_outcomes <- as.numeric(difftime(d$max_date, d$date_fhaart, units='days')/365.25) d.1000$fu_ignoring_outcomes <- as.numeric(difftime(d.1000$max_date, d.1000$date_fhaart, units='days')/365.25) # Need follow-up time for virologic failure d$total_followup_failure <- ifelse(!is.na(as.numeric(d$days_from_fhaart_to_vl_failure)),as.numeric(d$days_from_fhaart_to_vl_failure)/365.25, d$fu_ignoring_outcomes) d.1000$total_followup_failure.1000 <- ifelse(!is.na(as.numeric(d.1000$days_from_fhaart_to_vl_failure.1000)),as.numeric(d.1000$days_from_fhaart_to_vl_failure.1000)/365.25, d.1000$fu_ignoring_outcomes) d$total_followup_reg_change <- ifelse(!is.na(d$days_from_fhaart_to_second_line_start),d$days_from_fhaart_to_second_line_start/365.25, d$fu_ignoring_outcomes) d.1000$total_followup_reg_change <- ifelse(!is.na(d.1000$days_from_fhaart_to_second_line_start),d.1000$days_from_fhaart_to_second_line_start/365.25, d.1000$fu_ignoring_outcomes) d$total_followup_new <- as.numeric(difftime(d$max_date_outcome_new,d$date_fhaart,units='days')/365.25) d.1000$total_followup.1000_new <- as.numeric(difftime(d.1000$max_date_outcome.1000_new,d.1000$date_fhaart,units='days')/365.25) d$total_followup_failure_new <- ifelse(!is.na(as.numeric(d$days_from_fhaart_to_vl_failure_new)),as.numeric(d$days_from_fhaart_to_vl_failure_new)/365.25, ifelse(!is.na(d$year_gap_date),as.numeric((d$year_gap_date - d$date_fhaart)/365.25), d$fu_ignoring_outcomes)) d.1000$total_followup_failure.1000_new <- ifelse(!is.na(as.numeric(d.1000$days_from_fhaart_to_vl_failure.1000_new)),as.numeric(d.1000$days_from_fhaart_to_vl_failure.1000_new)/365.25, ifelse(!is.na(d.1000$year_gap_date),as.numeric((d.1000$year_gap_date - d.1000$date_fhaart)/365.25), d.1000$fu_ignoring_outcomes)) d$total_followup_reg_change_new <- ifelse(!is.na(d$days_from_fhaart_to_second_line_start_new),d$days_from_fhaart_to_second_line_start_new/365.25, ifelse(!is.na(d$year_gap_date),as.numeric((d$year_gap_date - d$date_fhaart)/365.25),d$fu_ignoring_outcomes)) d.1000$total_followup_reg_change_new <- ifelse(!is.na(d.1000$days_from_fhaart_to_second_line_start_new),d.1000$days_from_fhaart_to_second_line_start_new/365.25, ifelse(!is.na(d.1000$year_gap_date), as.numeric((d.1000$year_gap_date - d.1000$date_fhaart)/365.25),d.1000$fu_ignoring_outcomes)) # Need to also include year of HAART start as a covariate # Ignoring 1000 criterion for defining failure on 2nd regimen year. <- apply(subset(d, select=c(date_fhaart)), 1, FUN=function(y){ unlist(strsplit(as.character(y), split='-'))[1] }) d$year_fhaart <- as.numeric(year.) # Including 1000 criterion for defining failure on 2nd regimen year. <- apply(subset(d.1000, select=c(date_fhaart)), 1, FUN=function(y){ unlist(strsplit(as.character(y), split='-'))[1] }) d.1000$year_fhaart <- as.numeric(year.) # Removing those who had total_followup == 0.0 (ie., just record of age_fhaart) no.subjects.w.0.fu <- dim(d.1000[which(d.1000$site.factor != 'haiti'),])[1] d <- subset(d, total_followup > 0) d.1000 <- subset(d.1000, total_followup.1000 > 0) # Cumulative incidence # Need to create a 3-level censoring variable that is 0=No event, 1=Third line regimen or failure of 2nd line, 2=Death d$censor <- ifelse(d$composite_outcome == 0 & d$overall_death_indicator == 0,0, ifelse(d$composite_outcome == 1,1, ifelse(d$composite_outcome == 0 & d$overall_death_indicator == 1,2,NA))) d$censor_failure <- ifelse(d$overall_first_failure == 0 & d$overall_death_indicator == 0,0, ifelse(d$overall_first_failure == 1,1, ifelse(d$overall_first_failure == 0 & d$overall_death_indicator == 1,2,NA))) d$censor_reg_change <- ifelse(d$overall_second_line == 0 & d$overall_death_indicator == 0,0, ifelse(d$overall_second_line == 1,1, ifelse(d$overall_second_line == 0 & d$overall_death_indicator == 1,2,NA))) d.1000$censor.1000 <- ifelse(d.1000$composite_outcome.1000 == 0 & d.1000$overall_death_indicator == 0,0, ifelse(d.1000$composite_outcome.1000 == 1,1, ifelse(d.1000$composite_outcome.1000 == 0 & d.1000$overall_death_indicator == 1,2,NA))) d.1000$censor_failure.1000 <- ifelse(d.1000$overall_first_failure.1000 == 0 & d.1000$overall_death_indicator == 0,0, ifelse(d.1000$overall_first_failure.1000 == 1,1, ifelse(d.1000$overall_first_failure.1000 == 0 & d.1000$overall_death_indicator == 1,2,NA))) d.1000$censor_reg_change <- ifelse(d.1000$overall_second_line == 0 & d.1000$overall_death_indicator == 0,0, ifelse(d.1000$overall_second_line == 1,1, ifelse(d.1000$overall_second_line == 0 & d.1000$overall_death_indicator == 1,2,NA))) d$censor_new <- ifelse(d$composite_outcome_new == 0 & d$overall_death_indicator_new == 0,0, ifelse(d$composite_outcome_new == 1,1, ifelse(d$composite_outcome_new == 0 & d$overall_death_indicator_new == 1,2,NA))) d$censor_failure_new <- ifelse(d$overall_first_failure_new == 0 & d$overall_death_indicator_new == 0,0, ifelse(d$overall_first_failure_new == 1,1, ifelse(d$overall_first_failure_new == 0 & d$overall_death_indicator_new == 1,2,NA))) d$censor_reg_change_new <- ifelse(d$overall_second_line_new == 0 & d$overall_death_indicator_new == 0,0, ifelse(d$overall_second_line_new == 1,1, ifelse(d$overall_second_line_new == 0 & d$overall_death_indicator_new == 1,2,NA))) d.1000$censor.1000_new <- ifelse(d.1000$composite_outcome.1000_new == 0 & d.1000$overall_death_indicator_new == 0,0, ifelse(d.1000$composite_outcome.1000_new == 1,1, ifelse(d.1000$composite_outcome.1000_new == 0 & d.1000$overall_death_indicator_new == 1,2,NA))) d.1000$censor_failure.1000_new <- ifelse(d.1000$overall_first_failure.1000_new == 0 & d.1000$overall_death_indicator_new == 0,0, ifelse(d.1000$overall_first_failure.1000_new == 1,1, ifelse(d.1000$overall_first_failure.1000_new == 0 & d.1000$overall_death_indicator_new == 1,2,NA))) d.1000$censor_reg_change_new <- ifelse(d.1000$overall_second_line_new == 0 & d.1000$overall_death_indicator_new == 0,0, ifelse(d.1000$overall_second_line_new == 1,1, ifelse(d.1000$overall_second_line_new == 0 & d.1000$overall_death_indicator_new == 1,2,NA))) # Table 1 summary tbl1.summ <- summary(site.factor ~ male.factor + age_fhaart + transmission_risk_w_unknown_level.factor + pre_artaids_all + aids_at_fhaart.factor + baseline_cd4 + baseline_log_rna + overall_death_indicator.factor + fhaart.factor + nrti_type.factor + fu_ignoring_outcomes + fu_nnrti, data=all_visits,subset=!duplicated(unique_id) & fu_ignoring_outcomes > 0,method='reverse',overall=TRUE) tbl1.summ$labels <- c('Sex','Age of first HAART', 'Possible route of infection', 'AIDS at first HAART (all)', 'AIDS at first HAART (clinical)', 'Baseline CD4','Baseline log viral load', 'Death','Initial HAART','NRTI in initial HAART', 'Follow-up, ignoring outcomes','Follow-up for NNRTI first HAART, ignoring outcomes') # latex(tbl1.summ, file='', caption='Description of the cohort by site.',exclude1=FALSE,long=TRUE,digits=2,where='!h',longtable=TRUE,landscape=TRUE,size='tiny') # Deleting Haiti from main data d.1000$composite_outcome.1000.factor <- factor(d.1000$composite_outcome.1000, levels=c(0,1), labels=c('No','Yes')) d.1000$overall_first_failure.1000.factor <- factor(d.1000$overall_first_failure.1000, levels=c(0,1), labels=c('No','Yes')) d.1000$overall_second_line.factor <- factor(d.1000$overall_second_line, levels=c(0,1), labels=c('No','Yes')) second_line_summ. <- summary(composite_outcome.1000.factor ~ age_fhaart + male.factor + transmission_risk_cat.factor + pre_artaids_all + aids_at_fhaart.factor + baseline_cd4 + baseline_log_rna + site.factor + overall_death_indicator.factor + overall_first_failure.1000.factor + days_from_fhaart_to_vl_failure.1000 + overall_second_line.factor + days_from_fhaart_to_second_line_start, data=d.1000,subset=site != 'haiti',method='reverse',overall=TRUE) second_line_summ.$labels <- c('Age at first HAART','Sex', 'Possible route of infection', 'AIDS at first HAART (all)', 'AIDS at first HAART (clinical)', 'CD4 at start of first HAART', 'HIV1-RNA at start of first HAART (log-transformed)','Site', 'Died','Experienced virologic failure', 'Days from first HAART to failure', 'Major regimen change','Days from HAART start to 2nd line failure', 'Days from first HAART to major regimen change') # latex(second_line_summ., file='', exclude1=FALSE,long=TRUE,longtable=TRUE,caption='Descriptive statistics across the composite outcome (excluding Haiti).',digits=2, landscape=TRUE) d.1000 <- upData(d.1000, drop=c('composite_outcome.1000.factor', 'overall_first_failure.1000.factor', 'overall_second_line.factor')) # \begin{figure} # \caption{Incidence of virologic failure by site} cum.inc.failure.argentina <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.failure.brazil <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.failure.chile <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.failure.honduras <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.failure.mexico <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.failure.peru <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) cum.inc.failure.overall <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge.site <- range(c(cum.inc.failure.argentina[[1]]$est[which(cum.inc.failure.argentina[[1]]$time <= 7.0)],cum.inc.failure.brazil[[1]]$est[which(cum.inc.failure.brazil[[1]]$time <= 7.0)], cum.inc.failure.chile[[1]]$est[which(cum.inc.failure.chile[[1]]$time <= 7.0)],cum.inc.failure.honduras[[1]]$est[which(cum.inc.failure.honduras[[1]]$time <= 7.0)], cum.inc.failure.mexico[[1]]$est[which(cum.inc.failure.mexico[[1]]$time <= 7.0)],cum.inc.failure.peru[[1]]$est[which(cum.inc.failure.peru[[1]]$time <= 7.0)], cum.inc.failure.overall[[1]]$est[which(cum.inc.failure.overall[[1]]$time <= 7.0)])) # pdf('Incidence_of_virologic_failure.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.argentina[[1]]$est[which(cum.inc.failure.argentina[[1]]$time <= 6.0)] ~ cum.inc.failure.argentina[[1]]$time[which(cum.inc.failure.argentina[[1]]$time <= 6.0)], type='l',axes=FALSE,xlab='Years from cART initiation',ylab='Cumulative incidence',ylim=y.rge.site,col=2,xlim=c(0,6)) lines(cum.inc.failure.brazil[[1]]$est[which(cum.inc.failure.brazil[[1]]$time <= 6.0)] ~ cum.inc.failure.brazil[[1]]$time[which(cum.inc.failure.brazil[[1]]$time <= 6.0)],col=3) lines(cum.inc.failure.chile[[1]]$est[which(cum.inc.failure.chile[[1]]$time <= 6.0)] ~ cum.inc.failure.chile[[1]]$time[which(cum.inc.failure.chile[[1]]$time <= 6.0)],col=4) lines(cum.inc.failure.honduras[[1]]$est[which(cum.inc.failure.honduras[[1]]$time <= 6.0)] ~ cum.inc.failure.honduras[[1]]$time[which(cum.inc.failure.honduras[[1]]$time <= 6.0)],col=6) lines(cum.inc.failure.mexico[[1]]$est[which(cum.inc.failure.mexico[[1]]$time <= 6.0)] ~ cum.inc.failure.mexico[[1]]$time[which(cum.inc.failure.mexico[[1]]$time <= 6.0)],col=7) lines(cum.inc.failure.peru[[1]]$est[which(cum.inc.failure.peru[[1]]$time <= 6.0)] ~ cum.inc.failure.peru[[1]]$time[which(cum.inc.failure.peru[[1]]$time <= 6.0)],col=8) lines(cum.inc.failure.overall[[1]]$est[which(cum.inc.failure.overall[[1]]$time <= 6.0)] ~ cum.inc.failure.overall[[1]]$time[which(cum.inc.failure.overall[[1]]$time <= 6.0)],col='black',lwd=2) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru','Overall'),fill=c(2:3,4,6:8,1),bty='n',cex=0.75) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure by site, ignoring the single VL > 1000 criteria.} cum.inc.failure.argentina.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'argentina')], fstatus=d$censor_failure[which(d$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.failure.brazil.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'brazil')], fstatus=d$censor_failure[which(d$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.failure.chile.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'chile')], fstatus=d$censor_failure[which(d$site.factor == 'chile')], cencode='0',rho=0) cum.inc.failure.honduras.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'honduras')], fstatus=d$censor_failure[which(d$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.failure.mexico.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'mexico')], fstatus=d$censor_failure[which(d$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.failure.peru.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor == 'peru')], fstatus=d$censor_failure[which(d$site.factor == 'peru')], cencode='0',rho=0) cum.inc.failure.overall.no1000 <- cuminc(ftime=d$total_followup_failure[which(d.1000$site.factor != 'haiti')], fstatus=d$censor_failure[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge.site <- range(c(cum.inc.failure.argentina[[1]]$est[which(cum.inc.failure.argentina[[1]]$time <= 7.0)],cum.inc.failure.brazil[[1]]$est[which(cum.inc.failure.brazil[[1]]$time <= 7.0)], cum.inc.failure.chile[[1]]$est[which(cum.inc.failure.chile[[1]]$time <= 7.0)],cum.inc.failure.honduras[[1]]$est[which(cum.inc.failure.honduras[[1]]$time <= 7.0)], cum.inc.failure.mexico[[1]]$est[which(cum.inc.failure.mexico[[1]]$time <= 7.0)],cum.inc.failure.peru[[1]]$est[which(cum.inc.failure.peru[[1]]$time <= 7.0)], cum.inc.failure.overall[[1]]$est[which(cum.inc.failure.overall[[1]]$time <= 7.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.argentina.no1000[[1]]$est[which(cum.inc.failure.argentina.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.argentina.no1000[[1]]$time[which(cum.inc.failure.argentina.no1000[[1]]$time <= 6.0)], type='l',axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site,col=2) lines(cum.inc.failure.brazil.no1000[[1]]$est[which(cum.inc.failure.brazil.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.brazil.no1000[[1]]$time[which(cum.inc.failure.brazil.no1000[[1]]$time <= 6.0)],col=3) lines(cum.inc.failure.chile.no1000[[1]]$est[which(cum.inc.failure.chile.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.chile.no1000[[1]]$time[which(cum.inc.failure.chile.no1000[[1]]$time <= 6.0)],col=4) lines(cum.inc.failure.honduras.no1000[[1]]$est[which(cum.inc.failure.honduras.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.honduras.no1000[[1]]$time[which(cum.inc.failure.honduras.no1000[[1]]$time <= 6.0)],col=6) lines(cum.inc.failure.mexico.no1000[[1]]$est[which(cum.inc.failure.mexico.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.mexico.no1000[[1]]$time[which(cum.inc.failure.mexico.no1000[[1]]$time <= 6.0)],col=7) lines(cum.inc.failure.peru.no1000[[1]]$est[which(cum.inc.failure.peru.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.peru.no1000[[1]]$time[which(cum.inc.failure.peru.no1000[[1]]$time <= 6.0)],col=8) lines(cum.inc.failure.overall.no1000[[1]]$est[which(cum.inc.failure.overall.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.overall.no1000[[1]]$time[which(cum.inc.failure.overall.no1000[[1]]$time <= 6.0)],col='black',lwd=2) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru','Overall'),fill=c(2:3,4,6:8,1),bty='n',cex=0.75) # \end{figure} # \begin{figure} # \caption{Incidence of major regimen change by site} cum.inc.reg_chg.argentina <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.reg_chg.brazil <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.reg_chg.chile <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.reg_chg.haiti <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.reg_chg.honduras <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.reg_chg.mexico <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.reg_chg.peru <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) cum.inc.reg_chg.overall.m.haiti <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) cum.inc.reg_chg.overall <- cuminc(ftime=d.1000$total_followup_reg_change, fstatus=d.1000$censor_reg_change, cencode='0',rho=0) y.rge.site <- range(c(cum.inc.reg_chg.argentina[[1]]$est[which(cum.inc.reg_chg.argentina[[1]]$time <= 6.0)], cum.inc.reg_chg.brazil[[1]]$est[which(cum.inc.reg_chg.brazil[[1]]$time <= 6.0)], cum.inc.reg_chg.chile[[1]]$est[which(cum.inc.reg_chg.chile[[1]]$time <= 6.0)], cum.inc.reg_chg.haiti[[1]]$est[which(cum.inc.reg_chg.haiti[[1]]$time <= 6.0)], cum.inc.reg_chg.honduras[[1]]$est[which(cum.inc.reg_chg.honduras[[1]]$time <= 6.0)], cum.inc.reg_chg.mexico[[1]]$est[which(cum.inc.reg_chg.mexico[[1]]$time <= 6.0)], cum.inc.reg_chg.peru[[1]]$est[which(cum.inc.reg_chg.peru[[1]]$time <= 6.0)], cum.inc.reg_chg.overall.m.haiti[[1]]$est[which(cum.inc.reg_chg.overall.m.haiti[[1]]$time <= 6.0)], cum.inc.reg_chg.overall[[1]]$est[which(cum.inc.reg_chg.overall[[1]]$time <= 6.0)])) # pdf('Incidence_of_major_regimen_change.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.reg_chg.argentina[[1]]$est[which(cum.inc.reg_chg.argentina[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.argentina[[1]]$time[which(cum.inc.reg_chg.argentina[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from cART initiation',ylab='Cumulative incidence',ylim=y.rge.site,col=2,xlim=c(0,6)) lines(cum.inc.reg_chg.brazil[[1]]$est[which(cum.inc.reg_chg.brazil[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.brazil[[1]]$time[which(cum.inc.reg_chg.brazil[[1]]$time <= 6.0)],col=3) lines(cum.inc.reg_chg.chile[[1]]$est[which(cum.inc.reg_chg.chile[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.chile[[1]]$time[which(cum.inc.reg_chg.chile[[1]]$time <= 6.0)],col=4) lines(cum.inc.reg_chg.haiti[[1]]$est[which(cum.inc.reg_chg.haiti[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.haiti[[1]]$time[which(cum.inc.reg_chg.haiti[[1]]$time <= 6.0)],col=5) lines(cum.inc.reg_chg.honduras[[1]]$est[which(cum.inc.reg_chg.honduras[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.honduras[[1]]$time[which(cum.inc.reg_chg.honduras[[1]]$time <= 6.0)],col=6) lines(cum.inc.reg_chg.mexico[[1]]$est[which(cum.inc.reg_chg.mexico[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.mexico[[1]]$time[which(cum.inc.reg_chg.mexico[[1]]$time <= 6.0)],col=7) lines(cum.inc.reg_chg.peru[[1]]$est[which(cum.inc.reg_chg.peru[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.peru[[1]]$time[which(cum.inc.reg_chg.peru[[1]]$time <= 6.0)],lty=7,col=8) lines(cum.inc.reg_chg.overall.m.haiti[[1]]$est[which(cum.inc.reg_chg.overall.m.haiti[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.overall.m.haiti[[1]]$time[which(cum.inc.reg_chg.overall.m.haiti[[1]]$time <= 6.0)],lty=1,col=1, lwd=2) lines(cum.inc.reg_chg.overall[[1]]$est[which(cum.inc.reg_chg.overall[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.overall[[1]]$time[which(cum.inc.reg_chg.overall[[1]]$time <= 6.0)],lty=2,col=1,lwd=2) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Haiti','Honduras','Mexico','Peru','Overall (minus Haiti)','Overall (with Haiti)'),bty='n',cex=0.75,lty=c(rep(1,8),2),col=c(2,3,4,5,6,7,8,1,1)) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure or major regimen change by site} cum.inc.argentina <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.brazil <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.chile <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.haiti <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.honduras <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.mexico <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.peru <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor.1000[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) cum.inc.overall <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor.1000[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge.site <- range(c(cum.inc.argentina[[1]]$est[which(cum.inc.argentina[[1]]$time <= 6.0)], cum.inc.brazil[[1]]$est[which(cum.inc.brazil[[1]]$time <= 6.0)], cum.inc.chile[[1]]$est[which(cum.inc.chile[[1]]$time <= 6.0)], cum.inc.haiti[[1]]$est[which(cum.inc.haiti[[1]]$time <= 6.0)], cum.inc.honduras[[1]]$est[which(cum.inc.honduras[[1]]$time <= 6.0)], cum.inc.mexico[[1]]$est[which(cum.inc.mexico[[1]]$time <= 6.0)], cum.inc.peru[[1]]$est[which(cum.inc.peru[[1]]$time <= 6.0)], cum.inc.overall[[1]]$est[which(cum.inc.overall[[1]]$time <= 6.0)])) # pdf('Incidence_of_combined_outcome.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.argentina[[1]]$est[which(cum.inc.argentina[[1]]$time <= 6.0)] ~ cum.inc.argentina[[1]]$time[which(cum.inc.argentina[[1]]$time <= 6.0)], type='l',axes=FALSE, xlab='Years from cART initiation',ylab='Cumulative incidence',ylim=y.rge.site,col=2,xlim=c(0,6)) lines(cum.inc.brazil[[1]]$est[which(cum.inc.brazil[[1]]$time <= 6.0)] ~ cum.inc.brazil[[1]]$time[which(cum.inc.brazil[[1]]$time <= 6.0)],col=3) lines(cum.inc.chile[[1]]$est[which(cum.inc.chile[[1]]$time <= 6.0)] ~ cum.inc.chile[[1]]$time[which(cum.inc.chile[[1]]$time <= 6.0)],col=4) # lines(cum.inc.haiti[[1]]$est ~ cum.inc.haiti[[1]]$time,col=5) lines(cum.inc.honduras[[1]]$est[which(cum.inc.honduras[[1]]$time <= 6.0)] ~ cum.inc.honduras[[1]]$time[which(cum.inc.honduras[[1]]$time <= 6.0)],col=6) lines(cum.inc.mexico[[1]]$est[which(cum.inc.mexico[[1]]$time <= 6.0)] ~ cum.inc.mexico[[1]]$time[which(cum.inc.mexico[[1]]$time <= 6.0)],col=7) lines(cum.inc.peru[[1]]$est[which(cum.inc.peru[[1]]$time <= 6.0)] ~ cum.inc.peru[[1]]$time[which(cum.inc.peru[[1]]$time <= 6.0)],lty=7,col=8) lines(cum.inc.overall[[1]]$est[which(cum.inc.overall[[1]]$time <= 6.0)] ~ cum.inc.overall[[1]]$time[which(cum.inc.overall[[1]]$time <= 6.0)],col=1,lwd=2) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru','Overall'),bty='n',cex=0.75,fill=c(2,3,4,6,7,8,1)) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure, ignoring single VL > 1000 criteria, or major regimen change by site} cum.inc.argentina.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'argentina')], fstatus=d$censor[which(d$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.brazil.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'brazil')], fstatus=d$censor[which(d$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.chile.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'chile')], fstatus=d$censor[which(d$site.factor == 'chile')], cencode='0',rho=0) cum.inc.haiti.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'haiti')], fstatus=d$censor[which(d$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.honduras.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'honduras')], fstatus=d$censor[which(d$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.mexico.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'mexico')], fstatus=d$censor[which(d$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.peru.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor == 'peru')], fstatus=d$censor[which(d$site.factor == 'peru')], cencode='0',rho=0) cum.inc.overall.no1000 <- cuminc(ftime=d$total_followup[which(d.1000$site.factor != 'haiti')], fstatus=d$censor[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge.site.no1000 <- range(c(cum.inc.argentina.no1000[[1]]$est[which(cum.inc.argentina.no1000[[1]]$time <= 6.0)], cum.inc.brazil.no1000[[1]]$est[which(cum.inc.brazil.no1000[[1]]$time <= 6.0)], cum.inc.chile.no1000[[1]]$est[which(cum.inc.chile.no1000[[1]]$time <= 6.0)], cum.inc.haiti.no1000[[1]]$est[which(cum.inc.haiti.no1000[[1]]$time <= 6.0)], cum.inc.honduras.no1000[[1]]$est[which(cum.inc.honduras.no1000[[1]]$time <= 6.0)], cum.inc.mexico.no1000[[1]]$est[which(cum.inc.mexico.no1000[[1]]$time <= 6.0)], cum.inc.peru.no1000[[1]]$est[which(cum.inc.peru.no1000[[1]]$time <= 6.0)], cum.inc.overall.no1000[[1]]$est[which(cum.inc.overall.no1000[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.argentina.no1000[[1]]$est[which(cum.inc.argentina.no1000[[1]]$time <= 6.0)] ~ cum.inc.argentina.no1000[[1]]$time[which(cum.inc.argentina.no1000[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site,col=2) lines(cum.inc.brazil.no1000[[1]]$est[which(cum.inc.brazil.no1000[[1]]$time <= 6.0)] ~ cum.inc.brazil.no1000[[1]]$time[which(cum.inc.brazil.no1000[[1]]$time <= 6.0)],col=3) lines(cum.inc.chile.no1000[[1]]$est[which(cum.inc.chile.no1000[[1]]$time <= 6.0)] ~ cum.inc.chile.no1000[[1]]$time[which(cum.inc.chile.no1000[[1]]$time <= 6.0)],col=4) lines(cum.inc.honduras.no1000[[1]]$est[which(cum.inc.honduras.no1000[[1]]$time <= 6.0)] ~ cum.inc.honduras.no1000[[1]]$time[which(cum.inc.honduras.no1000[[1]]$time <= 6.0)],col=6) lines(cum.inc.mexico.no1000[[1]]$est[which(cum.inc.mexico.no1000[[1]]$time <= 6.0)] ~ cum.inc.mexico.no1000[[1]]$time[which(cum.inc.mexico.no1000[[1]]$time <= 6.0)],col=7) lines(cum.inc.peru.no1000[[1]]$est[which(cum.inc.peru.no1000[[1]]$time <= 6.0)] ~ cum.inc.peru.no1000[[1]]$time[which(cum.inc.peru.no1000[[1]]$time <= 6.0)],lty=7,col=8) lines(cum.inc.overall.no1000[[1]]$est[which(cum.inc.overall.no1000[[1]]$time <= 6.0)] ~ cum.inc.overall.no1000[[1]]$time[which(cum.inc.overall.no1000[[1]]$time <= 6.0)],col=1,lwd=2) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru','Overall'),bty='n',cex=0.75,fill=c(2,3,4,6,7,8,1)) # \end{figure} # \begin{figure} # \caption{Cumulative incidence of virologic failure, major regimen change, and the composite endpoint excluding Haiti for all outcomes.} cum.inc.failure <- cuminc(ftime=d.1000$total_followup_failure.1000[which(d.1000$site.factor != 'haiti')],fstatus=d.1000$censor_failure.1000[which(d.1000$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.reg_change <- cuminc(ftime=d.1000$total_followup_reg_change, fstatus=d.1000$censor_reg_change,cencode='0',rho=0) cum.inc.reg_change_no_haiti <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.reg_change_only_haiti <- cuminc(ftime=d.1000$total_followup_reg_change[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor_reg_change[which(d.1000$site.factor == 'haiti')],cencode='0',rho=0) cum.inc.composite_outcome <- cuminc(ftime=d.1000$total_followup.1000[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor.1000[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge <- range(c(cum.inc.failure[[1]]$est[which(cum.inc.failure[[1]]$time <= 6.0)], cum.inc.reg_change_no_haiti[[1]]$est[which(cum.inc.reg_change_no_haiti[[1]]$time <= 6.0)], cum.inc.composite_outcome[[1]]$est[which(cum.inc.composite_outcome[[1]]$time <= 6.0)])) # pdf('Cumulative_incidence_all_outcomes_excluding_Haiti.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure[[1]]$est[which(cum.inc.failure[[1]]$time <= 6.0)] ~ cum.inc.failure[[1]]$time[which(cum.inc.failure[[1]]$time <= 6.0)], type='l',axes=FALSE, xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge,lty=1) lines(cum.inc.reg_change_no_haiti[[1]]$est[which(cum.inc.reg_change_no_haiti[[1]]$time <= 6.0)] ~ cum.inc.reg_change_no_haiti[[1]]$time[which(cum.inc.reg_change_no_haiti[[1]]$time <= 6.0)], lty=2) lines(cum.inc.composite_outcome[[1]]$est[which(cum.inc.composite_outcome[[1]]$time <= 6.0)] ~ cum.inc.composite_outcome[[1]]$time[which(cum.inc.composite_outcome[[1]]$time <= 6.0)], lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Failure','Major regimen change','Composite outcome'),lty=c(1,2,3),cex=0.75) # dev.off() # \end{figure} # \begin{figure} # \caption{Cumulative incidence of virologic failure, ignoring the single VL > 1000 criteria, major regimen change, and the composite endpoint.} cum.inc.failure.no1000 <- cuminc(ftime=d$total_followup_failure[which(d$site.factor != 'haiti')],fstatus=d$censor_failure[which(d$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.reg_change.no1000 <- cuminc(ftime=d$total_followup_reg_change, fstatus=d$censor_reg_change,cencode='0',rho=0) cum.inc.reg_change.no1000_no_haiti <- cuminc(ftime=d$total_followup_reg_change[which(d$site.factor != 'haiti')], fstatus=d$censor_reg_change[which(d$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.composite_outcome.no1000 <- cuminc(ftime=d$total_followup[which(d$site.factor != 'haiti')], fstatus=d$censor[which(d$site.factor != 'haiti')], cencode='0',rho=0) y.rge.no1000 <- range(c(cum.inc.failure.no1000[[1]]$est[which(cum.inc.failure.no1000[[1]]$time <= 6.0)], cum.inc.reg_change.no1000_no_haiti[[1]]$est[which(cum.inc.reg_change.no1000_no_haiti[[1]]$time <= 6.0)], cum.inc.composite_outcome.no1000[[1]]$est[which(cum.inc.composite_outcome.no1000[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.no1000[[1]]$est[which(cum.inc.failure.no1000[[1]]$time <= 6.0)] ~ cum.inc.failure.no1000[[1]]$time[which(cum.inc.failure.no1000[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.no1000,lty=1) lines(cum.inc.reg_change.no1000_no_haiti[[1]]$est[which(cum.inc.reg_change.no1000_no_haiti[[1]]$time <= 6.0)] ~ cum.inc.reg_change.no1000_no_haiti[[1]]$time[which(cum.inc.reg_change.no1000_no_haiti[[1]]$time <= 6.0)], lty=2) lines(cum.inc.composite_outcome.no1000[[1]]$est[which(cum.inc.composite_outcome.no1000[[1]]$time <= 6.0)] ~ cum.inc.composite_outcome.no1000[[1]]$time[which(cum.inc.composite_outcome.no1000[[1]]$time <= 6.0)], lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Failure','Major regimen change','Composite outcome'),lty=c(1,2,3),cex=0.75) # \end{figure} # Failure junk <- subset(d.1000, site.factor != 'haiti') num_in_fu_failure_1 <- length(which(junk$total_followup_failure.1000 > 1.0)) tbl_1_failure <- table(junk$censor_failure.1000[which(junk$total_followup_failure.1000 <= 1.0)]) num_w_evt_failure_1 <- tbl_1_failure[2] num_w_cr_failure_1 <- tbl_1_failure[3] num_end_fu_failure_1 <- tbl_1_failure[1] for(i in 2:floor(max(junk$total_followup_failure.1000))){ assign(x=paste('num_in_fu_failure_',i,sep=''), value=length(which(junk$total_followup_failure.1000 > i))) tbl <- table(junk$censor_failure.1000[which(junk$total_followup_failure.1000 <= i & junk$total_followup_failure.1000 > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_failure_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_failure_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_failure_',i,sep=''), value=tmp) } # Major regimen change num_in_fu_reg_change_1 <- length(which(d.1000$total_followup_reg_change > 1.0)) tbl_1_reg_change <- table(d.1000$censor_reg_change[which(d.1000$total_followup_reg_change <= 1.0)]) num_w_evt_reg_change_1 <- tbl_1_reg_change[2] num_w_cr_reg_change_1 <- tbl_1_reg_change[3] num_end_fu_reg_change_1 <- tbl_1_reg_change[1] for(i in 2:floor(max(d.1000$total_followup_reg_change))){ assign(x=paste('num_in_fu_reg_change_',i,sep=''), value=length(which(d.1000$total_followup_reg_change > i))) tbl <- table(d.1000$censor_reg_change[which(d.1000$total_followup_reg_change <= i & d.1000$total_followup_reg_change > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_reg_change_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_reg_change_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_reg_change_',i,sep=''), value=tmp) } # Composite outcome # Switched from d.1000 to junk to eliminate Haiti subjects from this num_in_fu_composite_outcome_1 <- length(which(junk$total_followup.1000 > 1.0)) tbl_1_composite_outcome <- table(junk$censor.1000[which(junk$total_followup.1000 <= 1.0)]) num_w_evt_composite_outcome_1 <- tbl_1_composite_outcome[2] num_w_cr_composite_outcome_1 <- tbl_1_composite_outcome[3] num_end_fu_composite_outcome_1 <- tbl_1_composite_outcome[1] for(i in 2:floor(max(junk$total_followup.1000))){ assign(x=paste('num_in_fu_composite_outcome_',i,sep=''), value=length(which(junk$total_followup.1000 > i))) tbl <- table(junk$censor.1000[which(junk$total_followup.1000 <= i & junk$total_followup.1000 > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_composite_outcome_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_composite_outcome_',i,sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_composite_outcome_',i,sep=''), value=tmp) } # Finding cumulative incidence and 95% CI for each curve in Figure 4 at 1, 3, and 5 years after HAART initiation # Failure failure_inc_1 <- format(round(cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 1.0)][1],3),nsmall=3) failure_lci_1 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 1.0)][1] failure_uci_1 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 1.0)][1] failure_ci_1 <- paste('(',round(failure_lci_1,3),', ',round(failure_uci_1,3),')',sep='') failure_inc_3 <- format(round(cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 3.0)][1],3),nsmall=3) failure_lci_3 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 3.0)][1] failure_uci_3 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 3.0)][1] failure_ci_3 <- paste('(',round(failure_lci_3,3),', ',round(failure_uci_3,3),')',sep='') failure_inc_5 <- format(round(cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 5.0)][1],3),nsmall=3) failure_lci_5 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 5.0)][1] failure_uci_5 <- cum.inc.failure[[1]]$est[which(round(cum.inc.failure[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure[[1]]$var)[which(round(cum.inc.failure[[1]]$time,1) == 5.0)][1] failure_ci_5 <- paste('(',round(failure_lci_5,3),', ',round(failure_uci_5,3),')',sep='') # Failure, ignoring single VL > 1000 criteria failure_inc_1.no1000 <- format(round(cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 1.0)][1],4),nsmall=3) failure_lci_1.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var[which(round(cum.inc.failure.no1000[[1]]$time,1) == 1.0)][1]) failure_uci_1.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var [which(round(cum.inc.failure.no1000[[1]]$time,1) == 1.0)][1]) failure_ci_1.no1000 <- paste('(',round(failure_lci_1.no1000,3),', ',round(failure_uci_1.no1000,4),')',sep='') failure_inc_3.no1000 <- format(round(cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 3.0)][1],3),nsmall=3) failure_lci_3.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var [which(round(cum.inc.failure.no1000[[1]]$time,1) == 3.0)][1]) failure_uci_3.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var [which(round(cum.inc.failure.no1000[[1]]$time,1) == 3.0)][1]) failure_ci_3.no1000 <- paste('(',round(failure_lci_3.no1000,3),', ',round(failure_uci_3.no1000,3),')',sep='') failure_inc_5.no1000 <- format(round(cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 5.0)][1],3),nsmall=3) failure_lci_5.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var [which(round(cum.inc.failure.no1000[[1]]$time,1) == 5.0)][1]) failure_uci_5.no1000 <- cum.inc.failure.no1000[[1]]$est[which(round(cum.inc.failure.no1000[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000[[1]]$var [which(round(cum.inc.failure.no1000[[1]]$time,1) == 5.0)][1]) failure_ci_5.no1000 <- paste('(',round(failure_lci_5.no1000,3),', ',round(failure_uci_5.no1000,3),')',sep='') # Major regimen change # With all sites reg_change_inc_1 <- format(round(cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_1 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 1.0)][1] reg_change_uci_1 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 1.0)][1] reg_change_ci_1 <- paste('(',round(reg_change_lci_1,3),', ',round(reg_change_uci_1,3),')',sep='') reg_change_inc_3 <- format(round(cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_3 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 3.0)][1] reg_change_uci_3 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 3.0)][1] reg_change_ci_3 <- paste('(',round(reg_change_lci_3,3),', ',round(reg_change_uci_3,3),')',sep='') reg_change_inc_5 <- format(round(cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_5 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 5.0)][1] reg_change_uci_5 <- cum.inc.reg_change[[1]]$est[which(round(cum.inc.reg_change[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change[[1]]$var)[which(round(cum.inc.reg_change[[1]]$time,1) == 5.0)][1] reg_change_ci_5 <- paste('(',round(reg_change_lci_5,3),', ',round(reg_change_uci_5,3),')',sep='') # Without Haiti reg_change_inc_no_haiti_1 <- format(round(cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_no_haiti_1 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 1.0)][1] reg_change_uci_no_haiti_1 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 1.0)][1] reg_change_ci_no_haiti_1 <- paste('(',round(reg_change_lci_no_haiti_1,3),', ',round(reg_change_uci_no_haiti_1,3),')',sep='') reg_change_inc_no_haiti_3 <- format(round(cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_no_haiti_3 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 3.0)][1] reg_change_uci_no_haiti_3 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 3.0)][1] reg_change_ci_no_haiti_3 <- paste('(',round(reg_change_lci_no_haiti_3,3),', ',round(reg_change_uci_no_haiti_3,3),')',sep='') reg_change_inc_no_haiti_5 <- format(round(cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_no_haiti_5 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 5.0)][1] reg_change_uci_no_haiti_5 <- cum.inc.reg_change_no_haiti[[1]]$est[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti[[1]]$var)[which(round(cum.inc.reg_change_no_haiti[[1]]$time,1) == 5.0)][1] reg_change_ci_no_haiti_5 <- paste('(',round(reg_change_lci_no_haiti_5,3),', ',round(reg_change_uci_no_haiti_5,3),')',sep='') # With only Haiti reg_change_inc_only_haiti_1 <- format(round(cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_only_haiti_1 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 1.0)][1] reg_change_uci_only_haiti_1 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 1.0)][1] reg_change_ci_only_haiti_1 <- paste('(',round(reg_change_lci_only_haiti_1,3),', ',round(reg_change_uci_only_haiti_1,3),')',sep='') reg_change_inc_only_haiti_3 <- format(round(cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_only_haiti_3 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 3.0)][1] reg_change_uci_only_haiti_3 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 3.0)][1] reg_change_ci_only_haiti_3 <- paste('(',round(reg_change_lci_only_haiti_3,3),', ',round(reg_change_uci_only_haiti_3,3),')',sep='') reg_change_inc_only_haiti_5 <- format(round(cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_only_haiti_5 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 5.0)][1] reg_change_uci_only_haiti_5 <- cum.inc.reg_change_only_haiti[[1]]$est[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti[[1]]$var)[which(round(cum.inc.reg_change_only_haiti[[1]]$time,1) == 5.0)][1] reg_change_ci_only_haiti_5 <- paste('(',round(reg_change_lci_only_haiti_5,3),', ',round(reg_change_uci_only_haiti_5,3),')',sep='') # Composite outcome composite_outcome_inc_1 <- format(round(cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 1.0)][1],3),nsmall=3) composite_outcome_lci_1 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 1.0)][1] composite_outcome_uci_1 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 1.0)][1] composite_outcome_ci_1 <- paste('(',round(composite_outcome_lci_1,3),', ',round(composite_outcome_uci_1,3),')',sep='') composite_outcome_inc_3 <- format(round(cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 3.0)][1],3),nsmall=3) composite_outcome_lci_3 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 3.0)][1] composite_outcome_uci_3 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 3.0)][1] composite_outcome_ci_3 <- paste('(',round(composite_outcome_lci_3,3),', ',round(composite_outcome_uci_3,3),')',sep='') composite_outcome_inc_5 <- format(round(cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 5.0)][1],3),nsmall=3) composite_outcome_lci_5 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 5.0)][1] composite_outcome_uci_5 <- cum.inc.composite_outcome[[1]]$est[which(round(cum.inc.composite_outcome[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome[[1]]$var)[which(round(cum.inc.composite_outcome[[1]]$time,1) == 5.0)][1] composite_outcome_ci_5 <- paste('(',round(composite_outcome_lci_5,3),', ',round(composite_outcome_uci_5,3),')',sep='') # Composite outcome, ignoring single VL > 1000 criteria composite_outcome_inc_1.no1000 <- format(round(cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 1.0)][1],3),nsmall=3) composite_outcome_lci_1.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 1.0)][1] composite_outcome_uci_1.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 1.0)][1] composite_outcome_ci_1.no1000 <- paste('(',round(composite_outcome_lci_1.no1000,3),', ',round(composite_outcome_uci_1.no1000,3),')',sep='') composite_outcome_inc_3.no1000 <- format(round(cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 3.0)][1],3),nsmall=3) composite_outcome_lci_3.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 3.0)][1] composite_outcome_uci_3.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 3.0)][1] composite_outcome_ci_3.no1000 <- paste('(',round(composite_outcome_lci_3.no1000,3),', ',round(composite_outcome_uci_3.no1000,3),')',sep='') composite_outcome_inc_5.no1000 <- format(round(cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 5.0)][1],3),nsmall=3) composite_outcome_lci_5.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 5.0)][1] composite_outcome_uci_5.no1000 <- cum.inc.composite_outcome.no1000[[1]]$est[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000[[1]]$var)[which(round(cum.inc.composite_outcome.no1000[[1]]$time,1) == 5.0)][1] composite_outcome_ci_5.no1000 <- paste('(',round(composite_outcome_lci_5.no1000,3),', ',round(composite_outcome_uci_5.no1000,3),')',sep='') fu_df <- data.frame('Year'=rep(seq(1,12),times=3), 'Had event'=c(num_w_evt_failure_1,num_w_evt_failure_2,num_w_evt_failure_3,num_w_evt_failure_4, num_w_evt_failure_5,num_w_evt_failure_6,num_w_evt_failure_7,num_w_evt_failure_8, num_w_evt_failure_9,num_w_evt_failure_10,num_w_evt_failure_11,num_w_evt_failure_12, num_w_evt_reg_change_1,num_w_evt_reg_change_2,num_w_evt_reg_change_3,num_w_evt_reg_change_4, num_w_evt_reg_change_5,num_w_evt_reg_change_6,num_w_evt_reg_change_7,num_w_evt_reg_change_8, num_w_evt_reg_change_9,num_w_evt_reg_change_10,num_w_evt_reg_change_11,num_w_evt_reg_change_12, num_w_evt_composite_outcome_1,num_w_evt_composite_outcome_2,num_w_evt_composite_outcome_3,num_w_evt_composite_outcome_4, num_w_evt_composite_outcome_5,num_w_evt_composite_outcome_6,num_w_evt_composite_outcome_7,num_w_evt_composite_outcome_8, num_w_evt_composite_outcome_9,num_w_evt_composite_outcome_10,num_w_evt_composite_outcome_11,num_w_evt_composite_outcome_12), 'Had competing risk'=c(num_w_cr_failure_1,num_w_cr_failure_2,num_w_cr_failure_3,num_w_cr_failure_4, num_w_cr_failure_5,num_w_cr_failure_6,num_w_cr_failure_7,num_w_cr_failure_8, num_w_cr_failure_9,num_w_cr_failure_10,num_w_cr_failure_11,num_w_cr_failure_12, num_w_cr_reg_change_1,num_w_cr_reg_change_2,num_w_cr_reg_change_3,num_w_cr_reg_change_4, num_w_cr_reg_change_5,num_w_cr_reg_change_6,num_w_cr_reg_change_7,num_w_cr_reg_change_8, num_w_cr_reg_change_9,num_w_cr_reg_change_10,num_w_cr_reg_change_11,num_w_cr_reg_change_12, num_w_cr_composite_outcome_1,num_w_cr_composite_outcome_2,num_w_cr_composite_outcome_3,num_w_cr_composite_outcome_4, num_w_cr_composite_outcome_5,num_w_cr_composite_outcome_6,num_w_cr_composite_outcome_7,num_w_cr_composite_outcome_8, num_w_cr_composite_outcome_9,num_w_cr_composite_outcome_10,num_w_cr_composite_outcome_11,num_w_cr_composite_outcome_12), 'Censored'=c(num_end_fu_failure_1,num_end_fu_failure_2,num_end_fu_failure_3,num_end_fu_failure_4, num_end_fu_failure_5,num_end_fu_failure_6,num_end_fu_failure_7,num_end_fu_failure_8, num_end_fu_failure_9,num_end_fu_failure_10,num_end_fu_failure_11,num_end_fu_failure_12, num_end_fu_reg_change_1,num_end_fu_reg_change_2,num_end_fu_reg_change_3,num_end_fu_reg_change_4, num_end_fu_reg_change_5,num_end_fu_reg_change_6,num_end_fu_reg_change_7,num_end_fu_reg_change_8, num_end_fu_reg_change_9,num_end_fu_reg_change_10,num_end_fu_reg_change_11,num_end_fu_reg_change_12, num_end_fu_composite_outcome_1,num_end_fu_composite_outcome_2,num_end_fu_composite_outcome_3,num_end_fu_composite_outcome_4, num_end_fu_composite_outcome_5,num_end_fu_composite_outcome_6,num_end_fu_composite_outcome_7,num_end_fu_composite_outcome_8, num_end_fu_composite_outcome_9,num_end_fu_composite_outcome_10,num_end_fu_composite_outcome_11,num_end_fu_composite_outcome_12), 'Still in follow-up'=c(num_in_fu_failure_1,num_in_fu_failure_2,num_in_fu_failure_3,num_in_fu_failure_4, num_in_fu_failure_5,num_in_fu_failure_6,num_in_fu_failure_7,num_in_fu_failure_8, num_in_fu_failure_9,num_in_fu_failure_10,num_in_fu_failure_11,num_in_fu_failure_12, num_in_fu_reg_change_1,num_in_fu_reg_change_2,num_in_fu_reg_change_3,num_in_fu_reg_change_4, num_in_fu_reg_change_5,num_in_fu_reg_change_6,num_in_fu_reg_change_7,num_in_fu_reg_change_8, num_in_fu_reg_change_9,num_in_fu_reg_change_10,num_in_fu_reg_change_11,num_in_fu_reg_change_12, num_in_fu_composite_outcome_1,num_in_fu_composite_outcome_2,num_in_fu_composite_outcome_3,num_in_fu_composite_outcome_4, num_in_fu_composite_outcome_5,num_in_fu_composite_outcome_6,num_in_fu_composite_outcome_7,num_in_fu_composite_outcome_8, num_in_fu_composite_outcome_9,num_in_fu_composite_outcome_10,num_in_fu_composite_outcome_11,num_in_fu_composite_outcome_12)) fu_df <- as.matrix(fu_df) colnames(fu_df) <- c('Year','Had event','Had competing risk','Censored','Still in follow-up') # latex(fu_df, file='', rgroup=c('Virologic failure','Major regimen change','Composite outcome'),n.rgroup=c(12,12,12),rowname=c(''),rowlabel='',where='!h',longtable=TRUE,caption='Number of subjects having the event, the competing risk, or were censored by year of follow-up. Virologic failure and the composite outcome do not include subjects from Haiti.') sub_d.1000 <- subset(d.1000, select=c(age_fhaart, male.factor, aids_at_fhaart,aids_at_fhaart.factor, fhaart_regimen_class.factor,nrti_type.factor, baseline_cd4, baseline_log_rna, transmission_risk_cat.factor,transmission_risk_larger_cat.factor, total_followup.1000,total_followup_failure.1000,total_followup_reg_change, composite_outcome.1000,site,site.factor,year_fhaart, days_from_fhaart_to_vl_failure.1000,overall_first_failure.1000, days_from_fhaart_to_vl_failure,overall_first_failure,max_date,date_fhaart, days_from_fhaart_to_second_line_start,overall_second_line, year_gap_date)) # Changed the following two lines to exclude subjects from Haiti followup.median.1000 <- format(round(median(d.1000$total_followup.1000[which(d.1000$site.factor != 'haiti')]),2),nsmall=2) followup.iqr.1000 <- format(round(quantile(d.1000$total_followup.1000[which(d.1000$site.factor != 'haiti')],probs=c(0.25,0.75)),2),nsmall=2) t2evt_no_haiti <- time_to_evt_function(df=sub_d.1000,include_haiti=FALSE,which_outcome='composite',only_haiti=FALSE,reduced=FALSE,censor=FALSE) t2evt_w_haiti <- time_to_evt_function(df=sub_d.1000,include_haiti=TRUE,which_outcome='composite',only_haiti=FALSE,reduced=TRUE,censor=FALSE) # Need outcome for failure t2evt_no_haiti_failure <- time_to_evt_function(df=sub_d.1000,include_haiti=FALSE,which_outcome='failure',only_haiti=FALSE,reduced=FALSE,censor=FALSE) # Need outcome for major regimen change t2evt_haiti_regimen_change <- time_to_evt_function(df=sub_d.1000, include_haiti=TRUE,which_outcome='regimen',only_haiti=FALSE,reduced=TRUE,censor=FALSE) t2evt_no_haiti_regimen_change <- time_to_evt_function(df=sub_d.1000, include_haiti=FALSE,which_outcome='regimen',only_haiti=FALSE,reduced=FALSE,censor=FALSE) t2evt_only_haiti_regimen_change <- time_to_evt_function(df=sub_d.1000, include_haiti=TRUE,which_outcome='regimen',only_haiti=TRUE,reduced=TRUE,censor=FALSE) t2evt_no_haiti_reduced_regimen_change <- time_to_evt_function(df=sub_d.1000, include_haiti=FALSE,which_outcome='regimen',only_haiti=FALSE,reduced=TRUE,censor=FALSE) sub_d <- subset(d, select=c(age_fhaart, male.factor, aids_at_fhaart,aids_at_fhaart.factor, fhaart_regimen_class.factor,fhaart_regimen_class.factor,nrti_type.factor, baseline_cd4, baseline_log_rna, transmission_risk_cat.factor,transmission_risk_larger_cat.factor, total_followup,total_followup_failure,total_followup_reg_change, composite_outcome,site,site.factor,year_fhaart, days_from_fhaart_to_vl_failure,overall_first_failure, days_from_fhaart_to_vl_failure,overall_first_failure,max_date,date_fhaart, days_from_fhaart_to_second_line_start,overall_second_line)) followup.median <- format(round(median(d$total_followup),2),nsmall=2) followup.iqr <- format(round(quantile(d$total_followup,probs=c(0.25,0.75)),2),nsmall=2) t2evt_no_haiti.no1000 <- time_to_evt_function(df=sub_d,include_haiti=FALSE,which_outcome='composite.no1000',only_haiti=FALSE,reduced=FALSE,censor=FALSE) # Need outcome for failure t2evt_no_haiti_failure.no1000 <- time_to_evt_function(df=sub_d,include_haiti=FALSE,which_outcome='failure.no1000',only_haiti=FALSE,reduced=FALSE,censor=FALSE) haiti.df <- t2evt_w_haiti[['primary.df']][-(28:36),] # latex(haiti.df, file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and adjusting for age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,where='!h',size='small',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) haiti_imputed.df <- t2evt_w_haiti[['primary.df.impute']][-(28:36),] # latex(haiti_imputed.df, file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and adjusting for age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_haiti_regimen_change[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome, stratified by site, and including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_haiti_regimen_change[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome, stratified by site, and including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_regimen_change[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome, stratified by site, and excluding Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, log viral load, and probable route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_regimen_change[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome, stratified by site, and excluding Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, log viral load, and probable route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_only_haiti_regimen_change[['primary.df']],file='',caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome and including only data from Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'), n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_only_haiti_regimen_change[['primary.df.impute']],file='',caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome and including only data from Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'), n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_reduced_regimen_change[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome and excluding Haiti from the analysis. This analysis has a reduced list of covariates from the typical analysis excluding Haiti: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_reduced_regimen_change[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome and excluding Haiti from the analysis. This analysis has a reduced list of covariates from the typical analysis excluding Haiti: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti.no1000[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti.no1000[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure.no1000[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure.no1000[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # Creating Table 2 for manuscript table2 <- data.frame('Covariate'=t2evt_no_haiti[['primary.df.impute']][,1], t2evt_no_haiti_failure[['primary.df.impute']][,5:7], t2evt_no_haiti_regimen_change[['primary.df.impute']][,5:7], t2evt_no_haiti[['primary.df.impute']][,5:7]) table2 <- as.matrix(table2) colnames(table2) <- c('Covariate',rep(c('HR','95\\% CI','P'),times=3)) # latex(table2, file='', landscape=TRUE,cgroup=c('','Failure','Major regimen change','Composite outcome'),n.cgroup=c(1,3,3,3),caption='Cox regression using imputed data and stratified by site and excluding Haiti from the analysis.',exclude1=FALSE,long=TRUE,longtable=TRUE) # Creating Table 3 for manuscript table3 <- data.frame('Covariate'=t2evt_only_haiti_regimen_change[['primary.df.impute']][,1], t2evt_only_haiti_regimen_change[['primary.df.impute']][,5:7], t2evt_no_haiti_reduced_regimen_change[['primary.df.impute']][,5:7]) table3 <- as.matrix(table3) colnames(table3) <- c('Covariate',rep(c('HR','95\\% CI','P'),times=2)) # latex(table3, file='', landscape=TRUE,cgroup=c('','Only Haiti','Excluding Haiti'),n.cgroup=c(1,3,3),caption='Cox regression using imputed data and stratified by site running the reduced model with only Haiti and excluding Haiti.',longtable=TRUE) # Median frequency of viral load measurements per site -- # of viral load measurements dividied by the follow-up time after HAART initiation per person. Then compute median by site. # First need to add total_followup to rna rna$total_followup.1000 <- d.1000[match(rna$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] rna$site.factor <- art[match(rna$unique_id,art$unique_id,nomatch=NA),'site.factor'] tmp <- subset(rna, !is.na(total_followup.1000) & site.factor != 'haiti' & !is.na(rna_v)) tmp$unique_id.factor <- factor(tmp$unique_id, levels=unique(d.1000$unique_id)) num_vl_measurements <- tapply(tmp$rna_v, INDEX=tmp$unique_id, FUN=function(y){length(y[which(!is.na(y))])}) num_vl_measurements.df <- data.frame('unique_id'=names(num_vl_measurements), 'num_vl_measurements'=num_vl_measurements) tmp$num_vl_per_person <- num_vl_measurements.df[match(tmp$unique_id,num_vl_measurements.df$unique_id,nomatch=NA),'num_vl_measurements'] tmp$vl_frequency <- tmp$num_vl_per_person/tmp$total_followup.1000 all_visits$vl_frequency <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'vl_frequency'] all_visits$num_vl_per_person <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'num_vl_per_person'] vl_frequency_by_site <- tapply(all_visits$vl_frequency[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) vl_frequency_overall <- quantile(all_visits$vl_frequency[!duplicated(all_visits$unique_id)], probs=c(0.25,0.50,0.75),na.rm=TRUE) vl_per_person_by_site <- tapply(all_visits$num_vl_per_person[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) # Median frequency of CD4 measurements per site -- # of CD4 measurements dividied by the follow-up time after HAART initiation per person. Then compute median by site. # First need to add total_followup to cd4 cd4$total_followup.1000 <- d.1000[match(cd4$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] tmp <- subset(cd4, !is.na(total_followup.1000) & !is.na(cd4_v)) tmp$unique_id.factor <- factor(tmp$unique_id, levels=unique(d.1000$unique_id)) num_cd4_measurements <- tapply(tmp$cd4_v, INDEX=tmp$unique_id, FUN=function(y){length(y[which(!is.na(y))])}) num_cd4_measurements.df <- data.frame('unique_id'=names(num_cd4_measurements), 'num_cd4_measurements'=num_cd4_measurements) tmp$num_cd4_per_person <- num_cd4_measurements.df[match(tmp$unique_id,num_cd4_measurements.df$unique_id,nomatch=NA),'num_cd4_measurements'] tmp$cd4_frequency <- tmp$num_cd4_per_person/tmp$total_followup.1000 all_visits$cd4_frequency <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'cd4_frequency'] cd4_frequency_by_site <- tapply(all_visits$cd4_frequency[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) cd4_frequency_overall <- quantile(all_visits$cd4_frequency[!duplicated(all_visits$unique_id)], probs=c(0.25,0.50,0.75),na.rm=TRUE) all_visits$num_cd4_per_person <- tmp[match(all_visits$unique_id,cd4$unique_id,nomatch=NA),'num_cd4_per_person'] cd4_per_person_by_site <- tapply(all_visits$num_cd4_per_person[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) # For median follow-up time (IQR), calculating for the entire all_visits cohort (6598) and then for those in d.1000 (6510 -- excludes the 88 with 0 follow-up) all_visits$total_followup.1000 <- d.1000[match(all_visits$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] all_visits$total_followup.1000 <- ifelse(is.na(all_visits$total_followup.1000),0,all_visits$total_followup.1000) entire_fu <- describe(all_visits$total_followup.1000[!duplicated(all_visits$unique_id)]) d.1000_fu <- describe(d.1000$total_followup.1000[!duplicated(d.1000$unique_id)]) # Creating df for CD4 and VL frequencies by site df_site <- data.frame('Site'=names(cd4_frequency_by_site), 'By year'=round(cd4_frequency_by_site,2), 'Overall'=round(cd4_per_person_by_site,2), 'By year'=round(vl_frequency_by_site,2), 'Overall'=round(vl_per_person_by_site,2)) df_site <- as.matrix(df_site) colnames(df_site) <- c('Site','By year','Overall','By year','Overall') df_site[,4] <- ifelse(is.na(df_site[,4]),'',df_site[,4]) df_site[,5] <- ifelse(is.na(df_site[,5]),'',df_site[,5]) overall_cd4 <- format(round(cd4_frequency_overall[2],2),nsmall=2) overall_vl <- format(round(vl_frequency_overall[2],2),nsmall=2) # latex(df_site, rowname=NULL, file='', caption='Frequency of CD4 and viral load measurements by site.',where='!h',cgroup=c('','CD4','HIV-1 RNA'),n.cgroup=c(1,2,2), insert.bottom=paste('Median rate of CD4 measurement frequency=',overall_cd4,'; Median rate of HIV-1 RNA measurement frequency=',overall_vl,'.',sep='')) # Bryan needs the proportion of those with completely missing rna values # Some are not in the rna table: not.in.rna <- basic$unique_id[which(basic$unique_id %nin% rna$unique_id)] jk <- tapply(rna$rna_v, INDEX=rna$unique_id, FUN=function(y){as.numeric(all(is.na(y)))}) no.rna <- names(jk)[which(jk == 1)] missing.rna <- unique(c(not.in.rna, no.rna)) missing.rna.df <- subset(basic, unique_id %in% missing.rna) # Also need to look at whether those who failed a first line but didn't start a second line had any suppression information in ART after failing first line ids_w_failure_but_no_2nd_line <- unique(all_visits$unique_id[which(all_visits$overall_first_failure.1000 == 1 & all_visits$overall_second_line == 0)]) tr <- subset(art, unique_id %in% ids_w_failure_but_no_2nd_line & art_sd >= date_first_failure.1000) tr2 <- subset(art, unique_id %in% ids_w_failure_but_no_2nd_line) num_ids_w_no_regimens_past_first_failure <- length(which(unique(tr2$unique_id) %nin% unique(tr$unique_id))) ids_w_regimen_past_first_failure <- unique(tr$unique_id) tr[which(tr$unique_id == 'peru.465'),c('unique_id','art_id','art_sd','date_fhaart')] rna$date_first_failure.1000 <- all_visits[match(rna$unique_id, all_visits$unique_id,nomatch=NA),'date_first_failure.1000'] sub_rna <- subset(rna, rna_d >= date_first_failure.1000 & unique_id %in% tr$unique_id) sub_rna <- upData(sub_rna, drop='date_fhaart') sub_rna <- upData(sub_rna, rename=c('date_first_failure.1000'='date_fhaart')) testing_suppression <- first_failure_virologic_failure_function(switch='all',sub_rna,detectable_level=400) sub_rna <- subset(rna, unique_id %in% ids_w_failure_but_no_2nd_line) sub_rna <- subset(sub_rna, rna_d > date_first_failure.1000) testing_suppression2 <- first_failure_virologic_failure_function(switch='all',sub_rna, detectable_level=400) # # \subsection{Censoring VF that occurred after any gaps of at least a year in viral load measurements.} tbl_by_tbl_new <- with(tmp_2nd_line, table(overall_first_failure.1000_new.factor,overall_second_line_new.factor)) tmp_2nd_line_no_haiti$overall_death_indicator_new.factor <- all_visits[match(tmp_2nd_line_no_haiti$unique_id, all_visits$unique_id, nomatch=NA),'overall_death_indicator_new.factor'] second_line_summ._new <- summary(composite_outcome.1000_new.factor ~ age_fhaart + male.factor + transmission_risk_cat.factor + aids_at_fhaart.factor + baseline_cd4 + baseline_log_rna + site.factor + overall_death_indicator_new.factor + overall_first_failure.1000_new.factor + days_from_fhaart_to_vl_failure.1000_new + overall_second_line_new.factor + days_from_fhaart_to_second_line_start_new, data=tmp_2nd_line_no_haiti,method='reverse',overall=TRUE) second_line_summ._new$labels <- c('Age at first HAART','Sex', 'Possible route of infection', 'AIDS at first HAART', 'CD4 at start of first HAART', 'HIV1-RNA at start of first HAART (log-transformed)','Site', 'Died (accounting for censoring)', 'Experienced virologic failure (accounting for censoring)', 'Days from first HAART to failure (accounting for censoring)', 'Major regimen change (accounting for censoring)', 'Days from HAART start to 2nd line failure (accounting for censoring)', 'Days from first HAART to major regimen change (accounting for censoring)') # latex(second_line_summ._new, file='', exclude1=FALSE,long=TRUE,longtable=TRUE,caption='Descriptive statistics across the composite outcome (excluding Haiti) and accounting for censoring.',digits=2, landscape=TRUE) # \subsubsection{Time to event analysis (accounting for censoring)} #------------------------------------------------------------------------------------------------# # Time-to-event analysis # #------------------------------------------------------------------------------------------------# # \begin{figure} # \caption{Incidence of virologic failure by site, accounting for censoring.} cum.inc.failure.argentina_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.failure.brazil_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.failure.chile_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.failure.honduras_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.failure.mexico_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.failure.peru_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) y.rge.site_new <- range(c(cum.inc.failure.argentina_new[[1]]$est[which(cum.inc.failure.argentina_new[[1]]$time <= 6.0)], cum.inc.failure.brazil_new[[1]]$est[which(cum.inc.failure.brazil_new[[1]]$time <= 6.0)], cum.inc.failure.chile_new[[1]]$est[which(cum.inc.failure.chile_new[[1]]$time <= 6.0)], cum.inc.failure.honduras_new[[1]]$est[which(cum.inc.failure.honduras_new[[1]]$time <= 6.0)], cum.inc.failure.mexico_new[[1]]$est[which(cum.inc.failure.mexico_new[[1]]$time <= 6.0)], cum.inc.failure.peru_new[[1]]$est[which(cum.inc.failure.peru_new[[1]]$time <= 6.0)])) # pdf('Figure1_June132014.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.argentina_new[[1]]$est[which(cum.inc.failure.argentina_new[[1]]$time <= 6.0)] ~ cum.inc.failure.argentina_new[[1]]$time[which(cum.inc.failure.argentina_new[[1]]$time <= 6.0)], type='l',axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site_new,col=2,xlim=c(0,6)) lines(cum.inc.failure.brazil_new[[1]]$est[which(cum.inc.failure.brazil_new[[1]]$time <= 6.0)] ~ cum.inc.failure.brazil_new[[1]]$time[which(cum.inc.failure.brazil_new[[1]]$time <= 6.0)],col=3) lines(cum.inc.failure.chile_new[[1]]$est[which(cum.inc.failure.chile_new[[1]]$time <= 6.0)] ~ cum.inc.failure.chile_new[[1]]$time[which(cum.inc.failure.chile_new[[1]]$time <= 6.0)],col=4) lines(cum.inc.failure.honduras_new[[1]]$est[which(cum.inc.failure.honduras_new[[1]]$time <= 6.0)] ~ cum.inc.failure.honduras_new[[1]]$time[which(cum.inc.failure.honduras_new[[1]]$time <= 6.0)],col=6) lines(cum.inc.failure.mexico_new[[1]]$est[which(cum.inc.failure.mexico_new[[1]]$time <= 6.0)] ~ cum.inc.failure.mexico_new[[1]]$time[which(cum.inc.failure.mexico_new[[1]]$time <= 6.0)],col=7) lines(cum.inc.failure.peru_new[[1]]$est[which(cum.inc.failure.peru_new[[1]]$time <= 6.0)] ~ cum.inc.failure.peru_new[[1]]$time[which(cum.inc.failure.peru_new[[1]]$time <= 6.0)],col=8) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru'),fill=c(2:3,4,6:8),bty='n',cex=0.75) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure by site, ignoring the single VL > 1000 criteria, accounting for censoring.} cum.inc.failure.argentina.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'argentina')], fstatus=d$censor_failure_new[which(d$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.failure.brazil.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'brazil')], fstatus=d$censor_failure_new[which(d$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.failure.chile.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'chile')], fstatus=d$censor_failure_new[which(d$site.factor == 'chile')], cencode='0',rho=0) cum.inc.failure.honduras.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'honduras')], fstatus=d$censor_failure_new[which(d$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.failure.mexico.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'mexico')], fstatus=d$censor_failure_new[which(d$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.failure.peru.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor == 'peru')], fstatus=d$censor_failure_new[which(d$site.factor == 'peru')], cencode='0',rho=0) y.rge.site.no1000_new <- range(c(cum.inc.failure.argentina.no1000_new[[1]]$est[which(cum.inc.failure.argentina.no1000_new[[1]]$time <= 6.0)], cum.inc.failure.brazil.no1000_new[[1]]$est[which(cum.inc.failure.brazil.no1000_new[[1]]$time <= 6.0)], cum.inc.failure.chile.no1000_new[[1]]$est[which(cum.inc.failure.chile.no1000_new[[1]]$time <= 6.0)], cum.inc.failure.honduras.no1000_new[[1]]$est[which(cum.inc.failure.honduras.no1000_new[[1]]$time <= 6.0)], cum.inc.failure.mexico.no1000_new[[1]]$est[which(cum.inc.failure.mexico.no1000_new[[1]]$time <= 6.0)], cum.inc.failure.peru.no1000_new[[1]]$est[which(cum.inc.failure.peru.no1000_new[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.argentina.no1000_new[[1]]$est[which(cum.inc.failure.argentina.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.argentina.no1000_new[[1]]$time[which(cum.inc.failure.argentina.no1000_new[[1]]$time <= 6.0)], type='l',axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site_new,col=2) lines(cum.inc.failure.brazil.no1000_new[[1]]$est[which(cum.inc.failure.brazil.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.brazil.no1000_new[[1]]$time[which(cum.inc.failure.brazil.no1000_new[[1]]$time <= 6.0)],col=3) lines(cum.inc.failure.chile.no1000_new[[1]]$est[which(cum.inc.failure.chile.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.chile.no1000_new[[1]]$time[which(cum.inc.failure.chile.no1000_new[[1]]$time <= 6.0)],col=4) lines(cum.inc.failure.honduras.no1000_new[[1]]$est[which(cum.inc.failure.honduras.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.honduras.no1000_new[[1]]$time[which(cum.inc.failure.honduras.no1000_new[[1]]$time <= 6.0)],col=6) lines(cum.inc.failure.mexico.no1000_new[[1]]$est[which(cum.inc.failure.mexico.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.mexico.no1000_new[[1]]$time[which(cum.inc.failure.mexico.no1000_new[[1]]$time <= 6.0)],col=7) lines(cum.inc.failure.peru.no1000_new[[1]]$est[which(cum.inc.failure.peru.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.peru.no1000_new[[1]]$time[which(cum.inc.failure.peru.no1000_new[[1]]$time <= 6.0)],col=8) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru'),fill=c(2:3,4,6:8),bty='n',cex=0.75) # \end{figure} # \begin{figure} # \caption{Incidence of major regimen change by site, accounting for censoring.} cum.inc.reg_chg.argentina_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.reg_chg.brazil_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.reg_chg.chile_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.reg_chg.haiti_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.reg_chg.honduras_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.reg_chg.mexico_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.reg_chg.peru_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) y.rge.site_new <- range(c(cum.inc.reg_chg.argentina_new[[1]]$est[which(cum.inc.reg_chg.argentina_new[[1]]$time <= 6.0)], cum.inc.reg_chg.brazil_new[[1]]$est[which(cum.inc.reg_chg.brazil_new[[1]]$time <= 6.0)], cum.inc.reg_chg.chile_new[[1]]$est[which(cum.inc.reg_chg.chile_new[[1]]$time <= 6.0)], cum.inc.reg_chg.haiti_new[[1]]$est[which(cum.inc.reg_chg.haiti_new[[1]]$time <= 6.0)], cum.inc.reg_chg.honduras_new[[1]]$est[which(cum.inc.reg_chg.honduras_new[[1]]$time <= 6.0)], cum.inc.reg_chg.mexico_new[[1]]$est[which(cum.inc.reg_chg.mexico_new[[1]]$time <= 6.0)], cum.inc.reg_chg.peru_new[[1]]$est[which(cum.inc.reg_chg.peru_new[[1]]$time <= 6.0)])) # pdf('Figure3_June132014.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.reg_chg.argentina_new[[1]]$est[which(cum.inc.reg_chg.argentina_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.argentina_new[[1]]$time[which(cum.inc.reg_chg.argentina_new[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site_new,col=2,xlim=c(0,6)) lines(cum.inc.reg_chg.brazil_new[[1]]$est[which(cum.inc.reg_chg.brazil_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.brazil_new[[1]]$time[which(cum.inc.reg_chg.brazil_new[[1]]$time <= 6.0)],col=3) lines(cum.inc.reg_chg.chile_new[[1]]$est[which(cum.inc.reg_chg.chile_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.chile_new[[1]]$time[which(cum.inc.reg_chg.chile_new[[1]]$time <= 6.0)],col=4) lines(cum.inc.reg_chg.haiti_new[[1]]$est[which(cum.inc.reg_chg.haiti_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.haiti_new[[1]]$time[which(cum.inc.reg_chg.haiti_new[[1]]$time <= 6.0)],col=5) lines(cum.inc.reg_chg.honduras_new[[1]]$est[which(cum.inc.reg_chg.honduras_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.honduras_new[[1]]$time[which(cum.inc.reg_chg.honduras_new[[1]]$time <= 6.0)],col=6) lines(cum.inc.reg_chg.mexico_new[[1]]$est[which(cum.inc.reg_chg.mexico_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.mexico_new[[1]]$time[which(cum.inc.reg_chg.mexico_new[[1]]$time <= 6.0)],col=7) lines(cum.inc.reg_chg.peru_new[[1]]$est[which(cum.inc.reg_chg.peru_new[[1]]$time <= 6.0)] ~ cum.inc.reg_chg.peru_new[[1]]$time[which(cum.inc.reg_chg.peru_new[[1]]$time <= 6.0)],lty=7,col=8) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Haiti','Honduras','Mexico','Peru'),bty='n',cex=0.75,fill=c(2:3,4,5:8)) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure or major regimen change by site, accounting for censoring.} cum.inc.argentina_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'argentina')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.brazil_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'brazil')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.chile_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'chile')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'chile')], cencode='0',rho=0) cum.inc.haiti_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.honduras_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'honduras')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.mexico_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'mexico')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.peru_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor == 'peru')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor == 'peru')], cencode='0',rho=0) y.rge.site_new <- range(c(cum.inc.argentina_new[[1]]$est[which(cum.inc.argentina_new[[1]]$time <= 6.0)], cum.inc.brazil_new[[1]]$est[which(cum.inc.brazil_new[[1]]$time <= 6.0)], cum.inc.chile_new[[1]]$est[which(cum.inc.chile_new[[1]]$time <= 6.0)], cum.inc.haiti_new[[1]]$est[which(cum.inc.haiti_new[[1]]$time <= 6.0)], cum.inc.honduras_new[[1]]$est[which(cum.inc.honduras_new[[1]]$time <= 6.0)], cum.inc.mexico_new[[1]]$est[which(cum.inc.mexico_new[[1]]$time <= 6.0)], cum.inc.peru_new[[1]]$est[which(cum.inc.peru_new[[1]]$time <= 6.0)])) # pdf('Figure4_June130214.pdf') par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.argentina_new[[1]]$est[which(cum.inc.argentina_new[[1]]$time <= 6.0)] ~ cum.inc.argentina_new[[1]]$time[which(cum.inc.argentina_new[[1]]$time <= 6.0)], type='l',axes=FALSE, xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site_new,col=2,xlim=c(0,6)) lines(cum.inc.brazil_new[[1]]$est[which(cum.inc.brazil_new[[1]]$time <= 6.0)] ~ cum.inc.brazil_new[[1]]$time[which(cum.inc.brazil_new[[1]]$time <= 6.0)],col=3) lines(cum.inc.chile_new[[1]]$est[which(cum.inc.chile_new[[1]]$time <= 6.0)] ~ cum.inc.chile_new[[1]]$time[which(cum.inc.chile_new[[1]]$time <= 6.0)],col=4) lines(cum.inc.honduras_new[[1]]$est[which(cum.inc.honduras_new[[1]]$time <= 6.0)] ~ cum.inc.honduras_new[[1]]$time[which(cum.inc.honduras_new[[1]]$time <= 6.0)],col=6) lines(cum.inc.mexico_new[[1]]$est[which(cum.inc.mexico_new[[1]]$time <= 6.0)] ~ cum.inc.mexico_new[[1]]$time[which(cum.inc.mexico_new[[1]]$time <= 6.0)],col=7) lines(cum.inc.peru_new[[1]]$est[which(cum.inc.peru_new[[1]]$time <= 6.0)] ~ cum.inc.peru_new[[1]]$time[which(cum.inc.peru_new[[1]]$time <= 6.0)],lty=7,col=8) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru'),bty='n',cex=0.75,fill=c(2,3,4,6,7,8)) # dev.off() # \end{figure} # \begin{figure} # \caption{Incidence of virologic failure, ignoring single VL > 1000 criteria, or major regimen change by site, accounting for censoring.} cum.inc.argentina.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'argentina')], fstatus=d$censor_new[which(d$site.factor == 'argentina')], cencode='0',rho=0) cum.inc.brazil.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'brazil')], fstatus=d$censor_new[which(d$site.factor == 'brazil')], cencode='0',rho=0) cum.inc.chile.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'chile')], fstatus=d$censor_new[which(d$site.factor == 'chile')], cencode='0',rho=0) cum.inc.haiti.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'haiti')], fstatus=d$censor_new[which(d$site.factor == 'haiti')], cencode='0',rho=0) cum.inc.honduras.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'honduras')], fstatus=d$censor_new[which(d$site.factor == 'honduras')], cencode='0',rho=0) cum.inc.mexico.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'mexico')], fstatus=d$censor_new[which(d$site.factor == 'mexico')], cencode='0',rho=0) cum.inc.peru.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor == 'peru')], fstatus=d$censor_new[which(d$site.factor == 'peru')], cencode='0',rho=0) y.rge.site.no1000_new <- range(c(cum.inc.argentina.no1000_new[[1]]$est[which(cum.inc.argentina.no1000_new[[1]]$time <= 6.0)], cum.inc.brazil.no1000_new[[1]]$est[which(cum.inc.brazil.no1000_new[[1]]$time <= 6.0)], cum.inc.chile.no1000_new[[1]]$est[which(cum.inc.chile.no1000_new[[1]]$time <= 6.0)], cum.inc.haiti.no1000_new[[1]]$est[which(cum.inc.haiti.no1000_new[[1]]$time <= 6.0)], cum.inc.honduras.no1000_new[[1]]$est[which(cum.inc.honduras.no1000_new[[1]]$time <= 6.0)], cum.inc.mexico.no1000_new[[1]]$est[which(cum.inc.mexico.no1000_new[[1]]$time <= 6.0)], cum.inc.peru.no1000_new[[1]]$est[which(cum.inc.peru.no1000_new[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.argentina.no1000_new[[1]]$est[which(cum.inc.argentina.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.argentina.no1000_new[[1]]$time[which(cum.inc.argentina.no1000_new[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.site_new,col=2) lines(cum.inc.brazil.no1000_new[[1]]$est[which(cum.inc.brazil.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.brazil.no1000_new[[1]]$time[which(cum.inc.brazil.no1000_new[[1]]$time <= 6.0)],col=3) lines(cum.inc.chile.no1000_new[[1]]$est[which(cum.inc.chile.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.chile.no1000_new[[1]]$time[which(cum.inc.chile.no1000_new[[1]]$time <= 6.0)],col=4) lines(cum.inc.honduras.no1000_new[[1]]$est[which(cum.inc.honduras.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.honduras.no1000_new[[1]]$time[which(cum.inc.honduras.no1000_new[[1]]$time <= 6.0)],col=6) lines(cum.inc.mexico.no1000_new[[1]]$est[which(cum.inc.mexico.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.mexico.no1000_new[[1]]$time[which(cum.inc.mexico.no1000_new[[1]]$time <= 6.0)],col=7) lines(cum.inc.peru.no1000_new[[1]]$est[which(cum.inc.peru.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.peru.no1000_new[[1]]$time[which(cum.inc.peru.no1000_new[[1]]$time <= 6.0)],lty=7,col=8) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Argentina','Brazil','Chile','Honduras','Mexico','Peru'),bty='n',cex=0.75,fill=c(2,3,4,6,7,8)) # \end{figure} # \begin{figure} # \caption{Cumulative incidence of virologic failure, major regimen change, and the composite endpoint excluding Haiti for all outcomes, accounting for censoring.} cum.inc.failure_new <- cuminc(ftime=d.1000$total_followup_failure.1000_new[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor_failure.1000_new[which(d.1000$site.factor != 'haiti')],cencode='0',rho=0) # cum.inc.failure <- cuminc(ftime=d.1000$total_followup_failure.1000,fstatus=d.1000$censor_failure.1000,cencode='0',rho=0) cum.inc.reg_change_new <- cuminc(ftime=d.1000$total_followup_reg_change_new, fstatus=d.1000$censor_reg_change_new,cencode='0',rho=0) cum.inc.reg_change_no_haiti_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.reg_change_only_haiti_new <- cuminc(ftime=d.1000$total_followup_reg_change_new[which(d.1000$site.factor == 'haiti')], fstatus=d.1000$censor_reg_change_new[which(d.1000$site.factor == 'haiti')],cencode='0',rho=0) cum.inc.composite_outcome_new <- cuminc(ftime=d.1000$total_followup.1000_new[which(d.1000$site.factor != 'haiti')], fstatus=d.1000$censor.1000_new[which(d.1000$site.factor != 'haiti')], cencode='0',rho=0) y.rge_new <- range(c(cum.inc.failure_new[[1]]$est[which(cum.inc.failure_new[[1]]$time <= 6.0)], cum.inc.reg_change_no_haiti_new[[1]]$est[which(cum.inc.reg_change_no_haiti_new[[1]]$time <= 6.0)], cum.inc.composite_outcome_new[[1]]$est[which(cum.inc.composite_outcome_new[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure_new[[1]]$est[which(cum.inc.failure_new[[1]]$time <= 6.0)] ~ cum.inc.failure_new[[1]]$time[which(cum.inc.failure_new[[1]]$time <= 6.0)], type='l',axes=FALSE, xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge_new,lty=1) lines(cum.inc.reg_change_no_haiti_new[[1]]$est[which(cum.inc.reg_change_no_haiti_new[[1]]$time <= 6.0)] ~ cum.inc.reg_change_no_haiti_new[[1]]$time[which(cum.inc.reg_change_no_haiti_new[[1]]$time <= 6.0)], lty=2) lines(cum.inc.composite_outcome_new[[1]]$est[which(cum.inc.composite_outcome_new[[1]]$time <= 6.0)] ~ cum.inc.composite_outcome_new[[1]]$time[which(cum.inc.composite_outcome_new[[1]]$time <= 6.0)], lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Failure','Major regimen change','Composite outcome'),lty=c(1,2,3),cex=0.75) # \end{figure} # \begin{figure} # \caption{Cumulative incidence of virologic failure, ignoring the single VL > 1000 criteria, major regimen change, and the composite endpoint, accounting for censoring.} cum.inc.failure.no1000_new <- cuminc(ftime=d$total_followup_failure_new[which(d$site.factor != 'haiti')], fstatus=d$censor_failure_new[which(d$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.reg_change.no1000_new <- cuminc(ftime=d$total_followup_reg_change_new, fstatus=d$censor_reg_change_new,cencode='0',rho=0) cum.inc.reg_change.no1000_no_haiti_new <- cuminc(ftime=d$total_followup_reg_change_new[which(d$site.factor != 'haiti')], fstatus=d$censor_reg_change_new[which(d$site.factor != 'haiti')],cencode='0',rho=0) cum.inc.composite_outcome.no1000_new <- cuminc(ftime=d$total_followup_new[which(d$site.factor != 'haiti')], fstatus=d$censor_new[which(d$site.factor != 'haiti')], cencode='0',rho=0) y.rge.no1000_new <- range(c(cum.inc.failure.no1000_new[[1]]$est[which(cum.inc.failure.no1000_new[[1]]$time <= 6.0)], cum.inc.reg_change.no1000_no_haiti_new[[1]]$est[which(cum.inc.reg_change.no1000_no_haiti_new[[1]]$time <= 6.0)], cum.inc.composite_outcome.no1000_new[[1]]$est[which(cum.inc.composite_outcome.no1000_new[[1]]$time <= 6.0)])) par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) plot(cum.inc.failure.no1000_new[[1]]$est[which(cum.inc.failure.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.failure.no1000_new[[1]]$time[which(cum.inc.failure.no1000_new[[1]]$time <= 6.0)], type='l', axes=FALSE,xlab='Years from HAART initiation',ylab='Cumulative incidence',ylim=y.rge.no1000_new,lty=1) lines(cum.inc.reg_change.no1000_no_haiti_new[[1]]$est[which(cum.inc.reg_change.no1000_no_haiti_new[[1]]$time <= 6.0)] ~ cum.inc.reg_change.no1000_no_haiti_new[[1]]$time[which(cum.inc.reg_change.no1000_no_haiti_new[[1]]$time <= 6.0)],lty=2) lines(cum.inc.composite_outcome.no1000_new[[1]]$est[which(cum.inc.composite_outcome.no1000_new[[1]]$time <= 6.0)] ~ cum.inc.composite_outcome.no1000_new[[1]]$time[which(cum.inc.composite_outcome.no1000_new[[1]]$time <= 6.0)], lty=3) axis(1) axis(2) box(bty='l') legend('topleft',legend=c('Failure','Major regimen change','Composite outcome'),lty=c(1,2,3),cex=0.75) # \end{figure} # Failure junk <- subset(d.1000, site.factor != 'haiti') num_in_fu_failure_1_new <- length(which(junk$total_followup_failure.1000_new > 1.0)) tbl_1_failure <- table(junk$censor_failure.1000_new[which(junk$total_followup_failure.1000_new <= 1.0)]) num_w_evt_failure_1_new <- tbl_1_failure[2] num_w_cr_failure_1_new <- tbl_1_failure[3] num_end_fu_failure_1_new <- tbl_1_failure[1] for(i in 2:floor(max(junk$total_followup_failure.1000_new))){ assign(x=paste('num_in_fu_failure_',i,'_new',sep=''), value=length(which(junk$total_followup_failure.1000_new > i))) tbl <- table(junk$censor_failure.1000_new[which(junk$total_followup_failure.1000_new <= i & junk$total_followup_failure.1000_new > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_failure_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_failure_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_failure_',i,'_new',sep=''), value=tmp) } # Major regimen change num_in_fu_reg_change_1_new <- length(which(d.1000$total_followup_reg_change_new > 1.0)) tbl_1_reg_change <- table(d.1000$censor_reg_change_new[which(d.1000$total_followup_reg_change_new <= 1.0)]) num_w_evt_reg_change_1_new <- tbl_1_reg_change[2] num_w_cr_reg_change_1_new <- tbl_1_reg_change[3] num_end_fu_reg_change_1_new <- tbl_1_reg_change[1] for(i in 2:floor(max(d.1000$total_followup_reg_change_new))){ assign(x=paste('num_in_fu_reg_change_',i,'_new',sep=''), value=length(which(d.1000$total_followup_reg_change_new > i))) tbl <- table(d.1000$censor_reg_change_new[which(d.1000$total_followup_reg_change_new <= i & d.1000$total_followup_reg_change_new > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_reg_change_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_reg_change_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_reg_change_',i,'_new',sep=''), value=tmp) } # Composite outcome # Switched from d.1000 to junk to eliminate Haiti subjects from this num_in_fu_composite_outcome_1_new <- length(which(junk$total_followup.1000_new > 1.0)) tbl_1_composite_outcome <- table(junk$censor.1000_new[which(junk$total_followup.1000_new <= 1.0)]) num_w_evt_composite_outcome_1_new <- tbl_1_composite_outcome[2] num_w_cr_composite_outcome_1_new <- tbl_1_composite_outcome[3] num_end_fu_composite_outcome_1_new <- tbl_1_composite_outcome[1] for(i in 2:floor(max(junk$total_followup.1000_new))){ assign(x=paste('num_in_fu_composite_outcome_',i,'_new',sep=''), value=length(which(junk$total_followup.1000_new > i))) tbl <- table(junk$censor.1000_new[which(junk$total_followup.1000_new <= i & junk$total_followup.1000_new > (i-1))]) tmp <- tbl[which(names(tbl) == '1')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_evt_composite_outcome_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '2')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_w_cr_composite_outcome_',i,'_new',sep=''), value=tmp) tmp <- tbl[which(names(tbl) == '0')] if(length(tmp) == 0){ tmp <- 0 } assign(x=paste('num_end_fu_composite_outcome_',i,'_new',sep=''), value=tmp) } # Finding cumulative incidence and 95% CI for each curve in Figure 4 at 1, 3, and 5 years after HAART initiation # Failure failure_inc_1_new <- format(round(cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) failure_lci_1_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 1.0)][1] failure_uci_1_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 1.0)][1] failure_ci_1_new <- paste('(',round(failure_lci_1_new,3),', ',round(failure_uci_1_new,3),')',sep='') failure_inc_3_new <- format(round(cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) failure_lci_3_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 3.0)][1] failure_uci_3_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 3.0)][1] failure_ci_3_new <- paste('(',round(failure_lci_3_new,3),', ',round(failure_uci_3_new,3),')',sep='') failure_inc_5_new <- format(round(cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) failure_lci_5_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 5.0)][1] failure_uci_5_new <- cum.inc.failure_new[[1]]$est[which(round(cum.inc.failure_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure_new[[1]]$var)[which(round(cum.inc.failure_new[[1]]$time,1) == 5.0)][1] failure_ci_5_new <- paste('(',round(failure_lci_5_new,3),', ',round(failure_uci_5_new,3),')',sep='') # Failure, ignoring single VL > 1000 criteria failure_inc_1.no1000_new <- format(round(cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 1.0)][1],4),nsmall=3) failure_lci_1.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 1.0)][1]) failure_uci_1.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var [which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 1.0)][1]) failure_ci_1.no1000_new <- paste('(',round(failure_lci_1.no1000_new,3),', ',round(failure_uci_1.no1000_new,4),')',sep='') failure_inc_3.no1000_new <- format(round(cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) failure_lci_3.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var [which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 3.0)][1]) failure_uci_3.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var [which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 3.0)][1]) failure_ci_3.no1000_new <- paste('(',round(failure_lci_3.no1000_new,3),', ',round(failure_uci_3.no1000_new,3),')',sep='') failure_inc_5.no1000_new <- format(round(cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) failure_lci_5.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var [which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 5.0)][1]) failure_uci_5.no1000_new <- cum.inc.failure.no1000_new[[1]]$est[which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.failure.no1000_new[[1]]$var [which(round(cum.inc.failure.no1000_new[[1]]$time,1) == 5.0)][1]) failure_ci_5.no1000_new <- paste('(',round(failure_lci_5.no1000_new,3),', ',round(failure_uci_5.no1000_new,3),')',sep='') # Major regimen change # With all sites reg_change_inc_1_new <- format(round(cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_1_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 1.0)][1] reg_change_uci_1_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 1.0)][1] reg_change_ci_1_new <- paste('(',round(reg_change_lci_1_new,3),', ',round(reg_change_uci_1_new,3),')',sep='') reg_change_inc_3_new <- format(round(cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_3_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 3.0)][1] reg_change_uci_3_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 3.0)][1] reg_change_ci_3_new <- paste('(',round(reg_change_lci_3_new,3),', ',round(reg_change_uci_3_new,3),')',sep='') reg_change_inc_5_new <- format(round(cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_5_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 5.0)][1] reg_change_uci_5_new <- cum.inc.reg_change_new[[1]]$est[which(round(cum.inc.reg_change_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_new[[1]]$var)[which(round(cum.inc.reg_change_new[[1]]$time,1) == 5.0)][1] reg_change_ci_5_new <- paste('(',round(reg_change_lci_5_new,3),', ',round(reg_change_uci_5_new,3),')',sep='') # Without Haiti reg_change_inc_no_haiti_1_new <- format(round(cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_no_haiti_1_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 1.0)][1] reg_change_uci_no_haiti_1_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 1.0)][1] reg_change_ci_no_haiti_1_new <- paste('(',round(reg_change_lci_no_haiti_1_new,3),', ',round(reg_change_uci_no_haiti_1_new,3),')',sep='') reg_change_inc_no_haiti_3_new <- format(round(cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_no_haiti_3_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 3.0)][1] reg_change_uci_no_haiti_3_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 3.0)][1] reg_change_ci_no_haiti_3_new <- paste('(',round(reg_change_lci_no_haiti_3_new,3),', ',round(reg_change_uci_no_haiti_3_new,3),')',sep='') reg_change_inc_no_haiti_5_new <- format(round(cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_no_haiti_5_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 5.0)][1] reg_change_uci_no_haiti_5_new <- cum.inc.reg_change_no_haiti_new[[1]]$est[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_no_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_no_haiti_new[[1]]$time,1) == 5.0)][1] reg_change_ci_no_haiti_5_new <- paste('(',round(reg_change_lci_no_haiti_5_new,3),', ',round(reg_change_uci_no_haiti_5_new,3),')',sep='') # With only Haiti reg_change_inc_only_haiti_1_new <- format(round(cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) reg_change_lci_only_haiti_1_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 1.0)][1] reg_change_uci_only_haiti_1_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 1.0)][1] reg_change_ci_only_haiti_1_new <- paste('(',round(reg_change_lci_only_haiti_1_new,3),', ',round(reg_change_uci_only_haiti_1_new,3),')',sep='') reg_change_inc_only_haiti_3_new <- format(round(cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) reg_change_lci_only_haiti_3_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 3.0)][1] reg_change_uci_only_haiti_3_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 3.0)][1] reg_change_ci_only_haiti_3_new <- paste('(',round(reg_change_lci_only_haiti_3_new,3),', ',round(reg_change_uci_only_haiti_3_new,3),')',sep='') reg_change_inc_only_haiti_5_new <- format(round(cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) reg_change_lci_only_haiti_5_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 5.0)][1] reg_change_uci_only_haiti_5_new <- cum.inc.reg_change_only_haiti_new[[1]]$est[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.reg_change_only_haiti_new[[1]]$var)[which(round(cum.inc.reg_change_only_haiti_new[[1]]$time,1) == 5.0)][1] reg_change_ci_only_haiti_5_new <- paste('(',round(reg_change_lci_only_haiti_5_new,3),', ',round(reg_change_uci_only_haiti_5_new,3),')',sep='') # Composite outcome composite_outcome_inc_1_new <- format(round(cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) composite_outcome_lci_1_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 1.0)][1] composite_outcome_uci_1_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 1.0)][1] composite_outcome_ci_1_new <- paste('(',round(composite_outcome_lci_1_new,3),', ',round(composite_outcome_uci_1_new,3),')',sep='') composite_outcome_inc_3_new <- format(round(cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) composite_outcome_lci_3_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 3.0)][1] composite_outcome_uci_3_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 3.0)][1] composite_outcome_ci_3_new <- paste('(',round(composite_outcome_lci_3_new,3),', ',round(composite_outcome_uci_3_new,3),')',sep='') composite_outcome_inc_5_new <- format(round(cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) composite_outcome_lci_5_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 5.0)][1] composite_outcome_uci_5_new <- cum.inc.composite_outcome_new[[1]]$est[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome_new[[1]]$var)[which(round(cum.inc.composite_outcome_new[[1]]$time,1) == 5.0)][1] composite_outcome_ci_5_new <- paste('(',round(composite_outcome_lci_5_new,3),', ',round(composite_outcome_uci_5_new,3),')',sep='') # Composite outcome, ignoring single VL > 1000 criteria composite_outcome_inc_1.no1000_new <- format(round(cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 1.0)][1],3),nsmall=3) composite_outcome_lci_1.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 1.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 1.0)][1] composite_outcome_uci_1.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 1.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 1.0)][1] composite_outcome_ci_1.no1000_new <- paste('(',round(composite_outcome_lci_1.no1000_new,3),', ',round(composite_outcome_uci_1.no1000_new,3),')',sep='') composite_outcome_inc_3.no1000_new <- format(round(cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 3.0)][1],3),nsmall=3) composite_outcome_lci_3.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 3.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 3.0)][1] composite_outcome_uci_3.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 3.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 3.0)][1] composite_outcome_ci_3.no1000_new <- paste('(',round(composite_outcome_lci_3.no1000_new,3),', ',round(composite_outcome_uci_3.no1000_new,3),')',sep='') composite_outcome_inc_5.no1000_new <- format(round(cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 5.0)][1],3),nsmall=3) composite_outcome_lci_5.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 5.0)][1] - qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 5.0)][1] composite_outcome_uci_5.no1000_new <- cum.inc.composite_outcome.no1000_new[[1]]$est[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 5.0)][1] + qnorm(0.975)*sqrt(cum.inc.composite_outcome.no1000_new[[1]]$var)[which(round(cum.inc.composite_outcome.no1000_new[[1]]$time,1) == 5.0)][1] composite_outcome_ci_5.no1000_new <- paste('(',round(composite_outcome_lci_5.no1000_new,3),', ',round(composite_outcome_uci_5.no1000_new,3),')',sep='') fu_df_new <- data.frame('Year'=rep(seq(1,12),times=3), 'Had event'=c(num_w_evt_failure_1_new,num_w_evt_failure_2_new,num_w_evt_failure_3_new,num_w_evt_failure_4_new, num_w_evt_failure_5_new,num_w_evt_failure_6_new,num_w_evt_failure_7_new,num_w_evt_failure_8_new, num_w_evt_failure_9_new,num_w_evt_failure_10_new,num_w_evt_failure_11_new,num_w_evt_failure_12_new, num_w_evt_reg_change_1_new,num_w_evt_reg_change_2_new,num_w_evt_reg_change_3_new,num_w_evt_reg_change_4_new, num_w_evt_reg_change_5_new,num_w_evt_reg_change_6_new,num_w_evt_reg_change_7_new,num_w_evt_reg_change_8_new, num_w_evt_reg_change_9_new,num_w_evt_reg_change_10_new,num_w_evt_reg_change_11_new,num_w_evt_reg_change_12_new, num_w_evt_composite_outcome_1_new,num_w_evt_composite_outcome_2_new,num_w_evt_composite_outcome_3_new,num_w_evt_composite_outcome_4_new, num_w_evt_composite_outcome_5_new,num_w_evt_composite_outcome_6_new,num_w_evt_composite_outcome_7_new,num_w_evt_composite_outcome_8_new, num_w_evt_composite_outcome_9_new,num_w_evt_composite_outcome_10_new,num_w_evt_composite_outcome_11_new,num_w_evt_composite_outcome_12_new), 'Had competing risk'=c(num_w_cr_failure_1_new,num_w_cr_failure_2_new,num_w_cr_failure_3_new,num_w_cr_failure_4_new, num_w_cr_failure_5_new,num_w_cr_failure_6_new,num_w_cr_failure_7_new,num_w_cr_failure_8_new, num_w_cr_failure_9_new,num_w_cr_failure_10_new,num_w_cr_failure_11_new,num_w_cr_failure_12_new, num_w_cr_reg_change_1_new,num_w_cr_reg_change_2_new,num_w_cr_reg_change_3_new,num_w_cr_reg_change_4_new, num_w_cr_reg_change_5_new,num_w_cr_reg_change_6_new,num_w_cr_reg_change_7_new,num_w_cr_reg_change_8_new, num_w_cr_reg_change_9_new,num_w_cr_reg_change_10_new,num_w_cr_reg_change_11_new,num_w_cr_reg_change_12_new, num_w_cr_composite_outcome_1_new,num_w_cr_composite_outcome_2_new,num_w_cr_composite_outcome_3_new,num_w_cr_composite_outcome_4_new, num_w_cr_composite_outcome_5_new,num_w_cr_composite_outcome_6_new,num_w_cr_composite_outcome_7_new,num_w_cr_composite_outcome_8_new, num_w_cr_composite_outcome_9_new,num_w_cr_composite_outcome_10_new,num_w_cr_composite_outcome_11_new,num_w_cr_composite_outcome_12_new), 'Censored'=c(num_end_fu_failure_1_new,num_end_fu_failure_2_new,num_end_fu_failure_3_new,num_end_fu_failure_4_new, num_end_fu_failure_5_new,num_end_fu_failure_6_new,num_end_fu_failure_7_new,num_end_fu_failure_8_new, num_end_fu_failure_9_new,num_end_fu_failure_10_new,num_end_fu_failure_11_new,num_end_fu_failure_12_new, num_end_fu_reg_change_1_new,num_end_fu_reg_change_2_new,num_end_fu_reg_change_3_new,num_end_fu_reg_change_4_new, num_end_fu_reg_change_5_new,num_end_fu_reg_change_6_new,num_end_fu_reg_change_7_new,num_end_fu_reg_change_8_new, num_end_fu_reg_change_9_new,num_end_fu_reg_change_10_new,num_end_fu_reg_change_11_new,num_end_fu_reg_change_12_new, num_end_fu_composite_outcome_1_new,num_end_fu_composite_outcome_2_new,num_end_fu_composite_outcome_3_new,num_end_fu_composite_outcome_4_new, num_end_fu_composite_outcome_5_new,num_end_fu_composite_outcome_6_new,num_end_fu_composite_outcome_7_new,num_end_fu_composite_outcome_8_new, num_end_fu_composite_outcome_9_new,num_end_fu_composite_outcome_10_new,num_end_fu_composite_outcome_11_new,num_end_fu_composite_outcome_12_new), 'Still in follow-up'=c(num_in_fu_failure_1_new,num_in_fu_failure_2_new,num_in_fu_failure_3_new,num_in_fu_failure_4_new, num_in_fu_failure_5_new,num_in_fu_failure_6_new,num_in_fu_failure_7_new,num_in_fu_failure_8_new, num_in_fu_failure_9_new,num_in_fu_failure_10_new,num_in_fu_failure_11_new,num_in_fu_failure_12_new, num_in_fu_reg_change_1_new,num_in_fu_reg_change_2_new,num_in_fu_reg_change_3_new,num_in_fu_reg_change_4_new, num_in_fu_reg_change_5_new,num_in_fu_reg_change_6_new,num_in_fu_reg_change_7_new,num_in_fu_reg_change_8_new, num_in_fu_reg_change_9_new,num_in_fu_reg_change_10_new,num_in_fu_reg_change_11_new,num_in_fu_reg_change_12_new, num_in_fu_composite_outcome_1_new,num_in_fu_composite_outcome_2_new,num_in_fu_composite_outcome_3_new,num_in_fu_composite_outcome_4_new, num_in_fu_composite_outcome_5_new,num_in_fu_composite_outcome_6_new,num_in_fu_composite_outcome_7_new,num_in_fu_composite_outcome_8_new, num_in_fu_composite_outcome_9_new,num_in_fu_composite_outcome_10_new,num_in_fu_composite_outcome_11_new,num_in_fu_composite_outcome_12_new)) fu_df_new <- as.matrix(fu_df_new) colnames(fu_df_new) <- c('Year','Had event','Had competing risk','Censored','Still in follow-up') # latex(fu_df_new, file='', rgroup=c('Virologic failure','Major regimen change','Composite outcome'),n.rgroup=c(12,12,12),rowname=c(''),rowlabel='',where='!h',longtable=TRUE,caption='Number of subjects having the event, the competing risk, or were censored by year of follow-up. Virologic failure and the composite outcome do not include subjects from Haiti. This is in the cohort in which gaps of a year or more in virologic measurements are censored.') sub_d.1000 <- subset(d.1000, total_followup.1000_new > 0,select=c(age_fhaart, male.factor, aids_at_fhaart,aids_at_fhaart.factor, fhaart_regimen_class.factor,nrti_type.factor, baseline_cd4, baseline_log_rna, transmission_risk_cat.factor,transmission_risk_larger_cat.factor, total_followup.1000_new,total_followup_failure.1000_new,total_followup_reg_change_new, composite_outcome.1000_new,site,site.factor,year_fhaart, days_from_fhaart_to_vl_failure.1000_new,overall_first_failure.1000_new, days_from_fhaart_to_vl_failure_new,overall_first_failure_new,max_date,date_fhaart, days_from_fhaart_to_second_line_start_new,overall_second_line_new, year_gap_date)) # Changed the following two lines to exclude subjects from Haiti followup.median.1000_new <- format(round(median(d.1000$total_followup.1000_new[which(d.1000$site.factor != 'haiti')]),2),nsmall=2) followup.iqr.1000_new <- format(round(quantile(d.1000$total_followup.1000_new[which(d.1000$site.factor != 'haiti')],probs=c(0.25,0.75)),2),nsmall=2) t2evt_no_haiti_new <- time_to_evt_function(df=sub_d.1000,include_haiti=FALSE,which_outcome='composite',only_haiti=FALSE,reduced=FALSE,censor=TRUE) t2evt_w_haiti_new <- time_to_evt_function(df=sub_d.1000,include_haiti=TRUE,which_outcome='composite',only_haiti=FALSE,reduced=TRUE,censor=TRUE) t2evt_no_haiti_failure_new <- time_to_evt_function(df=sub_d.1000,include_haiti=FALSE,which_outcome='failure',only_haiti=FALSE,reduced=FALSE,censor=TRUE) t2evt_haiti_regimen_change_new <- time_to_evt_function(df=sub_d.1000, include_haiti=TRUE,which_outcome='regimen',only_haiti=FALSE,reduced=TRUE,censor=TRUE) t2evt_no_haiti_regimen_change_new <- time_to_evt_function(df=sub_d.1000, include_haiti=FALSE,which_outcome='regimen',only_haiti=FALSE,reduced=FALSE,censor=TRUE) t2evt_only_haiti_regimen_change_new <- time_to_evt_function(df=sub_d.1000, include_haiti=TRUE,which_outcome='regimen',only_haiti=TRUE,reduced=TRUE,censor=TRUE) t2evt_no_haiti_reduced_regimen_change_new <- time_to_evt_function(df=sub_d.1000, include_haiti=FALSE,which_outcome='regimen',only_haiti=FALSE,reduced=TRUE,censor=TRUE) sub_d <- subset(d, total_followup_new > 0, select=c(age_fhaart, male.factor, aids_at_fhaart,aids_at_fhaart.factor, fhaart_regimen_class.factor,nrti_type.factor, baseline_cd4, baseline_log_rna, transmission_risk_cat.factor,transmission_risk_larger_cat.factor, total_followup_new,total_followup_failure_new, total_followup_reg_change_new, composite_outcome_new,site,site.factor,year_fhaart, days_from_fhaart_to_vl_failure_new,overall_first_failure_new, days_from_fhaart_to_vl_failure_new,overall_first_failure_new,max_date,date_fhaart, days_from_fhaart_to_second_line_start_new,overall_second_line_new)) followup.median <- format(round(median(d$total_followup),2),nsmall=2) followup.iqr <- format(round(quantile(d$total_followup,probs=c(0.25,0.75)),2),nsmall=2) t2evt_no_haiti.no1000_new <- time_to_evt_function(df=sub_d,include_haiti=FALSE,which_outcome='composite.no1000',only_haiti=FALSE,reduced=FALSE,censor=TRUE) t2evt_no_haiti_failure.no1000_new <- time_to_evt_function(df=sub_d,include_haiti=FALSE,which_outcome='failure.no1000',only_haiti=FALSE,reduced=FALSE,censor=TRUE) haiti.df_new <- t2evt_w_haiti_new[['primary.df']][-(28:36),] # latex(haiti.df_new, file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and adjusting for age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,where='!h',size='small',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) haiti_imputed.df_new <- t2evt_w_haiti_new[['primary.df.impute']][-(28:36),] # latex(haiti_imputed.df_new, file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and adjusting for age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_haiti_regimen_change_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome, stratified by site, and including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_haiti_regimen_change_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome, stratified by site, and including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_regimen_change_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome, stratified by site, and excluding Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, log viral load, and probable route of infection.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_regimen_change_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome, stratified by site, and excluding Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, log viral load, and probable route of infection.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_only_haiti_regimen_change_new[['primary.df']],file='',caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome and including only data from Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'), n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_only_haiti_regimen_change_new[['primary.df.impute']],file='',caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome and including only data from Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'), n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_reduced_regimen_change_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the major regimen change endpoint as the outcome and excluding Haiti from the analysis. This analysis has a reduced list of covariates from the typical analysis excluding Haiti: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_reduced_regimen_change_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the major regimen change endpoint as the outcome and excluding Haiti from the analysis. This analysis has a reduced list of covariates from the typical analysis excluding Haiti: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, and year of HAART.',rowname=NULL,size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti.no1000_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti.no1000_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the composite endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure.no1000_new[['primary.df']], file='', caption='Cox regression (unadjusted and adjusted) with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small', where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # latex(t2evt_no_haiti_failure.no1000_new[['primary.df.impute']], file='', caption='Cox regression (unadjusted and adjusted) using imputed data with the first failure endpoint as the outcome, stratified by site, and not including Haiti in the analysis. Covariates included in the model are: age, sex, AIDS at HAART initiation, CD4 at HAART initiation, type of initial regimen, year of HAART, HIV-1 RNA (log10-transformed), and route of infection. The virologic failure definition does not include the VL > 1000 criteria.',rowname=NULL, size='small',where='!h',cgroup=c('','Unadjusted','Adjusted'),n.cgroup=c(1,3,3),col.just=c('l','r','r','r||','r','r','r')) # Creating Table 2 for manuscript table2_new <- data.frame('Covariate'=t2evt_no_haiti_new[['primary.df.impute']][,1], t2evt_no_haiti_failure_new[['primary.df.impute']][,5:7], t2evt_no_haiti_regimen_change_new[['primary.df.impute']][,5:7], t2evt_no_haiti_new[['primary.df.impute']][,5:7]) table2_new <- as.matrix(table2_new) colnames(table2_new) <- c('Covariate',rep(c('HR','95\\% CI','P'),times=3)) # latex(table2_new, file='', landscape=TRUE,cgroup=c('','Failure','Major regimen change','Composite outcome'),n.cgroup=c(1,3,3,3),caption='Cox regression using imputed data and stratified by site and excluding Haiti from the analysis, accounting for censoring.',exclude1=FALSE,long=TRUE,longtable=TRUE) # Creating Table 3 for manuscript table3_new <- data.frame('Covariate'=t2evt_only_haiti_regimen_change_new[['primary.df.impute']][,1], t2evt_only_haiti_regimen_change_new[['primary.df.impute']][,5:7], t2evt_no_haiti_reduced_regimen_change_new[['primary.df.impute']][,5:7]) table3_new <- as.matrix(table3_new) colnames(table3_new) <- c('Covariate',rep(c('HR','95\\% CI','P'),times=2)) # latex(table3_new, file='', landscape=TRUE,cgroup=c('','Only Haiti','Excluding Haiti'),n.cgroup=c(1,3,3),caption='Cox regression using imputed data and stratified by site running the reduced model with only Haiti and excluding Haiti, accounting for censoring.',longtable=TRUE) # Median frequency of viral load measurements per site -- # of viral load measurements dividied by the follow-up time after HAART initiation per person. Then compute median by site. # First need to add total_followup to rna rna$total_followup.1000 <- d.1000[match(rna$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] rna$site.factor <- art[match(rna$unique_id,art$unique_id,nomatch=NA),'site.factor'] tmp <- subset(rna, !is.na(total_followup.1000) & site.factor != 'haiti' & !is.na(rna_v) & (is.na(year_gap_date) | is.na(rna_d) | (!is.na(rna_d) & !is.na(year_gap_date) & rna_d <= year_gap_date))) tmp$unique_id.factor <- factor(tmp$unique_id, levels=unique(d.1000$unique_id)) num_vl_measurements <- tapply(tmp$rna_v, INDEX=tmp$unique_id, FUN=function(y){length(y[which(!is.na(y))])}) num_vl_measurements.df <- data.frame('unique_id'=names(num_vl_measurements), 'num_vl_measurements'=num_vl_measurements) tmp$num_vl_per_person <- num_vl_measurements.df[match(tmp$unique_id,num_vl_measurements.df$unique_id,nomatch=NA),'num_vl_measurements'] tmp$vl_frequency <- tmp$num_vl_per_person/tmp$total_followup.1000 all_visits$vl_frequency <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'vl_frequency'] all_visits$num_vl_per_person <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'num_vl_per_person'] vl_frequency_by_site <- tapply(all_visits$vl_frequency[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) vl_frequency_overall <- quantile(all_visits$vl_frequency[!duplicated(all_visits$unique_id)], probs=c(0.25,0.50,0.75),na.rm=TRUE) vl_per_person_by_site <- tapply(all_visits$num_vl_per_person[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) # Median frequency of CD4 measurements per site -- # of CD4 measurements dividied by the follow-up time after HAART initiation per person. Then compute median by site. # First need to add total_followup to cd4 cd4$total_followup.1000 <- d.1000[match(cd4$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] tmp <- subset(cd4, !is.na(total_followup.1000) & !is.na(cd4_v) & (is.na(year_gap_date) | is.na(cd4_d) | (!is.na(cd4_d) & !is.na(year_gap_date) & cd4_d <= year_gap_date))) tmp$unique_id.factor <- factor(tmp$unique_id, levels=unique(d.1000$unique_id)) num_cd4_measurements <- tapply(tmp$cd4_v, INDEX=tmp$unique_id, FUN=function(y){length(y[which(!is.na(y))])}) num_cd4_measurements.df <- data.frame('unique_id'=names(num_cd4_measurements), 'num_cd4_measurements'=num_cd4_measurements) tmp$num_cd4_per_person <- num_cd4_measurements.df[match(tmp$unique_id,num_cd4_measurements.df$unique_id,nomatch=NA),'num_cd4_measurements'] tmp$cd4_frequency <- tmp$num_cd4_per_person/tmp$total_followup.1000 all_visits$cd4_frequency <- tmp[match(all_visits$unique_id,tmp$unique_id,nomatch=NA),'cd4_frequency'] cd4_frequency_by_site <- tapply(all_visits$cd4_frequency[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) cd4_frequency_overall <- quantile(all_visits$cd4_frequency[!duplicated(all_visits$unique_id)], probs=c(0.25,0.50,0.75),na.rm=TRUE) all_visits$num_cd4_per_person <- tmp[match(all_visits$unique_id,cd4$unique_id,nomatch=NA),'num_cd4_per_person'] cd4_per_person_by_site <- tapply(all_visits$num_cd4_per_person[!duplicated(all_visits$unique_id)],INDEX=all_visits$site.factor[!duplicated(all_visits$unique_id)],FUN=median,na.rm=TRUE) # For median follow-up time (IQR), calculating for the entire all_visits cohort (6598) and then for those in d.1000 (6510 -- excludes the 88 with 0 follow-up) all_visits$total_followup.1000 <- d.1000[match(all_visits$unique_id,d.1000$unique_id,nomatch=NA),'total_followup.1000'] all_visits$total_followup.1000 <- ifelse(is.na(all_visits$total_followup.1000),0,all_visits$total_followup.1000) entire_fu <- describe(all_visits$total_followup.1000[!duplicated(all_visits$unique_id)]) d.1000_fu <- describe(d.1000$total_followup.1000[!duplicated(d.1000$unique_id)]) # Creating df for CD4 and VL frequencies by site df_site_new <- data.frame('Site'=names(cd4_frequency_by_site), 'By year'=round(cd4_frequency_by_site,2), 'Overall'=round(cd4_per_person_by_site,2), 'By year'=round(vl_frequency_by_site,2), 'Overall'=round(vl_per_person_by_site,2)) df_site_new <- as.matrix(df_site_new) colnames(df_site_new) <- c('Site','By year','Overall','By year','Overall') df_site_new[,4] <- ifelse(is.na(df_site_new[,4]),'',df_site_new[,4]) df_site_new[,5] <- ifelse(is.na(df_site_new[,5]),'',df_site_new[,5]) overall_cd4 <- format(round(cd4_frequency_overall[2],2),nsmall=2) overall_vl <- format(round(vl_frequency_overall[2],2),nsmall=2) # latex(df_site_new, rowname=NULL, file='', caption='Frequency of CD4 and viral load measurements by site, accounting for censoring.',where='!h',cgroup=c('','CD4','HIV-1 RNA'),n.cgroup=c(1,2,2), # insert.bottom=paste('Median rate of CD4 measurement frequency=',overall_cd4,'; Median rate of HIV-1 RNA measurement frequency=',overall_vl,'.',sep='')) # Checking manuscript mw <- subset(rna, unique_id %in% unique(d.1000$unique_id[which(d.1000$overall_first_failure.1000 == 1)])) mw$date_first_failure.1000 <- d.1000[match(mw$unique_id, d.1000$unique_id,nomatch=NA),'date_first_failure.1000'] mw$overall_second_line <- d.1000[match(mw$unique_id, d.1000$unique_id,nomatch=NA),'overall_second_line'] mw <- subset(mw, rna_d > date_first_failure.1000 & overall_second_line == 0) # \section{Appendix - R/Package Versions} # \small version['version.string'] pack <- installed.packages() pack.out <- pack[,c('Package','Version','Priority','Depends')] pack.in.session <- (.packages()) pack.out <- data.frame(pack.out[pack.out[,1] %in% pack.in.session,])[,-1] pack.out[!pack.out$Priority %in% c('base','recommended'),-2,drop=FALSE]