################################################################## ####### This code performs the dynamic marginal structural model analysis rm(list=ls()) setwd("/Volumes/encrypted/cnics/analyses-2013-Nov") #------------------------# # Functions # #------------------------# #---------------------------------------------------------------------# # add.fun # #---------------------------------------------------------------------# # In creating monthly (time-dependent) data frame, records are inserted # that do not exist in the original data. This creates missing values for # CD4, viral loads, CD4 %, etc. Add.fun fills in these values with the # last non-missing observation and records the age at which that non-missing # value occurred. add.fun <- function(labage.1, labage.2){ var <- names(labage.1)[3] labage.1$order <- 1 labage.2$order <- 0 pullData <- rbind(labage.1,labage.2) pullData$age <- round(pullData$age, 4) pullData <- pullData[order(pullData$pid, pullData$age, pullData$order),] hasVal <- diff(c(0L,as.integer(as.factor(pullData$pid)))) != 0 | !is.na(pullData[,var]) | (pullData$order == 0) index <- rep(which(hasVal), times=table(cumsum(hasVal))) pullData[,c(var, paste("age", var, sep='.'))] <- pullData[index, c(var, "age")] return(pullData) } u.limit=401 l.limit=301 start.within=180/365.25 LOST=FALSE COVARIATES=FALSE SITE.YEAR=FALSE EXCLUDE.IDU=FALSE DEATH.ONLY=FALSE CENSOR.TIME=FALSE EXCLUDE.DEATHS=FALSE ADD2=FALSE EXCLUDE.NON.HAART.ART=TRUE TRUNC.LEVEL=0.99 GAP1YEAR.LTFU=TRUE setwd("/Volumes/encrypted/cnics/analyses-2013-Nov") load("cnics-hernan-stuff5c.Rda") load("ltfu_gap.Rda") ####### Analysis function that does pretty much everything hernan.analysis<-function(d2, u.limit, l.limit, filename, start.within=180/365.25, LOST=FALSE, COVARIATES=FALSE, SITE.YEAR=FALSE, EXCLUDE.IDU=FALSE, DEATH.ONLY=FALSE, CENSOR.TIME=FALSE, EXCLUDE.DEATHS=FALSE, ADD2=FALSE, EXCLUDE.NON.HAART.ART=TRUE, TRUNC.LEVEL=0.99, GAP1YEAR.LTFU=FALSE) { if (EXCLUDE.NON.HAART.ART) { ### Currently, this is the only option. At some point, I should censor but I'm not sure VU gave me this data. d2<-d2[d2$non.haart.art==0,] d2$non.haart.art<-NULL } if (ADD2==TRUE){ d2$last.age<-with(d2, ifelse(deceased==0, last.age+60/365.25, last.age)) } # Identifying age at which CD4 is less than ulimit by PID less.ulimit <- lapply(split(subset(d2, select=c(pid,age,cd4)), d2$pid), FUN=function(y){ y[min(which(y[,"cd4"] < u.limit),na.rm=TRUE),"age"]}) less.ulimit.df <- data.frame("pid"=names(less.ulimit), "age.less.ulimit"=unlist(less.ulimit)) less.ulimit.df$pid <- as.character(less.ulimit.df$pid) d2[,"age.less.ulimit"] <- less.ulimit.df[match(d2$pid,less.ulimit.df$pid,nomatch=NA), "age.less.ulimit"] d2 <- d2[order(d2$pid, d2$age),] # If a PID has a missing age.less.ulimit, then they did not have a CD4 below the upper # limit and should be deleted from the analysis d2 <- subset(d2, !is.na(age.less.ulimit)) # Time.in.care is the time between the first recorded clinic visit and the age a PIDs CD4 fell # below the upper limit d2$time.in.care <- d2$age.less.ulimit - d2$first.age ##### adding by Qi for ltfu defined as >1 year gap if(GAP1YEAR.LTFU){ d.ltfu.gap <- merge(d2, ltfu.gap.all, by.x="pid",by.y="id", all.x=TRUE) d.ltfu.gap$ltfu.age <- round(d.ltfu.gap$ltfu.age,3) tmp <- subset(d.ltfu.gap, (ltfu.age>=age.less.ulimit)) #### only need the first >1 year gap for each patient tmp<- tmp[order(tmp$pid, tmp$age, tmp$ltfu.age), ] tmp <- tmp[!duplicated(tmp[, c("pid", "age")]), ] tmp$last.age.1 <- pmin(tmp$last.age, tmp$ltfu.age) tmp <- tmp[!duplicated(tmp$pid),] tmp$last.age.1 <- pmin(tmp$last.age, tmp$ltfu.age) d3 <- merge(d2, tmp[, c("pid", "ltfu.age", "last.age.1")], by="pid", all.x=TRUE) d3$last.age.1 <- ifelse(is.na(d3$last.age.1), d3$last.age, d3$last.age.1) d3 <- subset(d3, age<=last.age.1) d3 <- d3[order(d3$pid, d3$age),] d3$deceased1 <- ifelse(is.na(d3$ltfu.age), d3$deceased, ifelse(round(d3$last.age.1, 3)=0,] #Doesn't appear to be needed (2014-03-07) # In the original data, each record at which an OI occurred had a non-missing age.oi. All other records # had a missing age.oi. For analysis purposes, those OIs that occurred between a PIDs first age and # age.less.ulimit were recorded as prior.oi and the age recorded. For those OIs that occurred after # age.less.ulimit, only the first occurrence of an OI will be included in the analysis and the variable # age.oi will repeat the age at which the first OI after age.less.ulimit occurred. The following # code addresses both of these issues. # Defining prior.oi (OI between first.age and age.less.ulimit) age.oi <- subset(d2, !is.na(age.oi), select=c(pid, age.oi, age.less.ulimit)) age.oi$prior.oi <- ifelse(round(age.oi$age.oi,4) <= round(age.oi$age.less.ulimit,4),1,0) ##### Added by Bryan to try to simplify prior.oi.df<-data.frame(pid=unique(age.oi$pid[age.oi$prior.oi==1]),prior.oi=1) d2<-merge(d2,prior.oi.df,by="pid",all.x=TRUE) d2 <- d2[order(d2$pid,d2$age),] d2$prior.oi <- ifelse(is.na(d2$prior.oi),0,d2$prior.oi) # Grabbing first incidence of OI after age.less.ulimit if more than one oi age.oi.af <- lapply(split(age.oi, age.oi$pid), FUN=function(y){ ifelse(all(y[,"prior.oi"]==1),NA,min(y[y[,"prior.oi"]==0,"age.oi"]))}) age.oi.af.df <- data.frame("pid"=unlist(names(age.oi.af)), "age.oi.af"=unlist(age.oi.af)) age.oi.af.df$pid <- as.character(age.oi.af.df$pid) age.oi.af.df<-age.oi.af.df[!is.na(age.oi.af.df$age.oi.af),] d2 <- merge(d2, age.oi.af.df, by="pid", all.x=TRUE) d2 <- d2[order(d2$pid,d2$age),] d2$age.oi<-d2$age.oi.af d2$age.oi.af<-NULL # Subsetting out PIDs with prior.oi d2 <- subset(d2, prior.oi != 1) unique.id<-unique(d2$pid) gt.llimit.df <- subset(d2, round(age,4) == round(age.less.ulimit,4), select=c(pid,age,age.less.ulimit,cd4)) gt.llimit.df$gt.llimit <- ifelse(gt.llimit.df$cd4 >= l.limit,1,0) d2[,"gt.llimit"] <- gt.llimit.df[match(d2$pid, gt.llimit.df$pid,nomatch=NA), "gt.llimit"] sub.d2 <- subset(d2, (age.fhaart >= age.less.ulimit | is.na(age.fhaart)) & gt.llimit == 1) sub.d2 <- sub.d2[order(sub.d2$pid,sub.d2$age),] sub.d2$gt.llimit<-NULL sub.unique.id <- unique(sub.d2$pid) d2 <- sub.d2 ####### I'm being bold, and changing the structure of the code here. I'm going to expand to months first and then I'll define ####### treatment groups. #### I'm going to keep all of the repeating covariates to avoid changing them later, but I'll grab baseline values tmp<-d2[!duplicated(d2$pid),] tmp$age.rx.start<-tmp$haart_before_first_visit<-NULL tmp$age<-tmp$oi<-tmp$last.age<-tmp$prior.oi<-tmp$haart<-tmp$fhaart<-tmp$age.approx<-NULL max.follow <- 12*max(tmp$tot.fu, na.rm=TRUE) months <- seq(0, ceiling(max.follow),by=1) tmp$month <- 0 tmp2<-data.frame(pid=rep(tmp$pid,each=length(months[-1])), site=rep(tmp$site,each=length(months[-1])), black=rep(tmp$black,each=length(months[-1])), white=rep(tmp$white,each=length(months[-1])), deceased=rep(tmp$deceased,each=length(months[-1])), age_at_death=rep(tmp$age_at_death,each=length(months[-1])), idu=rep(tmp$idu,each=length(months[-1])), age.oi=rep(tmp$age.oi,each=length(months[-1])), cd4=NA, male=rep(tmp$male,each=length(months[-1])), first.age=rep(tmp$first.age,each=length(months[-1])), first.year=rep(tmp$first.year,each=length(months[-1])), ltfu=rep(tmp$ltfu,each=length(months[-1])), age.fhaart=rep(tmp$age.fhaart,each=length(months[-1])), age.less.ulimit=rep(tmp$age.less.ulimit,each=length(months[-1])), time.in.care=rep(tmp$time.in.care,each=length(months[-1])), true.last.age=rep(tmp$true.last.age,each=length(months[-1])), tot.fu=rep(tmp$tot.fu,each=length(months[-1])), month=rep(months[-1],times=length(tmp$pid))) surv.df<-rbind(tmp,tmp2) surv.df <- surv.df[order(surv.df$pid, surv.df$month),] surv.df$age <- round(surv.df$age.less.ulimit + surv.df$month/12,4) ###### Defining month of first HAART using the larger dataframe locate.fhaart <- lapply(split(subset(surv.df, select=c(pid,age,age.fhaart)),surv.df$pid), FUN=function(y){ y[min(which(y[,"age"] >= y[,"age.fhaart"])),"age"] }) locate.fhaart.df <- data.frame("pid"=names(locate.fhaart)[!is.na(locate.fhaart)], "new.age.fhaart"=unlist(locate.fhaart[!is.na(locate.fhaart)])) surv.df[,"new.age.fhaart"] <- locate.fhaart.df[match(surv.df$pid,locate.fhaart.df$pid,nomatch=NA), "new.age.fhaart"] surv.df$fhaart <- ifelse(round(surv.df$age,4) == round(surv.df$new.age.fhaart,4) & !is.na(surv.df$new.age.fhaart),1, ifelse(is.na(surv.df$new.age.fhaart) | round(surv.df$age,4) != round(surv.df$new.age.fhaart,4),0,NA)) locate.oi <- lapply(split(subset(surv.df, select=c(pid,age,age.oi)),surv.df$pid), FUN=function(y){ y[min(which(y[,"age"] >= y[,"age.oi"])),"age"] }) locate.oi.df <- data.frame("pid"=names(locate.oi)[!is.na(locate.oi)], "new.age.oi"=unlist(locate.oi[!is.na(locate.oi)])) surv.df[,"new.age.oi"] <- locate.oi.df[match(surv.df$pid,locate.oi.df$pid,nomatch=NA), "new.age.oi"] surv.df$oi <- ifelse(round(surv.df$age,4) == round(surv.df$new.age.oi,4) & !is.na(surv.df$new.age.oi),1, ifelse(is.na(surv.df$new.age.oi) | round(surv.df$age,4) != round(surv.df$new.age.oi,4),0,NA)) ###### Getting rid of time after follow-up. I did this here, because some of the above code needed it. # # #surv.df <- subset(surv.df, round(age,4) <=round(true.last.age,4)) ### Commented out on 2014-03-07 # # # # # # surv.df$death <- ifelse(is.na(surv.df$age_at_death),0, # # ifelse((round(surv.df$age,4) == # # round(surv.df$true.last.age,4)) | # # (surv.df$age + 1/12 > surv.df$true.last.age),1,0)) # # # # # # surv.df$oi <- ifelse(is.na(surv.df$age.oi),0, # # ifelse((surv.df$age + 1/12) > surv.df$true.last.age & # # (surv.df$age.oi + 1/12) > surv.df$true.last.age,1, # surv.df$oi)) surv.df <- subset(surv.df, round(age,4) < round(true.last.age,4) + 1/12) ### Changed on 2014-03-07 surv.df$death <- ifelse(is.na(surv.df$age_at_death),0, ifelse((round(surv.df$age,4) == round(surv.df$true.last.age,4)) | (surv.df$age > surv.df$true.last.age),1,0)) surv.df$oi <- ifelse(is.na(surv.df$age.oi),0, ifelse((surv.df$age ) > surv.df$true.last.age & (surv.df$age.oi) > surv.df$true.last.age,1, surv.df$oi)) surv.df$month<-ifelse(surv.df$month==0&surv.df$death==1,1,surv.df$month) if (DEATH.ONLY) { surv.df$event <- ifelse(surv.df$death == 0, 0,1) } if (!DEATH.ONLY) { surv.df$event <- ifelse(surv.df$death == 0 & surv.df$oi == 0,0,1) } ###### Adding CD4 counts ref.cd4 <- data.frame("pid"=d2$pid[!duplicated(d2$pid)], # "ref.cd4"=d2$cd4[round(d2$age,4)==round(d2$age.less.ulimit,4)], #### Bryan changed because I'll use CD4 at first visit as covariate (This was probably a poor choice, I fix back later - 2014-03-10) "ref.cd4"=d2$cd4[!duplicated(d2$pid)], "age.ref.cd4"=d2$age.less.ulimit[!duplicated(d2$pid)], "age.fhaart"=d2$age.fhaart[!duplicated(d2$pid)]) surv.df <- merge(surv.df, ref.cd4, all=TRUE) surv.df <- surv.df[order(surv.df$pid,surv.df$age),] labage.1 <- surv.df[,c("pid","age","cd4")] labage.1$cd4<-NA labage.2 <- subset(d2, pid %in% surv.df$pid, select=c(pid,age,cd4)) names(labage.2) <- c("pid","age","cd4") labage.2 <- subset(labage.2, !is.na(cd4)) pullData <- add.fun(labage.1, labage.2) surv.df$cd4 <- pullData[pullData$order==1,"cd4"] surv.df$age <- round(surv.df$age, 4) surv.df$age.fhaart <- round(surv.df$age.fhaart,4) surv.df$age.cd4 <- pullData[pullData$order==1, "age.cd4"] surv.df$time.to.cd4 <- 12*(surv.df$age - surv.df$age.cd4) junk<-subset(surv.df[!duplicated(surv.df$pid),],select=c(pid,cd4)) ###### The wrong ref.cd4 was saved above. I'm fixing it here to be the CD4 when first going below upper limit. 2014-03-10 surv.df<-merge(surv.df,junk,all.x=TRUE) surv.df<-surv.df[order(surv.df$pid,surv.df$age),] ## Need to ensure that for cd4 that age.cd4 ## is less than age.fhaart when fhaart==1 # Adjusting those who started HAART and had a CD4 count in that month after # the age.fhaart which.cd4 <- which(round(surv.df$age.fhaart,4) < round(surv.df$age.cd4,4) & surv.df$fhaart == 1 & round(surv.df$age.cd4,4) > round(surv.df$age.less.ulimit,4)) surv.df$cd4[which.cd4] <- surv.df$cd4[(which.cd4-1)] surv.df$age.cd4[which.cd4] <- surv.df$age.cd4[(which.cd4-1)] surv.df$time.to.cd4[which.cd4] <- 12*(surv.df$age[which.cd4]-surv.df$age.cd4[which.cd4]) ### Line added 12/12/14 (should have minor impact) #### do the same thing for oi: adjusting those who had oi and CD4 count in the same month which.cd4.oi <- which(round(surv.df$age.oi,4) < round(surv.df$age.cd4,4) & surv.df$oi== 1 & round(surv.df$age.oi,4) > round(surv.df$age.less.ulimit,4)) surv.df$cd4[which.cd4.oi] <- surv.df$cd4[(which.cd4.oi-1)] surv.df$age.cd4[which.cd4.oi] <- surv.df$age.cd4[(which.cd4.oi-1)] surv.df$time.to.cd4[which.cd4.oi] <- 12*(surv.df$age[which.cd4.oi]-surv.df$age.cd4[which.cd4.oi]) ### Line added 12/12/14 (should have minor impact) if (!DEATH.ONLY){ surv.df<-surv.df[is.na(surv.df$new.age.oi) | surv.df$age<=surv.df$new.age.oi,] } surv.df$haart<-with(surv.df, ifelse(is.na(new.age.fhaart),0, ifelse(age=1) { # locate.cd4.threshold[i]<-which(surv.df$age.cd4==min(surv.df$age.cd4[surv.df$cd4start.within) | (haart==1 & age.fhaart-age.less.ulimitage.cd4.llimit+start.within & (age.fhaart>age.cd4.llimit+start.within | is.na(age.fhaart)) & deferred1==1,1,0)) surv.df$cens.start.early<-with(surv.df, ifelse(haart==1 & (age.fhaart=age.cd4.llimit & age=age.cd4.llimit & agew.trunc,w.trunc,d.analysis$weight) with(d.analysis, table(deferred, event)) with(d.analysis, mean(event[deferred==0])) with(d.analysis, mean(event[deferred==1])) chisq.test(with(d.analysis, table(deferred, event))) if (!COVARIATES & !SITE.YEAR){ ### Simple model mod<-glm(event~rcs(month,4) + deferred , weights=weight.trunc, data=d.analysis, family=binomial) } if (COVARIATES & !SITE.YEAR){ ### Model with covariates mod<-glm(event~rcs(month,4) + deferred + ref.cd4 + age.less.ulimit + male, weights=weight.trunc, family="binomial", data=d.analysis) } if (!COVARIATES & SITE.YEAR){ ### Model stratifying by site and year mod<-glm(event~rcs(month,4) + site + year.less.ulimit + deferred, weights=weight.trunc, family="binomial",data=d.analysis) } if (COVARIATES & SITE.YEAR){ ### Model with covariates and stratifying by site and year mod<-glm(event~rcs(month,4) + site + year.less.ulimit + deferred + ref.cd4 + age.less.ulimit + male, #### Model may be too rich weights=weight.trunc, family="binomial", data=d.analysis) } #mod<-glm(event~rcs(month,4)+site+year.less.ulimit + deferred + ref.cd4 + age.less.ulimit + male, # weights=weight.trunc, family="binomial", data=d.analysis) library(sandwich) se<-sqrt(sandwich(mod)["deferred","deferred"]) hr<-exp(summary(mod)$coeff["deferred","Estimate"]) lower<-exp(mod$coeff[names(mod$coeff)=="deferred"]-1.96*se) upper<-exp(mod$coeff[names(mod$coeff)=="deferred"]+1.96*se) ests<-c(hr,lower,upper) save(ests,u.limit,l.limit,LOST,COVARIATES,SITE.YEAR,EXCLUDE.IDU,DEATH.ONLY,CENSOR.TIME,EXCLUDE.DEATHS,start.within,file=filename) } length(unique(surv.df$pid)) length(unique(d.analysis$pid)) length(unique(surv.df$pid[surv.df$cens.no.start==1])) length(unique(surv.df$pid[surv.df$cens.start.early==1])) length(d.analysis$pid[d.analysis$deferred==0])/12 length(d.analysis$pid[d.analysis$deferred==1])/12 sum(d.analysis$event[d.analysis$deferred==0]) sum(d.analysis$event[d.analysis$deferred==1]) with(surv.df,summary(age.cd4.llimit[!duplicated(pid)]-age.ref.cd4[!duplicated(pid)])) with(d.analysis,sum(event[is.na(age.cd4.llimit)&deferred1==1])) ######################## Here I do the analyses #setwd("~/Dropbox/cnics/analyses-2013-Nov") load("cnics-hernan-stuff5c.Rda") d1<-d2[!duplicated(d2$pid),] mean(d1$male) mean(d1$black) summary(d1$idu) summary(d1$first.age) summary(d1$first.year) summary(d1$last.age-d1$first.age) mean(!is.na(d1$age.fhaart)) summary(d1$age.fhaart-d1$first.age)*12 #sum(!is.na(d2$age.oi)) with(d2,length(unique(pid[!is.na(age.oi)]))) sum(!is.na(d1$age_at_death)) #### Basic Analyses hernan.analysis(d2, u.limit=201, l.limit=101, filename="basic-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, filename="basic-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, filename="basic-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, filename="basic-500vs400.Rda") #### Basic Analyses using death as outcome hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, filename="death-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, filename="death-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, filename="death-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, filename="death-500vs400.Rda") #### Administratively Censoring after 6 years hernan.analysis(d2, u.limit=201, l.limit=101, CENSOR.TIME=TRUE, filename="censor6-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, CENSOR.TIME=TRUE, filename="censor6-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, CENSOR.TIME=TRUE, filename="censor6-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, CENSOR.TIME=TRUE, filename="censor6-500vs400.Rda") #### Excluding Deaths in first 2 weeks hernan.analysis(d2, u.limit=201, l.limit=101, EXCLUDE.DEATHS=TRUE, filename="exdeath-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, EXCLUDE.DEATHS=TRUE, filename="exdeath-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, EXCLUDE.DEATHS=TRUE, filename="exdeath-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, EXCLUDE.DEATHS=TRUE, filename="exdeath-500vs400.Rda") #### Excluding IDUs hernan.analysis(d2, u.limit=201, l.limit=101, EXCLUDE.IDU=TRUE, filename="exIDU-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, EXCLUDE.IDU=TRUE, filename="exIDU-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, EXCLUDE.IDU=TRUE, filename="exIDU-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, EXCLUDE.IDU=TRUE, filename="exIDU-500vs400.Rda") #### Window 3 months - RE-RUN hernan.analysis(d2, u.limit=201, l.limit=101, start.within=90/365.25, filename="win3-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, start.within=90/365.25, filename="win3-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, start.within=90/365.25, filename="win3-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, start.within=90/365.25, filename="win3-500vs400.Rda") #### Including Covariates hernan.analysis(d2, u.limit=201, l.limit=101, COVARIATES=TRUE, filename="cov-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, COVARIATES=TRUE, filename="cov-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, COVARIATES=TRUE, filename="cov-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, COVARIATES=TRUE, filename="cov-500vs400.Rda") #### LTFU in weights hernan.analysis(d2, u.limit=201, l.limit=101, LOST=TRUE, filename="lost-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, LOST=TRUE, filename="lost-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, LOST=TRUE, filename="lost-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, LOST=TRUE, filename="lost-500vs400.Rda") #### Stratifying by site and year - RE-RUN hernan.analysis(d2, u.limit=201, l.limit=101, SITE.YEAR=TRUE, filename="strat-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, SITE.YEAR=TRUE, filename="strat-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, SITE.YEAR=TRUE, filename="strat-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, SITE.YEAR=TRUE, filename="strat-500vs400.Rda") ######################################################################### ############# Repeating analyses with death as the outcome #### Censor for >1 year visit gap hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, GAP1YEAR.LTFU=TRUE, filename="gap1year-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, GAP1YEAR.LTFU=TRUE, filename="gap1year-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, GAP1YEAR.LTFU=TRUE, filename="gap1year-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, GAP1YEAR.LTFU=TRUE, filename="gap1year-500vs400-d.Rda") #### Administratively Censoring after 6 years hernan.analysis(d2, u.limit=201, l.limit=101, CENSOR.TIME=TRUE, DEATH.ONLY=TRUE, filename="censor6-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, CENSOR.TIME=TRUE, DEATH.ONLY=TRUE, filename="censor6-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, CENSOR.TIME=TRUE, DEATH.ONLY=TRUE, filename="censor6-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, CENSOR.TIME=TRUE, DEATH.ONLY=TRUE, filename="censor6-500vs400-d.Rda") #### Excluding Deaths in first 2 weeks hernan.analysis(d2, u.limit=201, l.limit=101, EXCLUDE.DEATHS=TRUE, DEATH.ONLY=TRUE, filename="exdeath-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, EXCLUDE.DEATHS=TRUE, DEATH.ONLY=TRUE, filename="exdeath-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, EXCLUDE.DEATHS=TRUE, DEATH.ONLY=TRUE, filename="exdeath-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, EXCLUDE.DEATHS=TRUE, DEATH.ONLY=TRUE, filename="exdeath-500vs400-d.Rda") #### Excluding IDUs hernan.analysis(d2, u.limit=201, l.limit=101, EXCLUDE.IDU=TRUE, DEATH.ONLY=TRUE, filename="exIDU-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, EXCLUDE.IDU=TRUE, DEATH.ONLY=TRUE, filename="exIDU-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, EXCLUDE.IDU=TRUE, DEATH.ONLY=TRUE, filename="exIDU-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, EXCLUDE.IDU=TRUE, DEATH.ONLY=TRUE, filename="exIDU-500vs400-d.Rda") #### Window 3 months hernan.analysis(d2, u.limit=201, l.limit=101, start.within=90/365.25, DEATH.ONLY=TRUE, filename="win3-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, start.within=90/365.25, DEATH.ONLY=TRUE, filename="win3-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, start.within=90/365.25, DEATH.ONLY=TRUE, filename="win3-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, start.within=90/365.25, DEATH.ONLY=TRUE, filename="win3-500vs400-d.Rda") #### Including Covariates hernan.analysis(d2, u.limit=201, l.limit=101, COVARIATES=TRUE, DEATH.ONLY=TRUE, filename="cov-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, COVARIATES=TRUE, DEATH.ONLY=TRUE, filename="cov-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, COVARIATES=TRUE, DEATH.ONLY=TRUE, filename="cov-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, COVARIATES=TRUE, DEATH.ONLY=TRUE, filename="cov-500vs400-d.Rda") #### LTFU in weights hernan.analysis(d2, u.limit=201, l.limit=101, LOST=TRUE, DEATH.ONLY=TRUE, filename="lost-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, LOST=TRUE, DEATH.ONLY=TRUE, filename="lost-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, LOST=TRUE, DEATH.ONLY=TRUE, filename="lost-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, LOST=TRUE, DEATH.ONLY=TRUE, filename="lost-500vs400-d.Rda") #### Stratifying by site and year hernan.analysis(d2, u.limit=201, l.limit=101, SITE.YEAR=TRUE, DEATH.ONLY=TRUE, filename="strat-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, SITE.YEAR=TRUE, DEATH.ONLY=TRUE, filename="strat-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, SITE.YEAR=TRUE, DEATH.ONLY=TRUE, filename="strat-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, SITE.YEAR=TRUE, DEATH.ONLY=TRUE, filename="strat-500vs400-d.Rda") #### Adding 2 months follow-up to those who didn't die hernan.analysis(d2, u.limit=201, l.limit=101, ADD2=TRUE, filename="add2-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, ADD2=TRUE, filename="add2-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, ADD2=TRUE, filename="add2-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, ADD2=TRUE, filename="add2-500vs400.Rda") #### Adding 2 months follow-up to those who didn't die hernan.analysis(d2, u.limit=201, l.limit=101, ADD2=TRUE, DEATH.ONLY=TRUE, filename="add2-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, ADD2=TRUE, DEATH.ONLY=TRUE, filename="add2-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, ADD2=TRUE, DEATH.ONLY=TRUE, filename="add2-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, ADD2=TRUE, DEATH.ONLY=TRUE, filename="add2-500vs400-d.Rda") ####### "Original" NA-ACCORD analyses hernan.analysis(d2,u.limit=501,l.limit=351,COVARIATES=TRUE,LOST=TRUE,SITE.YEAR=TRUE,DEATH.ONLY=TRUE,filename="orig-500vs350-d.Rda") hernan.analysis(d2,u.limit=801,l.limit=501,COVARIATES=TRUE,LOST=TRUE,SITE.YEAR=TRUE,DEATH.ONLY=TRUE,filename="orig-800vs500-d.Rda") ####### "Original" NA-ACCORD analyses, but also including ADE hernan.analysis(d2,u.limit=501,l.limit=351,COVARIATES=TRUE,LOST=TRUE,SITE.YEAR=TRUE,filename="orig-500vs350.Rda") hernan.analysis(d2,u.limit=801,l.limit=501,COVARIATES=TRUE,LOST=TRUE,SITE.YEAR=TRUE,filename="orig-800vs500.Rda") ######## Now playing around with truncation levels for basic death analysis #### No truncation hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, TRUNC.LEVEL=1, filename="trunc-1-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, TRUNC.LEVEL=1, filename="trunc-1-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, TRUNC.LEVEL=1, filename="trunc-1-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, TRUNC.LEVEL=1, filename="trunc-1-500vs400-d.Rda") #### Truncating past 0.995 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-500vs400-d.Rda") #### Truncating past 0.975 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-500vs400-d.Rda") #### Truncating past 0.95 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-500vs400-d.Rda") ######## Now playing around with truncation levels for basic death analysis using covariates #### No truncation hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-500vs400-d.Rda") #### Truncating past 0.995 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-500vs400-d.Rda") #### Truncating past 0.975 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-500vs400-d.Rda") #### Truncating past 0.95 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-500vs400-d.Rda") ######## Now playing around with truncation levels for basic death analysis using covariates and lost #### No truncation hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-500vs400-d.Rda") #### Truncating past 0.995 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-500vs400-d.Rda") #### Truncating past 0.99 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-500vs400-d.Rda") #### Truncating past 0.975 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-500vs400-d.Rda") #### Truncating past 0.95 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-200vs100-d.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-300vs200-d.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-400vs300-d.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=TRUE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-500vs400-d.Rda") ######## Now playing around with truncation levels for basic AIDS/death analysis using covariates #### No truncation hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=1, filename="trunc-1-500vs400.Rda") #### Truncating past 0.995 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.995, filename="trunc-995-500vs400.Rda") #### Truncating past 0.975 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.975, filename="trunc-975-500vs400.Rda") #### Truncating past 0.95 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, TRUNC.LEVEL=.95, filename="trunc-95-500vs400.Rda") ######## Now playing around with truncation levels for basic AIDS/death analysis using covariates and lost #### No truncation hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=1, filename="cov-lost-trunc-1-500vs400.Rda") #### Truncating past 0.995 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.995, filename="cov-lost-trunc-995-500vs400.Rda") #### Truncating past 0.99 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.99, filename="cov-lost-trunc-99-500vs400.Rda") #### Truncating past 0.975 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.975, filename="cov-lost-trunc-975-500vs400.Rda") #### Truncating past 0.95 quantile hernan.analysis(d2, u.limit=201, l.limit=101, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-200vs100.Rda") hernan.analysis(d2, u.limit=301, l.limit=201, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-300vs200.Rda") hernan.analysis(d2, u.limit=401, l.limit=301, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-400vs300.Rda") hernan.analysis(d2, u.limit=501, l.limit=401, DEATH.ONLY=FALSE, COVARIATES=TRUE, LOST=TRUE, TRUNC.LEVEL=.95, filename="cov-lost-trunc-95-500vs400.Rda") ######################## Looking at the results setwd("/Volumes/encrypted/cnics/analyses-2013-Nov") mi.answers<-function(filename){ load(paste(filename,".Rda",sep="")) est<-ests hr<-est[1] lower<-est[2] upper<-est[3] type<-strsplit(filename,"-")[[1]][1] type.long<-ifelse(type=="basic","Basic, ADE/death", ifelse(type=="censor6","Censor after 6 years", ifelse(type=="exdeath","Exclude deaths w/in 2 wks of ART", ifelse(type=="exIDU","Exclude IDU", ifelse(type=="cov","Adjust for sex, base. age and CD4", ifelse(type=="lost","Account for LTFU with weights", ifelse(type=="strat","Stratify by site and year", ifelse(type=="death","Basic, death", ifelse(type=="add2","Add 2 mos. to those w/out death", ifelse(type=="gap1year","Censor if visit gap > 1 year", ifelse(type=="win3","3 mo. CD4 grace period", type))))))))))) stuff<-c(hr,lower,upper,u.limit,l.limit,type,type.long) return(stuff) } ##### Too few events at 800 vs 500, so doesn't work very well. ans.o<-rbind(mi.answers("orig-500vs350-d"),mi.answers("orig-800vs500-d")) ans.o1<-rbind(mi.answers("orig-500vs350"),mi.answers("orig-800vs500")) save(ans.o,ans.o1,file="answers-orig-hernan.Rda") ans<-rbind(mi.answers("basic-200vs100"), mi.answers("basic-300vs200"), mi.answers("basic-400vs300"), mi.answers("basic-500vs400"), mi.answers("exdeath-200vs100"), mi.answers("exdeath-300vs200"), mi.answers("exdeath-400vs300"), mi.answers("exdeath-500vs400"), mi.answers("add2-200vs100"), mi.answers("add2-300vs200"), mi.answers("add2-400vs300"), mi.answers("add2-500vs400"), mi.answers("censor6-200vs100"), mi.answers("censor6-300vs200"), mi.answers("censor6-400vs300"), mi.answers("censor6-500vs400"), mi.answers("exIDU-200vs100"), mi.answers("exIDU-300vs200"), mi.answers("exIDU-400vs300"), mi.answers("exIDU-500vs400"), mi.answers("cov-200vs100"), mi.answers("cov-300vs200"), mi.answers("cov-400vs300"), mi.answers("cov-500vs400"), mi.answers("lost-200vs100"), mi.answers("lost-300vs200"), mi.answers("lost-400vs300"), mi.answers("lost-500vs400"), mi.answers("strat-200vs100"), mi.answers("strat-300vs200"), mi.answers("strat-400vs300"), mi.answers("strat-500vs400")) colnames(ans)=c("hr","lower","upper","u.limit","l.limit","type","type.long") ansd<-rbind(mi.answers("death-200vs100"), mi.answers("death-300vs200"), mi.answers("death-400vs300"), mi.answers("death-500vs400"), mi.answers("exdeath-200vs100-d"), mi.answers("exdeath-300vs200-d"), mi.answers("exdeath-400vs300-d"), mi.answers("exdeath-500vs400-d"), mi.answers("add2-200vs100-d"), mi.answers("add2-300vs200-d"), mi.answers("add2-400vs300-d"), mi.answers("add2-500vs400-d"), mi.answers("censor6-200vs100-d"), mi.answers("censor6-300vs200-d"), mi.answers("censor6-400vs300-d"), mi.answers("censor6-500vs400-d"), mi.answers("exIDU-200vs100-d"), mi.answers("exIDU-300vs200-d"), mi.answers("exIDU-400vs300-d"), mi.answers("exIDU-500vs400-d"), mi.answers("win3-200vs100-d"), mi.answers("win3-300vs200-d"), mi.answers("win3-400vs300-d"), mi.answers("win3-500vs400-d"), mi.answers("gap1year-200vs100-d"), mi.answers("gap1year-300vs200-d"), mi.answers("gap1year-400vs300-d"), mi.answers("gap1year-500vs400-d"), mi.answers("cov-200vs100-d"), mi.answers("cov-300vs200-d"), mi.answers("cov-400vs300-d"), mi.answers("cov-500vs400-d"), mi.answers("lost-200vs100-d"), mi.answers("lost-300vs200-d"), mi.answers("lost-400vs300-d"), mi.answers("lost-500vs400-d"), mi.answers("strat-200vs100-d"), mi.answers("strat-300vs200-d"), mi.answers("strat-400vs300-d"), mi.answers("strat-500vs400-d")) colnames(ansd)=c("hr","lower","upper","u.limit","l.limit","type","type.long") pdf("hernan-results.pdf",width=8,height=4) n.cat<-8 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<-log(0.4) par(mfrow=c(1,4), mar=c(5,1,3,1)) plot(c(-1.1,1.1),c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("ADE/death, 100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-301 plot(c(-1.1,1.1),c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("ADE/death, 200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-401 plot(c(-1.1,1.1),c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("ADE/death, 300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-501 plot(c(-1.1,1.1),c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("ADE/death, 400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } dev.off() pdf("hernan-results1.pdf",width=9,height=4) n.cat<-8 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- -1 xlimits=c(-0.5,1.1) par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],ans[,"type.long"][ans[,"u.limit"]==ulim][i],pos=4,cex=.9) } plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans[,"hr"][ans[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans[,"lower"][ans[,"u.limit"]==ulim][i],ans[,"upper"][ans[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ans[,"type"][ans[,"u.limit"]==ulim][i]) } dev.off() pdf("hernan-results1-d.pdf",width=9,height=4) n.cat<-10 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- -1 xlimits=c(-1.1,1.1) par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],ansd[,"type.long"][ansd[,"u.limit"]==ulim][i],pos=4,cex=.9) } plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ansd[,"hr"][ansd[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ansd[,"lower"][ansd[,"u.limit"]==ulim][i],ansd[,"upper"][ansd[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ansd[,"type"][ansd[,"u.limit"]==ulim][i]) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ansd[,"hr"][ansd[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ansd[,"lower"][ansd[,"u.limit"]==ulim][i],ansd[,"upper"][ansd[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ansd[,"type"][ansd[,"u.limit"]==ulim][i]) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ansd[,"hr"][ansd[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ansd[,"lower"][ansd[,"u.limit"]==ulim][i],ansd[,"upper"][ansd[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ansd[,"type"][ansd[,"u.limit"]==ulim][i]) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ansd[,"hr"][ansd[,"u.limit"]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ansd[,"lower"][ansd[,"u.limit"]==ulim][i],ansd[,"upper"][ansd[,"u.limit"]==ulim][i]))),rep(y.pos[i],2)) #text(text.pos,y.pos[i],ansd[,"type"][ansd[,"u.limit"]==ulim][i]) } dev.off() trunc.answers<-function(filename){ load(paste(filename,".Rda",sep="")) est<-ests hr<-est[1] lower<-est[2] upper<-est[3] type<-strsplit(filename,"-")[[1]][2] stuff<-c(hr,lower,upper,u.limit,l.limit,type) return(stuff) } ans.trunc1<-rbind(trunc.answers("trunc-95-200vs100-d"), trunc.answers("trunc-95-300vs200-d"), trunc.answers("trunc-95-400vs300-d"), trunc.answers("trunc-95-500vs400-d"), trunc.answers("trunc-975-200vs100-d"), trunc.answers("trunc-975-300vs200-d"), trunc.answers("trunc-975-400vs300-d"), trunc.answers("trunc-975-500vs400-d"), trunc.answers("trunc-995-200vs100-d"), trunc.answers("trunc-995-300vs200-d"), trunc.answers("trunc-995-400vs300-d"), trunc.answers("trunc-995-500vs400-d"), trunc.answers("trunc-1-200vs100-d"), trunc.answers("trunc-1-300vs200-d"), trunc.answers("trunc-1-400vs300-d"), trunc.answers("trunc-1-500vs400-d"))[,1:6] stuff<-rbind(mi.answers("cov-200vs100-d"), mi.answers("cov-300vs200-d"), mi.answers("cov-400vs300-d"), mi.answers("cov-500vs400-d"))[,1:5] stuff<-cbind(stuff,c(99,99,99,99)) ans.trunc<-rbind(ans.trunc1[1:8,],stuff,ans.trunc1[9:16,]) pdf("hernan-truncation-cov-d.pdf",width=9,height=4) n.cat<-5 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- 0 xlimits=c(-0.5,1.1) text.trunc<-c(95,97.5,99,99.5,"None") par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],text.trunc[i],pos=4,cex=1.2) } text(0.2,.95,"Truncation Percentile",cex=1.2) ulim<-201 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } dev.off() ####### Repeating but for ADE/death ans.trunc1<-rbind(trunc.answers("trunc-95-200vs100"), trunc.answers("trunc-95-300vs200"), trunc.answers("trunc-95-400vs300"), trunc.answers("trunc-95-500vs400"), trunc.answers("trunc-975-200vs100"), trunc.answers("trunc-975-300vs200"), trunc.answers("trunc-975-400vs300"), trunc.answers("trunc-975-500vs400"), trunc.answers("trunc-995-200vs100"), trunc.answers("trunc-995-300vs200"), trunc.answers("trunc-995-400vs300"), trunc.answers("trunc-995-500vs400"), trunc.answers("trunc-1-200vs100"), trunc.answers("trunc-1-300vs200"), trunc.answers("trunc-1-400vs300"), trunc.answers("trunc-1-500vs400"))[,1:6] stuff<-rbind(mi.answers("cov-200vs100"), mi.answers("cov-300vs200"), mi.answers("cov-400vs300"), mi.answers("cov-500vs400"))[,1:5] stuff<-cbind(stuff,c(99,99,99,99)) ans.trunc<-rbind(ans.trunc1[1:8,],stuff,ans.trunc1[9:16,]) pdf("hernan-truncation-cov.pdf",width=9,height=4) n.cat<-5 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- 0 xlimits=c(-0.5,1.1) text.trunc<-c(95,97.5,99,99.5,"None") par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],text.trunc[i],pos=4,cex=1.2) } text(0.2,.95,"Truncation Percentile",cex=1.2) ulim<-201 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } dev.off() ####### Repeating but for ADE/death ans.trunc<-rbind(trunc.answers("cov-lost-trunc-95-200vs100"), trunc.answers("cov-lost-trunc-95-300vs200"), trunc.answers("cov-lost-trunc-95-400vs300"), trunc.answers("cov-lost-trunc-95-500vs400"), trunc.answers("cov-lost-trunc-975-200vs100"), trunc.answers("cov-lost-trunc-975-300vs200"), trunc.answers("cov-lost-trunc-975-400vs300"), trunc.answers("cov-lost-trunc-975-500vs400"), trunc.answers("cov-lost-trunc-99-200vs100"), trunc.answers("cov-lost-trunc-99-300vs200"), trunc.answers("cov-lost-trunc-99-400vs300"), trunc.answers("cov-lost-trunc-99-500vs400"), trunc.answers("cov-lost-trunc-995-200vs100"), trunc.answers("cov-lost-trunc-995-300vs200"), trunc.answers("cov-lost-trunc-995-400vs300"), trunc.answers("cov-lost-trunc-995-500vs400"), trunc.answers("cov-lost-trunc-1-200vs100"), trunc.answers("cov-lost-trunc-1-300vs200"), trunc.answers("cov-lost-trunc-1-400vs300"), trunc.answers("cov-lost-trunc-1-500vs400"))[,1:6] pdf("hernan-truncation-cov-lost.pdf",width=9,height=4) n.cat<-5 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- 0 xlimits=c(-0.5,1.1) text.trunc<-c(95,97.5,99,99.5,"None") par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],text.trunc[i],pos=4,cex=1.2) } text(0.2,.95,"Truncation Percentile",cex=1.2) ulim<-201 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } dev.off() ####### Repeating but for death ans.trunc<-rbind(trunc.answers("cov-lost-trunc-95-200vs100-d"), trunc.answers("cov-lost-trunc-95-300vs200-d"), trunc.answers("cov-lost-trunc-95-400vs300-d"), trunc.answers("cov-lost-trunc-95-500vs400-d"), trunc.answers("cov-lost-trunc-975-200vs100-d"), trunc.answers("cov-lost-trunc-975-300vs200-d"), trunc.answers("cov-lost-trunc-975-400vs300-d"), trunc.answers("cov-lost-trunc-975-500vs400-d"), trunc.answers("cov-lost-trunc-99-200vs100-d"), trunc.answers("cov-lost-trunc-99-300vs200-d"), trunc.answers("cov-lost-trunc-99-400vs300-d"), trunc.answers("cov-lost-trunc-99-500vs400-d"), trunc.answers("cov-lost-trunc-995-200vs100-d"), trunc.answers("cov-lost-trunc-995-300vs200-d"), trunc.answers("cov-lost-trunc-995-400vs300-d"), trunc.answers("cov-lost-trunc-995-500vs400-d"), trunc.answers("cov-lost-trunc-1-200vs100-d"), trunc.answers("cov-lost-trunc-1-300vs200-d"), trunc.answers("cov-lost-trunc-1-400vs300-d"), trunc.answers("cov-lost-trunc-1-500vs400-d"))[,1:6] pdf("hernan-truncation-cov-lost-d.pdf",width=9,height=4) n.cat<-5 y.pos<-1-c(1:n.cat)/(n.cat+1) ulim<-201 text.pos<- 0 xlimits=c(-1.1,1.1) text.trunc<-c(95,97.5,99,99.5,"None") par(mfrow=c(1,5), mar=c(5,.5,3,.2)) plot(c(-1,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) for (i in 1:n.cat){ text(text.pos,y.pos[i],text.trunc[i],pos=4,cex=1.2) } text(0.2,.95,"Truncation Percentile",cex=1.2) ulim<-201 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("100 vs. 200 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-301 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("200 vs. 300 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-401 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("300 vs. 400 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } ulim<-501 plot(xlimits,c(0,1),type="n",xlab="Hazard Ratio",ylab="",axes=FALSE) title("400 vs. 500 (ref)") axis(1,at=c(log(0.33),log(0.5),0,log(2),log(3)),c(0.3,0.5,1,2,3)) abline(v=0,col=gray(.8)) points(log(as.numeric(ans.trunc[,1][ans.trunc[,4]==ulim])),y.pos,pch=3) for (i in 1:n.cat){ lines(log(as.numeric(c(ans.trunc[,2][ans.trunc[,4]==ulim][i],ans.trunc[,3][ans.trunc[,4]==ulim][i]))),rep(y.pos[i],2)) } dev.off()