output_mgmt_no_hiv <- function(mod, le5_ind){ if(le5_ind){ age_comp1 <- c(2.0,2.0,0.5) age_comp2 <- c(2.0, 2.0, 1.0) age_comp3 <- c(2.0, 2.0, 5.0) } else { age_comp1 <- c(9.0,9.0,6.0) age_comp2 <- c(9.0,9.0,12.0) age_comp3 <- c(9.0,9.0,15.0) } summ_mod_1 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2002),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') summ_mod_2 <- summary(mod, age_fhaart=age_comp2, male.factor='Female',baseline_cd4=c(350,350,200), year_fhaart=c(2006,2006,2004),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') summ_mod_3 <- summary(mod, age_fhaart=age_comp3, male.factor='Female',baseline_cd4=c(350,350,500), year_fhaart=c(2006,2006,2006),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') summ_mod_4 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2008),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') summ_mod_5 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2010),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') summ_mod_6 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2012),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical='not AIDS') order.wanted <- c(which(rownames(summ_mod_1) == 'age_fhaart'), which(rownames(summ_mod_1) == 'male.factor - Male:Female'), which(rownames(summ_mod_1) == 'baseline_cd4'), which(rownames(summ_mod_1) == 'year_fhaart'), which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - AIDS:not AIDS'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - Unknown:not AIDS')) order.wanted <- order.wanted + 1 hr <- c(summ_mod_1[order.wanted[1],'Effect'],summ_mod_2[order.wanted[1],'Effect'], summ_mod_3[order.wanted[1],'Effect'], summ_mod_1[order.wanted[2],'Effect'], summ_mod_1[order.wanted[3],'Effect'],summ_mod_2[order.wanted[3],'Effect'],summ_mod_3[order.wanted[3],'Effect'], summ_mod_1[order.wanted[4],'Effect'],summ_mod_2[order.wanted[4],'Effect'],summ_mod_4[order.wanted[4],'Effect'], summ_mod_5[order.wanted[4],'Effect'],summ_mod_6[order.wanted[4],'Effect'], summ_mod_1[order.wanted[5],'Effect'],summ_mod_1[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect']) hr <- format(round(hr,2),nsmall=2) lci <- c(summ_mod_1[order.wanted[1],'Lower 0.95'],summ_mod_2[order.wanted[1],'Lower 0.95'], summ_mod_3[order.wanted[1],'Lower 0.95'], summ_mod_1[order.wanted[2],'Lower 0.95'], summ_mod_1[order.wanted[3],'Lower 0.95'],summ_mod_2[order.wanted[3],'Lower 0.95'],summ_mod_3[order.wanted[3],'Lower 0.95'], summ_mod_1[order.wanted[4],'Lower 0.95'],summ_mod_2[order.wanted[4],'Lower 0.95'],summ_mod_4[order.wanted[4],'Lower 0.95'], summ_mod_5[order.wanted[4],'Lower 0.95'],summ_mod_6[order.wanted[4],'Lower 0.95'], summ_mod_1[order.wanted[5],'Lower 0.95'],summ_mod_1[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95']) lci <- format(round(lci,2),nsmall=2) uci <- c(summ_mod_1[order.wanted[1],'Upper 0.95'],summ_mod_2[order.wanted[1],'Upper 0.95'], summ_mod_3[order.wanted[1],'Upper 0.95'], summ_mod_1[order.wanted[2],'Upper 0.95'], summ_mod_1[order.wanted[3],'Upper 0.95'],summ_mod_2[order.wanted[3],'Upper 0.95'],summ_mod_3[order.wanted[3],'Upper 0.95'], summ_mod_1[order.wanted[4],'Upper 0.95'],summ_mod_2[order.wanted[4],'Upper 0.95'],summ_mod_4[order.wanted[4],'Upper 0.95'], summ_mod_5[order.wanted[4],'Upper 0.95'],summ_mod_6[order.wanted[4],'Upper 0.95'], summ_mod_1[order.wanted[5],'Upper 0.95'],summ_mod_1[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95']) uci <- format(round(uci,2),nsmall=2) ci <- paste('(',lci, ', ',uci, ')',sep='') mod_anova <- anova(mod) order.wanted.p <- c(which(rownames(mod_anova) == 'age_fhaart'), which(rownames(mod_anova) == 'male.factor'), which(rownames(mod_anova) == 'baseline_cd4'), which(rownames(mod_anova) == 'year_fhaart'), which(rownames(mod_anova) == 'fhaart_class_dichot.factor'), which(rownames(mod_anova) == 'pre_art_aids.clinical')) mod_anova_p <- mod_anova[order.wanted.p, 'P'] mod_anova_p <- ifelse(mod_anova_p < 0.001, '< 0.001', ifelse(mod_anova_p >= 0.001 & mod_anova_p < 0.01, round(mod_anova_p,3),round(mod_anova_p,2))) if(le5_ind){ age_label <- c('Baseline age (years)','~~~~~0.5','~~~~~1.0','~~~~~2.0 (ref)','~~~~~5.0') age_hr <- c(hr[1:2],'1.00',hr[3]) age_ci <- c(ci[1:2],'1.00',ci[3]) } else { age_label <- c('Baseline age (years)','~~~~~6.0','~~~~~9.0 (ref)','~~~~~12.0','~~~~~15.0') age_hr <- c(hr[1],'1.00',hr[2:3]) age_ci <- c(ci[1],'1.00',ci[2:3]) } df <- data.frame('Covariate'=c(age_label,'', 'Sex','~~~~~Male','~~~~~Female (ref)','', 'Baseline CD4 (square root transformed)','~~~~~100','~~~~~200','~~~~~350 (ref)','~~~~~500','', 'Year of first HAART','~~~~~2002','~~~~~2004','~~~~~2006 (ref)','~~~~~2008','~~~~~2010','~~~~~2012','', 'First HAART regimen class','~~~~~PI-only (ref)','~~~~~non-PI','', 'Clinical AIDS at baseline','~~~~~not AIDS (ref)','~~~~~AIDS','~~~~~Unknown'), 'HR'=c('',age_hr,'', '',hr[4],'1.00','', '',hr[5:6],'1.00',hr[7],'', '',hr[8:9],'1.00',hr[10:12],'', '','1.00',hr[13],'', '','1.00',hr[14:15]), '95\\% CI'=c('',age_ci,'', '',ci[4],'','', '',ci[5:6],'',ci[7],'', '',ci[8:9],'',ci[10:12],'', '','',ci[13],'', '','',ci[14:15]), 'P'=c(mod_anova_p[1],rep('',times=5), mod_anova_p[2],rep('',times=3), mod_anova_p[3],rep('',times=5), mod_anova_p[4],rep('',times=7), mod_anova_p[5],rep('',times=3), mod_anova_p[6],rep('',times=3)),check.names=FALSE) return(df) } output_mgmt_cd4_per <- function(mod, le5_ind){ if(le5_ind){ age_comp1 <- c(2.0,2.0,0.5) age_comp2 <- c(2.0, 2.0, 1.0) age_comp3 <- c(2.0, 2.0, 5.0) } else { age_comp1 <- c(9.0,9.0,6.0) age_comp2 <- c(9.0,9.0,12.0) age_comp3 <- c(9.0,9.0,15.0) } summ_mod_1 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,10), year_fhaart=c(2006,2006,2002),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_2 <- summary(mod, age_fhaart=age_comp2, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2004),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','Unknown')) summ_mod_3 <- summary(mod, age_fhaart=age_comp3, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2006),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_4 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2008),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_5 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2010),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_6 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2012),fhaart_class_dichot.factor='PI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) order.wanted <- c(which(rownames(summ_mod_1) == 'age_fhaart'), which(rownames(summ_mod_1) == 'male.factor - Male:Female'), which(rownames(summ_mod_1) == 'baseline_cd4_per'), which(rownames(summ_mod_1) == 'year_fhaart'), which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - AIDS:not AIDS'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - Unknown:not AIDS')) order.wanted <- order.wanted + 1 hr <- c(summ_mod_1[order.wanted[1],'Effect'],summ_mod_2[order.wanted[1],'Effect'], summ_mod_3[order.wanted[1],'Effect'], summ_mod_1[order.wanted[2],'Effect'], summ_mod_1[order.wanted[3],'Effect'],summ_mod_2[order.wanted[3],'Effect'], summ_mod_1[order.wanted[4],'Effect'],summ_mod_2[order.wanted[4],'Effect'],summ_mod_4[order.wanted[4],'Effect'], summ_mod_5[order.wanted[4],'Effect'],summ_mod_6[order.wanted[4],'Effect'], summ_mod_1[order.wanted[5],'Effect'],summ_mod_1[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect']) hr <- format(round(hr,2),nsmall=2,scientific=FALSE) lci <- c(summ_mod_1[order.wanted[1],'Lower 0.95'],summ_mod_2[order.wanted[1],'Lower 0.95'], summ_mod_3[order.wanted[1],'Lower 0.95'], summ_mod_1[order.wanted[2],'Lower 0.95'], summ_mod_1[order.wanted[3],'Lower 0.95'],summ_mod_2[order.wanted[3],'Lower 0.95'], summ_mod_1[order.wanted[4],'Lower 0.95'],summ_mod_2[order.wanted[4],'Lower 0.95'],summ_mod_4[order.wanted[4],'Lower 0.95'], summ_mod_5[order.wanted[4],'Lower 0.95'],summ_mod_6[order.wanted[4],'Lower 0.95'], summ_mod_1[order.wanted[5],'Lower 0.95'],summ_mod_1[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95']) lci <- format(round(lci,2),nsmall=2,scientific=FALSE) uci <- c(summ_mod_1[order.wanted[1],'Upper 0.95'],summ_mod_2[order.wanted[1],'Upper 0.95'], summ_mod_3[order.wanted[1],'Upper 0.95'], summ_mod_1[order.wanted[2],'Upper 0.95'], summ_mod_1[order.wanted[3],'Upper 0.95'],summ_mod_2[order.wanted[3],'Upper 0.95'], summ_mod_1[order.wanted[4],'Upper 0.95'],summ_mod_2[order.wanted[4],'Upper 0.95'],summ_mod_4[order.wanted[4],'Upper 0.95'], summ_mod_5[order.wanted[4],'Upper 0.95'],summ_mod_6[order.wanted[4],'Upper 0.95'], summ_mod_1[order.wanted[5],'Upper 0.95'],summ_mod_1[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95']) uci <- format(round(uci,2),nsmall=2,scientific=FALSE) ci <- paste('(',lci, ', ',uci, ')',sep='') mod_anova <- anova(mod) order.wanted.p <- c(which(rownames(mod_anova) == 'age_fhaart'), which(rownames(mod_anova) == 'male.factor'), which(rownames(mod_anova) == 'baseline_cd4_per'), which(rownames(mod_anova) == 'year_fhaart'), which(rownames(mod_anova) == 'fhaart_class_dichot.factor'), which(rownames(mod_anova) == 'pre_art_aids.clinical')) mod_anova_p <- mod_anova[order.wanted.p, 'P'] mod_anova_p <- ifelse(mod_anova_p < 0.001, '< 0.001', ifelse(mod_anova_p >= 0.001 & mod_anova_p < 0.01, round(mod_anova_p,3),round(mod_anova_p,2))) if(le5_ind){ age_label <- c('Baseline age (years)','~~~~~0.5','~~~~~1.0','~~~~~2.0 (ref)','~~~~~5.0') age_hr <- c(hr[1:2],'1.00',hr[3]) age_ci <- c(ci[1:2],'1.00',ci[3]) } else { age_label <- c('Baseline age (years)','~~~~~6.0','~~~~~9.0 (ref)','~~~~~12.0','~~~~~15.0') age_hr <- c(hr[1],'1.00',hr[2:3]) age_ci <- c(ci[1],'1.00',ci[2:3]) } df <- data.frame('Covariate'=c(age_label,'', 'Sex','~~~~~Male','~~~~~Female (ref)','', 'Baseline CD4 percent (square root transformed)','~~~~~10','~~~~~25 (ref)','~~~~~30','', 'Year of first HAART','~~~~~2002','~~~~~2004','~~~~~2006 (ref)','~~~~~2008','~~~~~2010','~~~~~2012','', 'First HAART regimen class','~~~~~PI-only (ref)','~~~~~non-PI','', 'Clinical AIDS at baseline','~~~~~not AIDS (ref)','~~~~~AIDS','~~~~~Unknown'), 'HR'=c('',age_hr,'', '',hr[4],'1.00','', '',hr[5],'1.00',hr[6],'', '',hr[7:8],'1.00',hr[9:11],'', '','1.00',hr[12],'', '','1.00',hr[13:14]), '95\\% CI'=c('',age_ci,'', '',ci[4],'','', '',ci[5],'',ci[6],'', '',ci[7:8],'',ci[9:11],'', '','',ci[12],'', '','',ci[13:14]), 'P'=c(mod_anova_p[1],rep('',times=5), mod_anova_p[2],rep('',times=3), mod_anova_p[3],rep('',times=4), mod_anova_p[4],rep('',times=7), mod_anova_p[5],rep('',times=3), mod_anova_p[6],rep('',times=3)),check.names=FALSE) return(df) } output_mgmt <- function(mod, le5_ind, vl_restrict){ if(le5_ind){ age_comp1 <- c(2.0,2.0,0.5) age_comp2 <- c(2.0, 2.0, 1.0) age_comp3 <- c(2.0, 2.0, 5.0) } else { age_comp1 <- c(9.0,9.0,6.0) age_comp2 <- c(9.0,9.0,12.0) age_comp3 <- c(9.0,9.0,15.0) } summ_mod_1 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2002),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_2 <- summary(mod, age_fhaart=age_comp2, male.factor='Female',baseline_cd4=c(350,350,200), year_fhaart=c(2006,2006,2004),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,6),pre_art_aids.clinical=c('not AIDS','not AIDS','Unknown')) summ_mod_3 <- summary(mod, age_fhaart=age_comp3, male.factor='Female',baseline_cd4=c(350,350,500), year_fhaart=c(2006,2006,2006),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_4 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2008),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_5 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2010),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_6 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4=c(350,350,100), year_fhaart=c(2006,2006,2012),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) if(vl_restrict == 0){ fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only')) } else { # fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class.factor - PI-only or Other:NNRTI-only')) fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only')) } order.wanted <- c(which(rownames(summ_mod_1) == 'age_fhaart'), which(rownames(summ_mod_1) == 'male.factor - Male:Female'), which(rownames(summ_mod_1) == 'baseline_cd4'), which(rownames(summ_mod_1) == 'year_fhaart'), fhaart_rows, which(rownames(summ_mod_1) == 'log10_baseline_rna'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - AIDS:not AIDS'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - Unknown:not AIDS')) order.wanted <- order.wanted + 1 if(vl_restrict == 0){ hr_rows <- c(summ_mod_1[order.wanted[5],'Effect'],summ_mod_1[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect'],summ_mod_2[order.wanted[7],'Effect'], summ_mod_1[order.wanted[8],'Effect']) lci_rows <- c(summ_mod_1[order.wanted[5],'Lower 0.95'],summ_mod_1[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95'],summ_mod_2[order.wanted[7],'Lower 0.95'], summ_mod_1[order.wanted[8],'Lower 0.95']) uci_rows <- c(summ_mod_1[order.wanted[5],'Upper 0.95'],summ_mod_1[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95'],summ_mod_2[order.wanted[7],'Upper 0.95'], summ_mod_1[order.wanted[8],'Upper 0.95']) } else { hr_rows <- c(summ_mod_1[order.wanted[5],'Effect'], summ_mod_1[order.wanted[6],'Effect'],summ_mod_2[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect'],summ_mod_1[order.wanted[8],'Effect']) lci_rows <- c(summ_mod_1[order.wanted[5],'Lower 0.95'], summ_mod_1[order.wanted[6],'Lower 0.95'],summ_mod_2[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95'],summ_mod_1[order.wanted[8],'Lower 0.95']) uci_rows <- c(summ_mod_1[order.wanted[5],'Upper 0.95'], summ_mod_1[order.wanted[6],'Upper 0.95'],summ_mod_2[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95'],summ_mod_1[order.wanted[8],'Upper 0.95']) } hr <- c(summ_mod_1[order.wanted[1],'Effect'],summ_mod_2[order.wanted[1],'Effect'], summ_mod_3[order.wanted[1],'Effect'], summ_mod_1[order.wanted[2],'Effect'], summ_mod_1[order.wanted[3],'Effect'],summ_mod_2[order.wanted[3],'Effect'],summ_mod_3[order.wanted[3],'Effect'], summ_mod_1[order.wanted[4],'Effect'],summ_mod_2[order.wanted[4],'Effect'],summ_mod_4[order.wanted[4],'Effect'], summ_mod_5[order.wanted[4],'Effect'],summ_mod_6[order.wanted[4],'Effect'], hr_rows) hr <- format(round(hr,2),nsmall=2,scientific=FALSE) lci <- c(summ_mod_1[order.wanted[1],'Lower 0.95'],summ_mod_2[order.wanted[1],'Lower 0.95'], summ_mod_3[order.wanted[1],'Lower 0.95'], summ_mod_1[order.wanted[2],'Lower 0.95'], summ_mod_1[order.wanted[3],'Lower 0.95'],summ_mod_2[order.wanted[3],'Lower 0.95'],summ_mod_3[order.wanted[3],'Lower 0.95'], summ_mod_1[order.wanted[4],'Lower 0.95'],summ_mod_2[order.wanted[4],'Lower 0.95'],summ_mod_4[order.wanted[4],'Lower 0.95'], summ_mod_5[order.wanted[4],'Lower 0.95'],summ_mod_6[order.wanted[4],'Lower 0.95'], lci_rows) lci <- format(round(lci,2),nsmall=2,scientific=FALSE) uci <- c(summ_mod_1[order.wanted[1],'Upper 0.95'],summ_mod_2[order.wanted[1],'Upper 0.95'], summ_mod_3[order.wanted[1],'Upper 0.95'], summ_mod_1[order.wanted[2],'Upper 0.95'], summ_mod_1[order.wanted[3],'Upper 0.95'],summ_mod_2[order.wanted[3],'Upper 0.95'],summ_mod_3[order.wanted[3],'Upper 0.95'], summ_mod_1[order.wanted[4],'Upper 0.95'],summ_mod_2[order.wanted[4],'Upper 0.95'],summ_mod_4[order.wanted[4],'Upper 0.95'], summ_mod_5[order.wanted[4],'Upper 0.95'],summ_mod_6[order.wanted[4],'Upper 0.95'], uci_rows) uci <- format(round(uci,2),nsmall=2) ci <- paste('(',lci, ', ',uci, ')',sep='') mod_anova <- anova(mod) order.wanted.p <- c(which(rownames(mod_anova) == 'age_fhaart'), which(rownames(mod_anova) == 'male.factor'), which(rownames(mod_anova) == 'baseline_cd4'), which(rownames(mod_anova) == 'year_fhaart'), which(rownames(mod_anova) == 'fhaart_class_dichot.factor'), which(rownames(mod_anova) == 'log10_baseline_rna'), which(rownames(mod_anova) == 'pre_art_aids.clinical')) mod_anova_p <- mod_anova[order.wanted.p, 'P'] mod_anova_p <- ifelse(mod_anova_p < 0.001, '< 0.001', ifelse(mod_anova_p >= 0.001 & mod_anova_p < 0.01, round(mod_anova_p,3),round(mod_anova_p,2))) if(le5_ind){ age_label <- c('Baseline age (years)','~~~~~0.5','~~~~~1.0','~~~~~2.0 (ref)','~~~~~5.0') age_hr <- c(hr[1:2],'1.00',hr[3]) age_ci <- c(ci[1:2],'',ci[3]) } else { age_label <- c('Baseline age (years)','~~~~~6.0','~~~~~9.0 (ref)','~~~~~12.0','~~~~~15.0') age_hr <- c(hr[1],'1.00',hr[2:3]) age_ci <- c(ci[1],'',ci[2:3]) } if(vl_restrict == 0){ fhaart_label <- c('First HAART regimen class','~~~~~PI-only (ref)','~~~~~non-PI') fhaart_hr <- c('1.00',hr[13],'', '',hr[14],'1.00',hr[15],'', '','1.00',hr[16:17]) fhaart_ci <- c('',ci[13],'', '',ci[14],'',ci[15],'', '','',ci[16:17]) fhaart_times <- 3 } else { fhaart_label <- c('First HAART regimen class','~~~~~NNRTI-only (ref)','~~~~~PI-only or Other') fhaart_hr <- c('1.00',hr[13],'', '',hr[14],'1.00',hr[15],'', '','1.00',hr[16:17]) fhaart_ci <- c('',ci[13],'', '',ci[14],'',ci[15],'', '','',ci[16:17]) fhaart_times <- 3 } df <- data.frame('Covariate'=c(age_label,'', 'Sex','~~~~~Male','~~~~~Female (ref)','', 'Baseline CD4 (square root transformed)','~~~~~100','~~~~~200','~~~~~350 (ref)','~~~~~500','', 'Year of first HAART','~~~~~2002','~~~~~2004','~~~~~2006 (ref)','~~~~~2008','~~~~~2010','~~~~~2012','', fhaart_label,'', 'Baseline HIV-1 RNA (log-transformed)','~~~~~4','~~~~~5 (ref)','~~~~~6','', 'Clinical AIDS at baseline','~~~~~not AIDS (ref)','~~~~~AIDS','~~~~~Unknown'), 'HR'=c('',age_hr,'', '',hr[4],'1.00','', '',hr[5:6],'1.00',hr[7],'', '',hr[8:9],'1.00',hr[10:12],'', '',fhaart_hr), '95\\% CI'=c('',age_ci,'', '',ci[4],'','', '',ci[5:6],'',ci[7],'', '',ci[8:9],'',ci[10:12],'', '',fhaart_ci), 'P'=c(mod_anova_p[1],rep('',times=5), mod_anova_p[2],rep('',times=3), mod_anova_p[3],rep('',times=5), mod_anova_p[4],rep('',times=7), mod_anova_p[5],rep('',times=fhaart_times), mod_anova_p[6],rep('',times=4), mod_anova_p[7],rep('',times=3)),check.names=FALSE) return(df) } output_mgmt_vl_cd4_per <- function(mod, le5_ind,vl_restrict){ if(le5_ind){ age_comp1 <- c(2.0,2.0,0.5) age_comp2 <- c(2.0, 2.0, 1.0) age_comp3 <- c(2.0, 2.0, 5.0) } else { age_comp1 <- c(9.0,9.0,6.0) age_comp2 <- c(9.0,9.0,12.0) age_comp3 <- c(9.0,9.0,15.0) } summ_mod_1 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,10), year_fhaart=c(2006,2006,2002),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_2 <- summary(mod, age_fhaart=age_comp2, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2004),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,6),pre_art_aids.clinical=c('not AIDS','not AIDS','Unknown')) summ_mod_3 <- summary(mod, age_fhaart=age_comp3, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2006),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_4 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2008),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_5 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2010),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) summ_mod_6 <- summary(mod, age_fhaart=age_comp1, male.factor='Female',baseline_cd4_per=c(25,25,30), year_fhaart=c(2006,2006,2012),fhaart_class_dichot.factor='PI-only', log10_baseline_rna=c(5,5,4),pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS')) if(vl_restrict == 0){ fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only')) } else { # fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class.factor - PI-only or Other:NNRTI-only')) fhaart_rows <- c(which(rownames(summ_mod_1) == 'fhaart_class_dichot.factor - non-PI:PI-only')) } order.wanted <- c(which(rownames(summ_mod_1) == 'age_fhaart'), which(rownames(summ_mod_1) == 'male.factor - Male:Female'), which(rownames(summ_mod_1) == 'baseline_cd4_per'), which(rownames(summ_mod_1) == 'year_fhaart'), fhaart_rows, which(rownames(summ_mod_1) == 'log10_baseline_rna'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - AIDS:not AIDS'), which(rownames(summ_mod_1) == 'pre_art_aids.clinical - Unknown:not AIDS')) order.wanted <- order.wanted + 1 if(vl_restrict == 0){ hr_rows <- c(summ_mod_1[order.wanted[5],'Effect'],summ_mod_1[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect'],summ_mod_2[order.wanted[7],'Effect'], summ_mod_1[order.wanted[8],'Effect']) lci_rows <- c(summ_mod_1[order.wanted[5],'Lower 0.95'],summ_mod_1[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95'],summ_mod_2[order.wanted[7],'Lower 0.95'], summ_mod_1[order.wanted[8],'Lower 0.95']) uci_rows <- c(summ_mod_1[order.wanted[5],'Upper 0.95'],summ_mod_1[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95'],summ_mod_2[order.wanted[7],'Upper 0.95'], summ_mod_1[order.wanted[8],'Upper 0.95']) } else { hr_rows <- c(summ_mod_1[order.wanted[5],'Effect'], summ_mod_1[order.wanted[6],'Effect'],summ_mod_2[order.wanted[6],'Effect'], summ_mod_1[order.wanted[7],'Effect'],summ_mod_1[order.wanted[8],'Effect']) lci_rows <- c(summ_mod_1[order.wanted[5],'Lower 0.95'], summ_mod_1[order.wanted[6],'Lower 0.95'],summ_mod_2[order.wanted[6],'Lower 0.95'], summ_mod_1[order.wanted[7],'Lower 0.95'],summ_mod_1[order.wanted[8],'Lower 0.95']) uci_rows <- c(summ_mod_1[order.wanted[5],'Upper 0.95'], summ_mod_1[order.wanted[6],'Upper 0.95'],summ_mod_2[order.wanted[6],'Upper 0.95'], summ_mod_1[order.wanted[7],'Upper 0.95'],summ_mod_1[order.wanted[8],'Upper 0.95']) } hr <- c(summ_mod_1[order.wanted[1],'Effect'],summ_mod_2[order.wanted[1],'Effect'], summ_mod_3[order.wanted[1],'Effect'], summ_mod_1[order.wanted[2],'Effect'], summ_mod_1[order.wanted[3],'Effect'],summ_mod_2[order.wanted[3],'Effect'], summ_mod_1[order.wanted[4],'Effect'],summ_mod_2[order.wanted[4],'Effect'],summ_mod_4[order.wanted[4],'Effect'], summ_mod_5[order.wanted[4],'Effect'],summ_mod_6[order.wanted[4],'Effect'], hr_rows) hr <- format(round(hr,2),nsmall=2,scientific=FALSE) lci <- c(summ_mod_1[order.wanted[1],'Lower 0.95'],summ_mod_2[order.wanted[1],'Lower 0.95'], summ_mod_3[order.wanted[1],'Lower 0.95'], summ_mod_1[order.wanted[2],'Lower 0.95'], summ_mod_1[order.wanted[3],'Lower 0.95'],summ_mod_2[order.wanted[3],'Lower 0.95'], summ_mod_1[order.wanted[4],'Lower 0.95'],summ_mod_2[order.wanted[4],'Lower 0.95'],summ_mod_4[order.wanted[4],'Lower 0.95'], summ_mod_5[order.wanted[4],'Lower 0.95'],summ_mod_6[order.wanted[4],'Lower 0.95'], lci_rows) lci <- format(round(lci,2),nsmall=2,scientific=FALSE) uci <- c(summ_mod_1[order.wanted[1],'Upper 0.95'],summ_mod_2[order.wanted[1],'Upper 0.95'], summ_mod_1[order.wanted[2],'Upper 0.95'], summ_mod_1[order.wanted[3],'Upper 0.95'],summ_mod_2[order.wanted[3],'Upper 0.95'],summ_mod_3[order.wanted[3],'Upper 0.95'], summ_mod_1[order.wanted[4],'Upper 0.95'],summ_mod_2[order.wanted[4],'Upper 0.95'],summ_mod_4[order.wanted[4],'Upper 0.95'], summ_mod_5[order.wanted[4],'Upper 0.95'],summ_mod_6[order.wanted[4],'Upper 0.95'], uci_rows) uci <- format(round(uci,2),nsmall=2) ci <- paste('(',lci, ', ',uci, ')',sep='') mod_anova <- anova(mod) order.wanted.p <- c(which(rownames(mod_anova) == 'age_fhaart'), which(rownames(mod_anova) == 'male.factor'), which(rownames(mod_anova) == 'baseline_cd4_per'), which(rownames(mod_anova) == 'year_fhaart'), which(rownames(mod_anova) == 'fhaart_class_dichot.factor'), which(rownames(mod_anova) == 'log10_baseline_rna'), which(rownames(mod_anova) == 'pre_art_aids.clinical')) mod_anova_p <- mod_anova[order.wanted.p, 'P'] mod_anova_p <- ifelse(mod_anova_p < 0.001, '< 0.001', ifelse(mod_anova_p >= 0.001 & mod_anova_p < 0.01, round(mod_anova_p,3),round(mod_anova_p,2))) if(le5_ind){ age_label <- c('Baseline age (years)','~~~~~0.5','~~~~~1.0','~~~~~2.0 (ref)','~~~~~5.0') age_hr <- c(hr[1:2],'1.00',hr[3]) age_ci <- c(ci[1:2],'1.00',ci[3]) } else { age_label <- c('Baseline age (years)','~~~~~6.0','~~~~~9.0 (ref)','~~~~~12.0','~~~~~15.0') age_hr <- c(hr[1],'1.00',hr[2:3]) age_ci <- c(ci[1],'1.00',ci[2:3]) } if(vl_restrict == 0){ fhaart_label <- c('First HAART regimen class','~~~~~PI-only (ref)','~~~~~non-PI') fhaart_hr <- c('1.00',hr[12],'', '',hr[13],'1.00',hr[14],'', '','1.00',hr[15:16]) fhaart_ci <- c('',ci[12],'', '',ci[13],'',ci[14],'', '','',ci[15:16]) fhaart_times <- 3 } else { fhaart_label <- c('First HAART regimen class','~~~~~PI-only (ref)','~~~~~non-PI') fhaart_hr <- c('1.00',hr[12],'', '',hr[13],'1.00',hr[14],'', '','1.00',hr[15:16]) fhaart_ci <- c('',ci[12],'', '',ci[13],'',ci[14],'', '','',ci[15:16]) fhaart_times <- 3 } df <- data.frame('Covariate'=c(age_label,'', 'Sex','~~~~~Male','~~~~~Female (ref)','', 'Baseline CD4 percent (square root transformed)','~~~~~10','~~~~~25 (ref)','~~~~~30','', 'Year of first HAART','~~~~~2002','~~~~~2004','~~~~~2006 (ref)','~~~~~2008','~~~~~2010','~~~~~2012','', fhaart_label,'', 'Baseline HIV-1 RNA (log-transformed)','~~~~~4','~~~~~5 (ref)','~~~~~6','', 'Clinical AIDS at baseline','~~~~~not AIDS (ref)','~~~~~AIDS','~~~~~Unknown'), 'HR'=c('',age_hr,'', '',hr[4],'1.00','', '',hr[5],'1.00',hr[6],'', '',hr[7:8],'1.00',hr[9:11],'', '',fhaart_hr), '95\\% CI'=c('',age_ci,'', '',ci[4],'','', '',ci[5],'',ci[6],'', '',ci[7:8],'',ci[9:11],'', '',fhaart_ci), 'P'=c(mod_anova_p[1],rep('',times=5), mod_anova_p[2],rep('',times=3), mod_anova_p[3],rep('',times=4), mod_anova_p[4],rep('',times=7), mod_anova_p[5],rep('',times=fhaart_times), mod_anova_p[6],rep('',times=4), mod_anova_p[7],rep('',times=3)),check.names=FALSE) return(df) } # Function to do imputation for various models of interest # Used two methods to calculate p-values -- a modified likelihood ratio test and a modified Wald test. # Analysis ended up focusing on the modified Wald test results although little to no difference in the p-values # and no difference in the significance impute.fun <- function(p.cd4, r.cd4, tmp_le_5){ cd4.beta <- cd4.beta2 <- cd4.beta3 <- cd4.beta4 <- cd4.beta5 <- cd4.beta6 <- NULL cd4.se <- cd4.se.2 <- cd4.se.3 <- cd4.se.4 <- cd4.se.5 <- cd4.se.6 <- NULL new.cd4.df <- matrix(data=NA,nrow=nrow(tmp_le_5),ncol=10) vcov.full.mod <- vector("list",10) vcov.mod.wo.int <- vector("list",10) set.seed(2) # m = # loops for the imputation model m <- 5 for(i in 1:m){ cd4.predicted <- p.cd4 + sample(r.cd4,1) new.cd4 <- ifelse(is.na(tmp_le_5$baseline_cd4_sqrt),cd4.predicted,tmp_le_5$baseline_cd4_sqrt) tmp_le_5$new.cd4 <- new.cd4 new.cd4.df[,i] <- tmp_le_5$new.cd4 ddist <<- datadist(tmp_le_5) options(datadist='ddist') mod.impute <- cph(Srv_le_5 ~ rcs(age_fhaart,3) + male.factor + rcs(new.cd4,3) + rcs(year_fhaart,3) + fhaart_class.factor + pre_art_aids.clinical + strat(site_center.factor), data=tmp_le_5) # Collecting betas from various comparisons summ0 <- summary(mod.impute) row_num <- seq(1,length(rownames(summ0))) age_comp1 <- c(2.0,2.0,0.5) age_comp2 <- c(2.0, 2.0, 1.0) age_comp3 <- c(2.0, 2.0, 5.0) summ1 <- summary(mod.impute, age_fhaart=age_comp1, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(100)), year_fhaart=c(2006,2006,2002),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS'))[which(row_num %% 2 == 1),c('Effect','S.E.')] summ2 <- summary(mod.impute, age_fhaart=age_comp2, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(200)), year_fhaart=c(2006,2006,2004),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','Unknown'))[which(row_num %% 2 == 1),c('Effect','S.E.')] summ3 <- summary(mod.impute, age_fhaart=age_comp3, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(500)), year_fhaart=c(2006,2006,2006),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS'))[which(row_num %% 2 == 1),c('Effect','S.E.')] summ4 <- summary(mod.impute, age_fhaart=age_comp1, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(100)), year_fhaart=c(2006,2006,2008),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS'))[which(row_num %% 2 == 1),c('Effect','S.E.')] summ5 <- summary(mod.impute, age_fhaart=age_comp1, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(100)), year_fhaart=c(2006,2006,2010),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS'))[which(row_num %% 2 == 1),c('Effect','S.E.')] summ6 <- summary(mod.impute, age_fhaart=age_comp1, male.factor='Female',new.cd4=c(sqrt(350),sqrt(350),sqrt(100)), year_fhaart=c(2006,2006,2012),fhaart_class.factor='NNRTI-only', pre_art_aids.clinical=c('not AIDS','not AIDS','AIDS'))[which(row_num %% 2 == 1),c('Effect','S.E.')] # For 1 df Wald tests and reporting overall beta's for BMI and sex that includes interactions, # etc, collecting beta's and SE's from summary function cd4.beta <- c(cd4.beta, summ1[,'Effect']) cd4.beta2 <- c(cd4.beta2, summ2[,'Effect']) cd4.beta3 <- c(cd4.beta3, summ3[,'Effect']) cd4.beta4 <- c(cd4.beta4, summ4[,'Effect']) cd4.beta5 <- c(cd4.beta5, summ5[,'Effect']) cd4.beta6 <- c(cd4.beta6, summ6[,'Effect']) cd4.se <- c(cd4.se, summ1[,'S.E.']) cd4.se.2 <- c(cd4.se.2, summ2[,'S.E.']) cd4.se.3 <- c(cd4.se.3, summ3[,'S.E.']) cd4.se.4 <- c(cd4.se.4, summ4[,'S.E.']) cd4.se.5 <- c(cd4.se.5, summ5[,'S.E.']) cd4.se.6 <- c(cd4.se.6, summ6[,'S.E.']) } beta.var.summary <- function(mod.beta, mod.se){ beta.overall <- tapply(as.numeric(mod.beta), INDEX=names(mod.beta), FUN=mean) btwn.var <- tapply(as.numeric(mod.beta), INDEX=names(mod.beta), FUN=var) mean.var <- tapply(as.numeric(mod.se)^2, INDEX=names(mod.beta), FUN=mean) var.overall <- mean.var + 11/10*btwn.var overall.lci <- beta.overall - qnorm(0.975)*sqrt(var.overall) overall.uci <- beta.overall + qnorm(0.975)*sqrt(var.overall) wald.st <- beta.overall/sqrt(var.overall) p.val <- (1-pchisq((wald.st)^2,1)) p.val <- ifelse(p.val < 0.001, '< 0.001',ifelse(p.val >= 0.001 & p.val < 0.01, format(round(p.val,3),nsmall=3),format(round(p.val,2),nsmall=2))) # Formatting beta.overall <- format(round(exp(beta.overall),2),nsmall=2) overall.ci <- paste("(",format(round(exp(overall.lci),2),nsmall=2),", ", format(round(exp(overall.uci),2),nsmall=2),")",sep="") return(list(beta.overall=beta.overall, overall.ci=overall.ci, p.val=p.val)) } # Original model (comparing CD4 350 to 100) orig.mod.info <- beta.var.summary(mod.beta=cd4.beta, mod.se=cd4.se) # '.2' model (comparing CD4 350 to 200) sec.mod.info <- beta.var.summary(mod.beta=cd4.beta2, mod.se=cd4.se.2) # '.3' model (comparing CD4 350 to 500) third.mod.info <- beta.var.summary(mod.beta=cd4.beta3, mod.se=cd4.se.3) fourth.mod.info <- beta.var.summary(mod.beta=cd4.beta4, mod.se=cd4.se.4) fifth.mod.info <- beta.var.summary(mod.beta=cd4.beta5, mod.se=cd4.se.5) sixth.mod.info <- beta.var.summary(mod.beta=cd4.beta6, mod.se=cd4.se.6) which.order <- c(which(names(orig.mod.info[['beta.overall']]) == 'age_fhaart') , which(names(orig.mod.info[['beta.overall']]) == 'male.factor - Male:Female'), which(names(orig.mod.info[['beta.overall']]) == 'new.cd4'), which(names(orig.mod.info[['beta.overall']]) == 'year_fhaart'), which(names(orig.mod.info[['beta.overall']]) == 'fhaart_class.factor - PI-only:NNRTI-only'), which(names(orig.mod.info[['beta.overall']]) == 'fhaart_class.factor - Other:NNRTI-only'), which(names(orig.mod.info[['beta.overall']]) == 'pre_art_aids.clinical - AIDS:not AIDS'), which(names(orig.mod.info[['beta.overall']]) == 'pre_art_aids.clinical - Unknown:not AIDS')) hr <- as.numeric(as.vector(orig.mod.info[['beta.overall']][which.order])) hr.ci <- as.vector(orig.mod.info[['overall.ci']][which.order]) p.val <- as.vector(orig.mod.info[['p.val']][which.order]) names(hr) <- names(orig.mod.info[['beta.overall']])[which.order] names(hr.ci) <- names(p.val) <- names(hr) hr2 <- as.numeric(as.vector(sec.mod.info[['beta.overall']][which.order])) hr.ci2 <- as.vector(sec.mod.info[['overall.ci']][which.order]) p.val2 <- as.vector(sec.mod.info[['p.val']][which.order]) names(hr2) <- names(hr.ci2) <- names(p.val2) <- names(hr) hr3 <- as.numeric(as.vector(third.mod.info[['beta.overall']][which.order])) hr.ci3 <- as.vector(third.mod.info[['overall.ci']][which.order]) p.val3 <- as.vector(third.mod.info[['p.val']][which.order]) names(hr3) <- names(hr.ci3) <- names(p.val3) <- names(hr) hr4 <- as.numeric(as.vector(fourth.mod.info[['beta.overall']][which.order])) hr.ci4 <- as.vector(fourth.mod.info[['overall.ci']][which.order]) p.val4 <- as.vector(fourth.mod.info[['p.val']][which.order]) names(hr4) <- names(hr.ci4) <- names(p.val4) <- names(hr) hr5 <- as.numeric(as.vector(fifth.mod.info[['beta.overall']][which.order])) hr.ci5 <- as.vector(fifth.mod.info[['overall.ci']][which.order]) p.val5 <- as.vector(fifth.mod.info[['p.val']][which.order]) names(hr5) <- names(hr.ci5) <- names(p.val5) <- names(hr) hr6 <- as.numeric(as.vector(sixth.mod.info[['beta.overall']][which.order])) hr.ci6 <- as.vector(sixth.mod.info[['overall.ci']][which.order]) p.val6 <- as.vector(sixth.mod.info[['p.val']][which.order]) names(hr6) <- names(hr.ci6) <- names(p.val6) <- names(hr) overall_hr <- c('',hr[which(names(hr) == 'age_fhaart')],hr2[which(names(hr2) == 'age_fhaart')],'1.00',hr3[which(names(hr3) == 'age_fhaart')],'', '',hr[which(names(hr) == 'male.factor - Male:Female')],'1.00','', '',hr[which(names(hr) == 'new.cd4')],hr2[which(names(hr2) == 'new.cd4')],'1.00',hr3[which(names(hr3) == 'new.cd4')],'', '',hr[which(names(hr) == 'year_fhaart')],hr2[which(names(hr2) == 'year_fhaart')],'1.00',hr4[which(names(hr4) == 'year_fhaart')],hr5[which(names(hr5) == 'year_fhaart')], hr6[which(names(hr6) == 'year_fhaart')],'', '','1.00',hr[which(names(hr) == 'fhaart_class.factor - PI-only:NNRTI-only')],hr[which(names(hr) == 'fhaart_class.factor - Other:NNRTI-only')],'', '','1.00',hr[which(names(hr) == 'pre_art_aids.clinical - AIDS:not AIDS')],hr[which(names(hr) == 'pre_art_aids.clinical - Unknown:not AIDS')]) overall_hr.ci <- c('',hr.ci[which(names(hr.ci) == 'age_fhaart')],hr.ci2[which(names(hr.ci2) == 'age_fhaart')],'',hr.ci3[which(names(hr.ci3) == 'age_fhaart')],'', '',hr.ci[which(names(hr.ci) == 'male.factor - Male:Female')],'','', '',hr.ci[which(names(hr.ci) == 'new.cd4')],hr.ci2[which(names(hr.ci2) == 'new.cd4')],'',hr.ci3[which(names(hr.ci3) == 'new.cd4')],'', '',hr.ci[which(names(hr.ci) == 'year_fhaart')],hr.ci2[which(names(hr.ci2) == 'year_fhaart')],'',hr.ci4[which(names(hr.ci4) == 'year_fhaart')],hr.ci5[which(names(hr.ci5) == 'year_fhaart')], hr.ci6[which(names(hr.ci6) == 'year_fhaart')],'', '','',hr.ci[which(names(hr.ci) == 'fhaart_class.factor - PI-only:NNRTI-only')],hr.ci[which(names(hr.ci) == 'fhaart_class.factor - Other:NNRTI-only')],'', '','',hr.ci[which(names(hr.ci) == 'pre_art_aids.clinical - AIDS:not AIDS')],hr.ci[which(names(hr.ci) == 'pre_art_aids.clinical - Unknown:not AIDS')]) return(list(overall_hr=overall_hr, overall_hr.ci=overall_hr.ci)) } #CCASAnet AIDS at treatment initiation definitions #Last Modified: 29JUL2014 #Last Modified by: M. Giganti ##################### #AIDS-defining events ##################### #We are interested in generating a variable that classifies whether a patient has an AIDS-defining event prior to treatment initiation. There are four primary variables in the database for assessing a patient's clinical status. The 'aids_y' variable in the basic table, the 'aids_cl_y' variable in the basic table, the 'whostage' variable in the visit table and the 'cdcstage variable in the visit table. We want to create a composite variable using all of these sources. (05MAY2014: 'whostage','cdcstage' variables no longer in basic table) #To generate such a composite variable, the work is divided into three subtasks: (1) ADE using clinical data in visit table, (2) ADE using aids_cl_y in variable and the corresponding dates (before or after treatment initiation?) basic table, (3) ADE using the aids_y variable and the corresponding dates (before or after treatment initiation?) in basic table. #For every patient interaction in the visit table, we can classify a patient as AIDS on not AIDS using the whostage and cdcstage variables. I used the function previously written by BS. Next, we need to limit our interactions to just those which occurred prior to patient treatment initiation. Thus, we exclude all observations which occurred more than 30 days AFTER the patient's initial treatment date. Using this subsetted dataset, I next look for unique patients who had at least one visit with a ADE. Create a variable that contains a list of these patients. Next, I subset the dataset to visits that occurred no more than 365 days BEFORE the patient's treatment initiation and no more than 30 days AFTER treatment initiation. In this subset, I look for patients with at least one confirmed not AIDS visit and create a variable that contains a list of these patients. I combine these two lists -- if a patient is in both groups, I override with the AIDS classification. It is important to generate the not AIDS list because we want to differentiate between patients with not AIDS and thus with unknown status. This provides ADE using clinical data in the visit table (1). #Future consideration: Should we try to identify patients with no visits preART and their first visit after ART suggests not AIDS? ###################################### #14MAY2014 updated definition for ADE ###################################### #Previously, we used aids_y (AIDSstatus.reported) + aids_cl_y (AIDSstatus.enroll) + clinical staging from visit table (AIDSstatus.visit). This composite variable was referred to as: 'preART.aids'. #It is important to note that aids_y and aids_cl_y have different meanings. aids_cl_y designates an AIDS defining event, while aids_y is AIDS defining events + CD4 < 200. Three sites (argentina, brazil, honduras) populated the aids_y variable only. #We decided to create two variables to summarize AIDS prior to treatment initiation: #One representing clinical aids only (preARTAIDS.clinical) #Step 1: Use clinical staging from visit table (AIDSstatus.visit) #Step 2: If clinical staging from visit table is unknown, use clinical staging from basic table (AIDSstatus.enroll) #Step 3: If clinical staging from visit table is unknown and patient is from Argentina/Honduras, use aids status from basic table (AIDSstatus.reported) #Step 4: If clinical staging from visit table is not AIDS and clinical staging from basic table suggests AIDS, update as AIDS #Step 5: If clinical staging from visit table is not AIDS patient is from Argentina/Honduras and aids status from basic table suggests AIDS, update as AIDS. #Step 6: Update unknowns to not AIDS for brazil patients with aids_y=0. #One representing clinical aids + CD4 < 200 (preARTAIDS.all) ##################################################################### #Important note: The final dataset contains AIDS status for patients who initiated ART only! ###################################### #29JUL2014 update ###################################### # There is a step where we override a patient's clinical AIDS status IF: (1) they were not AIDS in the visit table, (2) they self-reported AIDS in the basic table and (3) they were from argentina/honduras. # There was a mistake in the code programming the 2nd condition and thus no changes were made. This has been updated in this iteration. #Set working directory # setwd("~/Documents/CCASAnet/Data/") #Load the plyr package - ddply() function very efficient for finding min/max observations by patient require(plyr) # #Import data # basic<-read.csv("basic.csv",header=T) # art<-read.csv("art.csv",header=T) # follow<-read.csv("follow.csv",header=T) # visit<-read.csv("visit.csv",header=T) # cd4<-read.csv("cd4.csv",header=T) # lab_rna<-read.csv("lab_rna.csv",header=T) #Create a single patient identifier that combines site and patient ID #Note: formatC is there to add leading zeros so IDs stay numeric! SitePatientID<-function(df){ df$site.patient<-paste(df$site,"-",formatC(df$patient, width = 6, format = "d", flag = "0")) } #Create unique patient ID for each data table basic$site.patient<-SitePatientID(basic) art$site.patient<-SitePatientID(art) cd4$site.patient<-SitePatientID(cd4) visit$site.patient<-SitePatientID(visit) #Date reformatting #Reformat other dates from different data tables basic$baseline_dt<-as.Date(basic$baseline_d) art$startdate<-as.Date(art$art_sd) art$enddate<-as.Date(art$art_ed) cd4$cd4_dt<-as.Date(cd4$cd4_d) visit$visit_dt<-as.Date(visit$visit_d) #Workaround to create date variable #Problem because first observation is "" temp_aids_dt<-ifelse(basic$aids_d=="","1900-01-01",as.character(basic$aids_d)) temp_aids_dt<-as.Date(temp_aids_dt) basic$aids_dt<-ifelse(temp_aids_dt=="1900-01-01",NA,temp_aids_dt) temp_aids_dt<-NULL temp_aids_cl_dt<-ifelse(basic$aids_cl_d=="","1900-01-01",as.character(basic$aids_cl_d)) temp_aids_cl_dt<-as.Date(temp_aids_cl_dt) basic$aids_cl_dt<-ifelse(temp_aids_cl_dt=="1900-01-01",NA,temp_aids_cl_dt) temp_aids_cl_dt<-NULL #Identify the ART start date for each patient Artstartdate <- ddply(art, .(site.patient), summarize,FirstARTdate=min(startdate)) #Check for lowest CD4 count before ART initiation lab_temp<-merge(cd4,Artstartdate,by="site.patient") #Sort by patient,lab date lab_temp<-with(lab_temp,lab_temp[order(site.patient,cd4_dt),]) lab_temp$ART1lag<-lab_temp$cd4_dt-lab_temp$FirstARTdate #Nadir CD4 #Does not place lower bound on CD4 date lab_temp_t<-lab_temp[lab_temp$ART1lag<=7 & is.na(lab_temp$cd4_v)==F & is.na(lab_temp$ART1lag)==F,] nadirCD4data<-ddply(lab_temp_t,.(site.patient), summarize, nadirCD4=min(cd4_v), nadirCD4_lag=max(ART1lag[cd4_v==min(cd4_v)])) ################################################### #Substep 1: ADE using clinical data in visit table ################################################### #Function to combine whostage and cdc staging info assign.clincalstage<-function(table){ clinical.stage<-with(table, ifelse(whostage==""&cdcstage=="","", ifelse(whostage==""&cdcstage!="",as.character(cdcstage), ifelse(whostage!=''&cdcstage=="",paste("WHO ",as.character(whostage),sep=""), paste(whostage,cdcstage,sep=";"))))) return(clinical.stage) } #Function to classify staging as AIDS/not AIDS classify.stage<-function(clinical.stage) { ifelse(is.na(clinical.stage)| clinical.stage=="","Unknown", ifelse(clinical.stage=="A or B"| clinical.stage=="A"| clinical.stage=="A1"| clinical.stage=="A2"| clinical.stage=="A3"| clinical.stage=="B"| clinical.stage=="B1"| clinical.stage=="B2"| clinical.stage=="B3"| clinical.stage=="WHO 1"| clinical.stage=="WHO 2"| clinical.stage=="WHO 3"| clinical.stage=="WHO 2 or 3","not AIDS", ifelse(clinical.stage=="C"| clinical.stage=="C1"| clinical.stage=="C2"| clinical.stage=="C3"| clinical.stage=="WHO 4","AIDS","SOMETHING ELSE"))) } #Run combination and classification of clinical staging for visit table visit$clinical.stage<-assign.clincalstage(visit) visit$clinical.stage.cat<-classify.stage(visit$clinical.stage) #Merge visit info with ART start date #Come up with list of patients with ART start date and an entry in visit table visit_temp<-merge(visit[,c("site.patient","site","visit_dt","clinical.stage.cat")], Artstartdate,by="site.patient",all.y=TRUE) visit_temp$indicator<-1 visit_temp$daysfromART<-visit_temp$visit_dt-visit_temp$FirstARTdate visit.uniquepatients<-unique(visit_temp[,c("site.patient","site")]) #The window between a patient's last visit with an assessment of AIDS and their treatment intiation varies by patient #A patient's AIDS status is dynamic. Currently, we are classifying a patient as non-AIDS if their last visit before treatment was non-AIDS, regardless of the length of the window between the last visit date and treatment start date #The variable - PreART.window tries to limit how recent that visit is, to distinguish between unknown and non-AIDS. #Post ART window is the window following the ART initiation date that we use to assess status at time of initiation. PreART.window<-365 PostART.window<-30 #Find unique patients with any clinical documentation of AIDS before ART visit.preARTaids<-unique( with(visit_temp, visit_temp[which(daysfromART<=PostART.window & clinical.stage.cat=="AIDS"), c("site.patient","indicator")])) colnames(visit.preARTaids)<-c("site.patient","AIDS") #Find patients with non-AIDS status at time of ART initiation #Want to identify all cases where patient has documentation of 'not-AIDS' between the pre-ART window and post-ART window visit.preARTnoaids<-unique( with(visit_temp, visit_temp[which(daysfromART>=-PreART.window & daysfromART<=PostART.window & clinical.stage.cat=="not AIDS"),c("site.patient","indicator")])) colnames(visit.preARTnoaids)<-c("site.patient","notAIDS") visit_temp3<-merge(visit.uniquepatients ,visit.preARTaids,by="site.patient", all=T) visit_temp2<-merge(visit_temp3, visit.preARTnoaids,,by="site.patient",all=T) #Classify patients as AIDS if at least one pre-ART visit with ADE, classify as non-AIDS if all pre-ART visits are non-AIDS #Assign patients with no visits as unknown AIDS status for the 'visit' variable visit_temp2$AIDSstatus.visit<-with(visit_temp2, ifelse(AIDS==1 & is.na(AIDS)==FALSE,"AIDS", ifelse(notAIDS==1 & is.na(notAIDS)==FALSE,"not AIDS", "Unknown"))) AIDSstatus.visit<-visit_temp2[,c("site.patient","AIDSstatus.visit")] ####################################################################################################################################### #Substep 2: ADE using the aids_y variable and the corresponding dates (before or after treatment initiation?) in basic table ####################################################################################################################################### ####################################################################################################################################### #Substep 3: ADE using aids_cl_y in variable and the corresponding dates (before or after treatment initiation?) basic table ####################################################################################################################################### AIDSstatus.basic<-merge(basic[,c("site.patient","patient","site","baseline_dt","aids_y","aids_dt","aids_cl_y","aids_cl_dt")], Artstartdate,all.y=TRUE) #Classify AIDS status at treatment initiation using aids_y variable #WTD: How to classify non-AIDS patients who enroll after ART initiation #WTD: How to classify non-AIDS patients who enroll more than year before ART initiation #WTD: xx patients with aids_d after ART date -Unknown or not ART? AIDSstatus.basic$AIDSstatus.reported<- with(AIDSstatus.basic, ifelse(is.na(aids_y)=="TRUE" | aids_y==9, "Unknown", ifelse(aids_y==1 & is.na(aids_dt)=="FALSE" & aids_dt <= (FirstARTdate+30), "AIDS", ifelse(aids_y==1 & is.na(aids_dt)=="FALSE" & aids_dt > (FirstARTdate+30), "Unknown", ifelse(aids_y==1 & is.na(aids_dt)=="TRUE" & baseline_dt <= (FirstARTdate+30), "AIDS", ifelse(aids_y==1 & is.na(aids_dt)=="TRUE" & baseline_dt> (FirstARTdate+30), "Unknown", ifelse(aids_y==0 & baseline_dt<=(FirstARTdate+30) & baseline_dt>=(FirstARTdate-PreART.window), "not AIDS", ifelse(aids_y==0 & baseline_dt<=(FirstARTdate+30) & baseline_dt<(FirstARTdate-PreART.window), "Unknown", ifelse(aids_y==0 & baseline_dt>(FirstARTdate-30), "Unknown", "Unknown"))))))))) #Classify AIDS status at treatment initiation using aids_cl_y variables in basic AIDSstatus.basic$AIDSstatus.enroll<- with(AIDSstatus.basic, ifelse(is.na(aids_cl_y)=="TRUE" | aids_cl_y==9, "Unknown", ifelse(aids_cl_y==1 & is.na(aids_cl_dt)=="FALSE" & aids_cl_dt <= (FirstARTdate+30), "AIDS", ifelse(aids_cl_y==1 & is.na(aids_cl_dt)=="FALSE" & aids_cl_dt > (FirstARTdate+30), "Unknown", ifelse(aids_cl_y==1 & is.na(aids_cl_dt)=="TRUE" & baseline_dt <= (FirstARTdate+30), "AIDS", ifelse(aids_cl_y==1 & is.na(aids_cl_dt)=="TRUE" & baseline_dt> (FirstARTdate+30), "Unknown", ifelse(aids_cl_y==0 & baseline_dt<=(FirstARTdate+30) & baseline_dt>=(FirstARTdate-PreART.window), "not AIDS", ifelse(aids_cl_y==0 & baseline_dt<=(FirstARTdate+30) & baseline_dt<(FirstARTdate-PreART.window), "Unknown", ifelse(aids_cl_y==0 & baseline_dt>(FirstARTdate-30), "Unknown", "Unknown"))))))))) AIDSstatus.basic<-AIDSstatus.basic[,c("site.patient","patient","site","AIDSstatus.reported","AIDSstatus.enroll")] ######################################################################################## ###Clinical AIDS definition using clinical staging from visit table + aids_cl_y + aids_y ######################################################################################## AIDSstatus_temp<-merge(AIDSstatus.basic,AIDSstatus.visit,by=c("site.patient"),all=T) AIDSstatus<-merge(AIDSstatus_temp,nadirCD4data,by=c("site.patient"),all=T) AIDSstatus$nadirCD4cat<-with(AIDSstatus,ifelse(nadirCD4<=200 & is.na(nadirCD4)==F,"LT 200", ifelse(nadirCD4>200 & abs(nadirCD4_lag)<=180 & is.na(nadirCD4)==F & is.na(nadirCD4_lag)==F,"GT 200", "Unknown"))) #Step 1: Use clinical staging from visit table (AIDSstatus.visit) #Step 2: If clinical staging from visit table is unknown, use clinical staging from basic table (AIDSstatus.enroll) #Step 3: If clinical staging from visit table is unknown and patient is from Argentina/Honduras, use aids status from basic table (AIDSstatus.reported) #Step 4: If clinical staging from visit table is not AIDS and clinical staging from basic table suggests AIDS, update as AIDS #Step 5: If clinical staging from visit table is not AIDS patient is from Argentina/Honduras and aids status from basic table #Step 6: Update unknowns to not AIDS for brazil patients with aids_y=0 #Steps 1+ 2 AIDSstatus$preARTAIDS.clinical<- with(AIDSstatus, ifelse(AIDSstatus.visit=="Unknown", AIDSstatus.enroll, AIDSstatus.visit)) #Step 3 AIDSstatus$preARTAIDS.clinical<- with(AIDSstatus, ifelse(preARTAIDS.clinical=="Unknown" & site %in% c("argentina","honduras"), AIDSstatus.reported, preARTAIDS.clinical)) #Steps 4+ 5 - Replace non-AIDS with AIDS if AIDS for any visit AIDSstatus$preARTAIDS.clinical<-with(AIDSstatus, ifelse(preARTAIDS.clinical=="not AIDS" & AIDSstatus.enroll=="AIDS", "AIDS",preARTAIDS.clinical)) #29JUL2014 - there was a mistake in the code below - fixed AIDSstatus$preARTAIDS.clinical<-with(AIDSstatus, ifelse(preARTAIDS.clinical=="not AIDS" & AIDSstatus.reported=="AIDS" & site %in% c("argentina","honduras"), "AIDS",preARTAIDS.clinical)) #update unknowns to not AIDS for brazil patients with aids_y=0. AIDSstatus$preARTAIDS.clinical<- with(AIDSstatus, ifelse(preARTAIDS.clinical=="Unknown" & site %in% c("brazil") & AIDSstatus.reported=="not AIDS", AIDSstatus.reported, preARTAIDS.clinical)) ######################################################################################## ###AIDS definition using clinical staging from visit table + aids_cl_y + aids_y + nadir CD4 < 200 ######################################################################################## #Create composite variable of aids_y (AIDSstatus.reported) + aids_cl_y (AIDSstatus.enroll) + clinical staging from visit table (AIDSstatus.visit). #Previously referred to as: 'preART.aids' #Use data from visit table as primary source #Replace unknowns with data from updated clinical status variable in basic #Replace remaining unknowns with updated aids_y variable AIDSstatus$AIDSstatus.composite<- with(AIDSstatus, ifelse(AIDSstatus.visit=="Unknown", AIDSstatus.enroll, AIDSstatus.visit)) AIDSstatus$AIDSstatus.composite<- with(AIDSstatus, ifelse(AIDSstatus.composite=="Unknown", AIDSstatus.reported, AIDSstatus.composite)) #Replace non-AIDS with AIDS if AIDS for any visit AIDSstatus$AIDSstatus.composite<-with(AIDSstatus, ifelse(AIDSstatus.composite=="not AIDS" & (AIDSstatus.reported=="AIDS" | AIDSstatus.enroll=="AIDS"), "AIDS",AIDSstatus.composite)) #Merge nadir CD4 counts for each patient AIDSstatus$preARTAIDS.all<-with(AIDSstatus, ifelse(AIDSstatus.composite=="AIDS" & nadirCD4cat=="LT 200", "AIDS", ifelse(AIDSstatus.composite=="AIDS" & nadirCD4cat=="GT 200", "AIDS", ifelse(AIDSstatus.composite=="AIDS" & nadirCD4cat=="Unknown", "AIDS", ifelse(AIDSstatus.composite=="not AIDS" & nadirCD4cat=="LT 200", "AIDS", ifelse(AIDSstatus.composite=="not AIDS" & nadirCD4cat=="GT 200", "not AIDS", ifelse(AIDSstatus.composite=="not AIDS" & nadirCD4cat=="Unknown", "not AIDS", ifelse(AIDSstatus.composite=="Unknown" & nadirCD4cat=="LT 200", "AIDS", ifelse(AIDSstatus.composite=="Unknown" & nadirCD4cat=="GT 200", "Unknown", ifelse(AIDSstatus.composite=="Unknown" & nadirCD4cat=="Unknown", "Unknown", "Not ART")))))))))) #Final N x 4 matrix containing relevant patient info AIDSstatus<-AIDSstatus[,c("site","patient","preARTAIDS.clinical","preARTAIDS.all")] save(AIDSstatus,file="/Users/giganti/Documents/CCASAnet/Latino/CCASAnetAIDSstatus.Rda")