model.fun <- function(data,cov,type){ if(cov == "go"){ var. <- "artcc_origin_g.factor" } else if (cov == "nb"){ var. <- "native.born.factor" } else if (cov == "race"){ var. <- "combined.factor" } ddist <- datadist(data) # For sub-group analyses in adjusted models data.europe <- subset(data, european.clinic.factor == "European site") data.north_american <- subset(data, european.clinic.factor == "North American site") data.europe <- upData(data.europe, drop=c("rna_u","event_d7")) ddist.europe <- datadist(data.europe) data.north_american <- upData(data.north_american, drop=c("event_d5","event_d6","event6","event7","haart0_d", "event6_name","event7_name","event8_name")) ddist.north_american <- datadist(data.north_american) #-----------------------------# # Time to any ADE # # in the first year # #-----------------------------# if(type == "all"){ Srv.ade <- Surv(data$followup_time_first_year,data$ade_first_year_of_haart) Srv.ade.europe <- Surv(data.europe$followup_time_first_year, data.europe$ade_first_year_of_haart) Srv.ade.north_american <- Surv(data.north_american$followup_time_first_year, data.north_american$ade_first_year_of_haart) } else if(type == "non_tb"){ Srv.ade <- Surv(data$followup_time_non_tb_first_year, data$non_tb_first_year_of_haart) Srv.ade.europe <- Surv(data.europe$followup_time_non_tb_first_year, data.europe$non_tb_first_year_of_haart) Srv.ade.north_american <- Surv(data.north_american$followup_time_non_tb_first_year, data.north_american$non_tb_first_year_of_haart) } # Unadjusted predictors.unadj <- c(var.,"strat(cohort)") predictors2.unadj <- c(var.,"strata(cohort)") fmla <- as.formula(paste("Srv.ade ~ ",paste(predictors.unadj,collapse="+"),sep="")) options(datadist='ddist') any.ade.unadj <- cph(fmla, data=data) # Just to get p-value for logrank test of geographic origin test <- summary(coxph(Srv.ade ~ artcc_origin_g.factor, data=data)) fmla2 <- as.formula(paste("Srv.ade ~ ",paste(predictors2.unadj,collapse="+"),sep="")) any.ade.unadj2.p <- summary(coxph(fmla2, data=data))$coefficients[,"Pr(>|z|)"] any.ade.unadj2.p <- ifelse(any.ade.unadj2.p < 0.001,"< 0.001", ifelse(any.ade.unadj2.p >= 0.001 & any.ade.unadj2.p < 0.01, round(any.ade.unadj2.p,3),round(any.ade.unadj2.p,2))) if(cov == "go"){ summ.any.ade.unadj <- summary(any.ade.unadj, artcc_origin_g.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "nb"){ summ.any.ade.unadj <- summary(any.ade.unadj, native.born.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "race"){ summ.any.ade.unadj <- summary(any.ade.unadj, combined.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } any.ade.unadj.df <- round(summ.any.ade.unadj[rownames(summ.any.ade.unadj)==" Hazard Ratio", c("Effect","Lower 0.95","Upper 0.95")],2) if(is.vector(any.ade.unadj.df)){ any.ade.unadj.df <- t(as.matrix(any.ade.unadj.df)) } ci <- paste("(",any.ade.unadj.df[,"Lower 0.95"],", ",any.ade.unadj.df[,"Upper 0.95"],")",sep="") any.ade.unadj.df <- cbind(any.ade.unadj.df[,1],ci) any.ade.unadj.p <- anova(any.ade.unadj)[1,"P"] any.ade.unadj.p <- ifelse(any.ade.unadj.p < 0.001, "< 0.001", ifelse(any.ade.unadj.p < 0.01 & any.ade.unadj.p >= 0.001,round(any.ade.unadj.p,3), round(any.ade.unadj.p,2))) if(cov == "go"){ covariate <- c("Northern Africa and Middle East","Sub-Saharan Africa","Asia","Latin America") column. <- "Geographic Origin (ref=West)" } else if (cov == "nb"){ covariate <- "Non-native born" column. <- "Native born status" } else if (cov == "race"){ covariate <- c("Black","Hispanic","Asian","Arab","Indigenous") column. <- "Ethnicity/Race (ref=White)" } any.ade.unadj.df <- cbind("Covariate"=covariate,any.ade.unadj.df,any.ade.unadj2.p) colnames(any.ade.unadj.df) <- c(column.,"Hazard Ratio","95\\% CI","P") # Adjusted predictors.adj <- c(var.,"gender.factor","rcs(cd4_t0_corrected,4)","rcs(lrna_t0_corrected,3)", "rcs(age,3)","pi_based.factor","nnrti_based.factor", "rcs(haartyr,3)", "risk.factor","stagec.factor","strat(cohort)") predictors2.adj <- c(var.,"gender.factor","rcs(cd4_t0_corrected,4)","rcs(lrna_t0_corrected,3)", "rcs(age,3)","pi_based.factor","nnrti_based.factor", "rcs(haartyr,3)", "risk.factor","stagec.factor","strata(cohort)") fmla <- as.formula(paste("Srv.ade ~ ",paste(predictors.adj,collapse="+"))) fmla.europe <- as.formula(paste("Srv.ade.europe ~ ",paste(predictors.adj,collapse="+"))) fmla.north_american <- as.formula(paste("Srv.ade.north_american ~ ",paste(predictors.adj,collapse="+"))) options(datadist='ddist') any.ade.adj <- cph(fmla, data=data) summ.any.ade.adj <- summary.fun(mod=any.ade.adj,cov=cov,risk.levels=levels(data$risk.factor)) summ.any.ade.adj.format <- format.fun(summ.=summ.any.ade.adj, mod=any.ade.adj,cov=cov, risk.levels=levels(data$risk.factor), main.levels=levels(eval(parse(text=paste("data$",var.,sep=""))))) options(datadist='ddist.europe') any.ade.adj.europe <- cph(fmla.europe, data=data.europe) summ.any.ade.adj.europe <- summary.fun(mod=any.ade.adj.europe,cov=cov,risk.levels=levels(data.europe$risk.factor)) summ.any.ade.adj.europe.format <- format.fun(summ.=summ.any.ade.adj.europe, mod=any.ade.adj.europe,cov=cov, risk.levels=levels(data.europe$risk.factor), main.levels=levels(eval(parse(text=paste("data.europe$",var.,sep=""))))) options(datadist='ddist.north_american') any.ade.adj.north_american <- cph(fmla.north_american, data=data.north_american) summ.any.ade.adj.north_american <- summary.fun(mod=any.ade.adj.north_american,cov=cov, risk.levels=levels(data.north_american$risk.factor)) summ.any.ade.adj.north_american.format <- format.fun(summ.=summ.any.ade.adj.north_american, mod=any.ade.adj.north_american,cov=cov, risk.levels=levels(data.north_american$risk.factor), main.levels=levels(eval(parse(text=paste("data.north_american$",var.,sep=""))))) fmla.int <- as.formula(paste('Srv.ade ~ ',paste(c('native.born.factor*european.clinic.factor', 'strata(cohort)'),collapse='+'))) #-----------------------------# # Time to ADE or death # # in the first year # #-----------------------------# if(type == "all"){ Srv.ade.dth <- Surv(data$followup_time_first_year,data$ade_or_death_first_year_of_haart) Srv.ade.dth.europe <- Surv(data.europe$followup_time_first_year,data.europe$ade_or_death_first_year_of_haart) Srv.ade.dth.north_american <- Surv(data.north_american$followup_time_first_year, data.north_american$ade_or_death_first_year_of_haart) } else if(type == "non_tb"){ Srv.ade.dth <- Surv(data$followup_time_non_tb_first_year, data$non_tb_or_death_first_year_of_haart) Srv.ade.dth.europe <- Surv(data.europe$followup_time_non_tb_first_year, data.europe$non_tb_or_death_first_year_of_haart) Srv.ade.dth.north_american <- Surv(data.north_american$followup_time_non_tb_first_year, data.north_american$non_tb_or_death_first_year_of_haart) } # Unadjusted options(datadist='ddist') fmla <- as.formula(paste("Srv.ade.dth ~ ",paste(predictors.unadj,collapse="+"),sep="")) ade.or.dth.unadj <- cph(fmla, data=data) fmla2 <- as.formula(paste("Srv.ade.dth ~ ",paste(predictors2.unadj,collapse="+"),sep="")) ade.or.dth.unadj2.p <- summary(coxph(fmla2, data=data))$coefficients[,"Pr(>|z|)"] ade.or.dth.unadj2.p <- ifelse(ade.or.dth.unadj2.p < 0.001,"< 0.001", ifelse(ade.or.dth.unadj2.p >= 0.001 & ade.or.dth.unadj2.p < 0.01, round(ade.or.dth.unadj2.p,3),round(ade.or.dth.unadj2.p,2))) if(cov == "go"){ summ.ade.or.dth.unadj <- summary(ade.or.dth.unadj, artcc_origin_g.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "nb"){ summ.ade.or.dth.unadj <- summary(ade.or.dth.unadj, native.born.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "race"){ summ.ade.or.dth.unadj <- summary(ade.or.dth.unadj, combined.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } ade.or.dth.unadj.df <- round(summ.ade.or.dth.unadj[rownames(summ.any.ade.unadj)==" Hazard Ratio", c("Effect","Lower 0.95","Upper 0.95")],2) if(is.vector(ade.or.dth.unadj.df)){ ade.or.dth.unadj.df <- t(as.matrix(ade.or.dth.unadj.df)) } ci <- paste("(",ade.or.dth.unadj.df[,"Lower 0.95"],", ",ade.or.dth.unadj.df[,"Upper 0.95"],")",sep="") ade.or.dth.unadj.df <- cbind(ade.or.dth.unadj.df[,1],ci) ade.or.dth.unadj.p <- anova(ade.or.dth.unadj)[1,"P"] ade.or.dth.unadj.p <- ifelse(ade.or.dth.unadj.p < 0.001, "< 0.001", ifelse(ade.or.dth.unadj.p < 0.01 & ade.or.dth.unadj.p >= 0.001,round(ade.or.dth.unadj.p,3), round(ade.or.dth.unadj.p,2))) ade.or.dth.unadj.df <- cbind("Covariate"=covariate,ade.or.dth.unadj.df,ade.or.dth.unadj2.p) colnames(ade.or.dth.unadj.df) <- c(column.,"Hazard Ratio","95\\% CI","P") # Adjusted fmla <- as.formula(paste("Srv.ade.dth ~ ",paste(predictors.adj, collapse="+"))) fmla.europe <- as.formula(paste("Srv.ade.dth.europe ~ ",paste(predictors.adj, collapse="+"))) fmla.north_american <- as.formula(paste("Srv.ade.dth.north_american ~ ",paste(predictors.adj, collapse="+"))) options(datadist='ddist') ade.or.dth.adj <- cph(fmla, data=data) summ.ade.or.dth.adj <- summary.fun(mod=ade.or.dth.adj,cov=cov,risk.levels=levels(data$risk.factor)) summ.ade.or.dth.format <- format.fun(summ.=summ.ade.or.dth.adj, mod=ade.or.dth.adj, cov=cov, risk.levels=levels(data$risk.factor), main.levels=levels(eval(parse(text=paste("data$",var.,sep=""))))) summ.ade.or.dth.format[[1]] <- summ.ade.or.dth.format[[1]][,-1] options(datadist='ddist.europe') ade.or.dth.adj.europe <- cph(fmla.europe, data=data.europe) summ.ade.or.dth.adj.europe <- summary.fun(mod=ade.or.dth.adj.europe,cov=cov, risk.levels=levels(data.europe$risk.factor)) summ.ade.or.dth.europe.format <- format.fun(summ.=summ.ade.or.dth.adj.europe, mod=ade.or.dth.adj.europe, cov=cov, risk.levels=levels(data.europe$risk.factor), main.levels=levels(eval(parse(text=paste("data.europe$",var.,sep=""))))) summ.ade.or.dth.europe.format[[1]] <- summ.ade.or.dth.europe.format[[1]][,-1] options(datadist='ddist.north_american') ade.or.dth.adj.north_american <- cph(fmla.north_american, data=data.north_american) summ.ade.or.dth.adj.north_american <- summary.fun(mod=ade.or.dth.adj.north_american,cov=cov, risk.levels=levels(data.north_american$risk.factor)) summ.ade.or.dth.north_american.format <- format.fun(summ.=summ.ade.or.dth.adj.north_american, mod=ade.or.dth.adj.north_american, cov=cov, risk.levels=levels(data.north_american$risk.factor), main.levels=levels(eval(parse(text=paste("data.north_american$",var.,sep=""))))) summ.ade.or.dth.north_american.format[[1]] <- summ.ade.or.dth.north_american.format[[1]][,-1] #-----------------------------# # Time to L2FU # # in the first year # #-----------------------------# Srv.l2fu <- Surv(as.numeric(data$followup_time_l2fu_first_year),data$l2fu) # Unadjusted fmla <- as.formula(paste("Srv.l2fu ~ ",paste(predictors.unadj,collapse="+"),sep="")) options(datadist='ddist') l2fu.unadj <- cph(fmla, data=data) fmla2 <- as.formula(paste("Srv.l2fu ~ ",paste(predictors2.unadj,collapse="+"),sep="")) l2fu.unadj2.p <- summary(coxph(fmla2,data=data))$coefficients[,"Pr(>|z|)"] l2fu.unadj2.p <- ifelse(l2fu.unadj2.p < 0.001,"< 0.001", ifelse(l2fu.unadj2.p >= 0.001 & l2fu.unadj2.p < 0.01, round(l2fu.unadj2.p,3),round(l2fu.unadj2.p,2))) if(cov == "go"){ summ.l2fu.unadj <- summary(l2fu.unadj, artcc_origin_g.factor = levels(eval(parse(text=paste("data$",var.,sep=""))))[1]) } else if (cov == "nb"){ summ.l2fu.unadj <- summary(l2fu.unadj, native.born.factor = levels(eval(parse(text=paste("data$",var.,sep=""))))[1]) } else if (cov == "race"){ summ.l2fu.unadj <- summary(l2fu.unadj, combined.factor = levels(eval(parse(text=paste("data$",var.,sep=""))))[1]) } l2fu.unadj.df <- round(summ.l2fu.unadj[rownames(summ.l2fu.unadj)==" Hazard Ratio",c("Effect","Lower 0.95","Upper 0.95")],2) if(is.vector(l2fu.unadj.df)){ l2fu.unadj.df <- t(as.matrix(l2fu.unadj.df)) } ci <- paste("(",l2fu.unadj.df[,"Lower 0.95"],", ",l2fu.unadj.df[,"Upper 0.95"],")",sep="") l2fu.unadj.df <- cbind(l2fu.unadj.df[,1],ci) l2fu.unadj.p <- anova(l2fu.unadj)[1,"P"] l2fu.unadj.p <- ifelse(l2fu.unadj.p < 0.001, "< 0.001", ifelse(l2fu.unadj.p >= 0.001 & l2fu.unadj.p < 0.01,round(l2fu.unadj.p,3),round(l2fu.unadj.p,2))) l2fu.unadj.df <- cbind("Covariate"=covariate,l2fu.unadj.df,l2fu.unadj2.p) colnames(l2fu.unadj.df) <- c(column.,"Hazard Ratio","95\\% CI","P") # Adjusted fmla <- as.formula(paste("Srv.l2fu ~ ",paste(predictors.adj,collapse="+"))) l2fu.adj <- cph(fmla, data=data) summ.l2fu.adj <- summary.fun(mod=l2fu.adj, cov=cov, risk.levels=levels(data$risk.factor)) summ.l2fu.adj.format <- format.fun(summ.=summ.l2fu.adj, mod=l2fu.adj,cov=cov, risk.levels=levels(data$risk.factor), main.levels=levels(eval(parse(text=paste("data$",var.,sep=""))))) #-----------------------------# # Time to death # # in the first year # #-----------------------------# if(type == "all"){ Srv.dth <- Surv(data$followup_time_first_year,data$death_first_year_of_haart) } else if(type == "non_tb"){ Srv.dth <- Surv(data$followup_time_non_tb_first_year, data$death_first_year_of_haart) } # Unadjusted options(datadist='ddist') fmla <- as.formula(paste("Srv.dth ~ ",paste(predictors.unadj,collapse="+"),sep="")) dth.unadj <- cph(fmla, data=data) fmla2 <- as.formula(paste("Srv.dth ~ ",paste(predictors2.unadj,collapse="+"),sep="")) dth.unadj2.p <- summary(coxph(fmla2, data=data))$coefficients[,"Pr(>|z|)"] dth.unadj2.p <- ifelse(dth.unadj2.p < 0.001,"< 0.001", ifelse(dth.unadj2.p >= 0.001 & dth.unadj2.p < 0.01, round(dth.unadj2.p,3),round(dth.unadj2.p,2))) if(cov == "go"){ summ.dth.unadj <- summary(dth.unadj, artcc_origin_g.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "nb"){ summ.dth.unadj <- summary(dth.unadj, native.born.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } else if (cov == "race"){ summ.dth.unadj <- summary(dth.unadj, combined.factor = levels(eval(parse(text=paste("data$", var.,sep=""))))[1]) } dth.unadj.df <- round(summ.dth.unadj[rownames(summ.dth.unadj)==" Hazard Ratio", c("Effect","Lower 0.95","Upper 0.95")],2) if(is.vector(dth.unadj.df)){ dth.unadj.df <- t(as.matrix(dth.unadj.df)) } ci <- paste("(",dth.unadj.df[,"Lower 0.95"],", ",dth.unadj.df[,"Upper 0.95"],")",sep="") dth.unadj.df <- cbind(dth.unadj.df[,1],ci) dth.unadj.p <- anova(dth.unadj)[1,"P"] dth.unadj.p <- ifelse(dth.unadj.p < 0.001, "< 0.001", ifelse(dth.unadj.p < 0.01 & dth.unadj.p >= 0.001,round(dth.unadj.p,3), round(dth.unadj.p,2))) dth.unadj.df <- cbind("Covariate"=covariate,dth.unadj.df,dth.unadj2.p) colnames(dth.unadj.df) <- c(column.,"Hazard Ratio","95\\% CI","P") # Adjusted fmla <- as.formula(paste("Srv.dth ~ ",paste(predictors.adj, collapse="+"))) options(datadist='ddist') dth.adj <- cph(fmla, data=data) # Just to get p-value for logrank test of geographic origin test <- summary(coxph(Srv.dth ~ artcc_origin_g.factor, data=data)) summ.dth.adj <- summary.fun(mod=dth.adj,cov=cov,risk.levels=levels(data$risk.factor)) summ.dth.format <- format.fun(summ.=summ.dth.adj, mod=dth.adj, cov=cov, risk.levels=levels(data$risk.factor), main.levels=levels(eval(parse(text=paste("data$",var.,sep=""))))) summ.dth.format[[1]] <- summ.dth.format[[1]][,-1] return(list(any.ade.unadj.df=any.ade.unadj.df, ade.or.dth.unadj.df=ade.or.dth.unadj.df, any.ade.adj.df=summ.any.ade.adj.format[[1]], ade.or.dth.adj.df=summ.ade.or.dth.format[[1]], any.ade.unadj.p=any.ade.unadj.p, ade.or.dth.unadj.p=ade.or.dth.unadj.p, any.ade.adj.nonlinear.p=summ.any.ade.adj.format[[2]], ade.or.dth.nonlinear.p=summ.ade.or.dth.format[[2]], any.ade.adj=any.ade.adj, ade.or.dth.adj=ade.or.dth.adj, l2fu.unadj.df=l2fu.unadj.df, l2fu.adj.df = summ.l2fu.adj.format[[1]], l2fu.unadj.p=l2fu.unadj.p, l2fu.adj.nonlinear.p=summ.l2fu.adj.format[[2]], l2fu.adj=l2fu.adj, any.ade.adj.europe=any.ade.adj.europe, ade.or.dth.adj.europe=ade.or.dth.adj.europe, any.ade.adj.north_american=any.ade.adj.north_american, ade.or.dth.adj.north_american=ade.or.dth.adj.north_american, any.ade.adj.europe.df=summ.any.ade.adj.europe.format[[1]], ade.or.dth.adj.europe.df=summ.ade.or.dth.europe.format[[1]], any.ade.adj.north_american.df=summ.any.ade.adj.north_american.format[[1]], ade.or.dth.adj.north_american.df=summ.ade.or.dth.north_american.format[[1]], dth.adj.df=summ.dth.format[[1]], dth.adj.nonlinear.p=summ.dth.format[[2]], dth.adj=dth.adj, dth.unadj.df=dth.unadj.df)) }