rm(list=ls()) library(foreign) cl<-read.dta("when_to_start_03dec2009.dta") nade<-read.csv("w2s_non_ades_dx_18nov2009b.csv") tx<-read.csv("when_to_start_art_03dec2009.csv") oi<-read.csv("when_to_start_year1_oi_18nov2009.csv") prior.art<-read.csv("when_to_start_prior_art_20nov09.csv") prior.oi<-read.csv("when_to_start_prior_oi_19apr2010.csv") ## Treatment data select<-tx$FIRST_HAART==1&tx$REGIMEN_NBR==1 id<-tx$CFAR_PID[select] reg<-tx$REGIMEN[select] age.haart<-tx$AGE_AT_RX_START[select] trt<-data.frame(id,reg,age.haart) ## Baseline demographics and death/ltfu select<-!duplicated(cl$cfar_pid) id<-cl$cfar_pid[select] sex<-cl$sex[select] male<-sex-1 black<-cl$black[select] white<-cl$white[select] deceased<-cl$deceased[select] age.death<-cl$age_at_death[select] idu<-cl$ivdu[select] route<-cl$probable_route_of_infection[select] ltfu<-cl$lost_to_fu[select] age.first.visit<-cl$age_at_first_visit[select] age.last<-year.haart<-age.haart1<-NULL for (i in 1:length(id)){ age.last[i]<-max(cl$age[cl$cfar_pid==id[i]]) year.haart[i]<-NA if (sum(cl$cfar_pid==id[i]&cl$indicator_first_haart_0_1==1)>0) { year.haart[i]<-cl$year_of_first_haart[cl$cfar_pid==id[i]&cl$indicator_first_haart_0_1==1] age.haart1[i]<-cl$age_at_first_haart[cl$cfar_pid==id[i]&cl$indicator_first_haart_0_1==1] } } demo<-data.frame(id,male,black,white,age.haart1,age.first.visit,idu,route,year.haart,ltfu,age.last,deceased,age.death) ## Lab values choose<-!is.na(cl$age_at_lab) cd4<-cl$cd4_count[choose] pcd4<-cl$cd4_percent[choose] age.lab<-cl$age_at_lab[choose] vl<-cl$vl[choose] id<-cl$cfar_pid[choose] cd4.lab<-ifelse(!is.na(cd4),1,0) vl.lab<-ifelse(!is.na(vl),1,0) labs<-data.frame(id,age.lab,cd4.lab,vl.lab,cd4,pcd4,vl) d1<-merge(demo,labs,by.x="id",by.y="id",all=TRUE) d2<-merge(d1,trt,by.x="id",by.y="id",all=TRUE) diff<-d2$age.haart-d2$age.lab k<-.5 ### 3 months = 1/4 year; 6 months= 1/2 year. stuff<-cd4.base<-pcd4.base<-vl.base<-NULL unique.id<-unique(d2$id) for (i in 1:length(unique.id)) { stuff[i]<-cd4.base[i]<-vl.base[i]<-pcd4.base<-NA if (sum(d2$id==unique.id[i]&diff>0&diff0) { stuff[i]<-min(diff[d2$id==unique.id[i]&diff>0&diff0&diff0) { stuff[i]<-min(diff[d2$id==unique.id[i]&diff>0&diff0,] ###### OIs oi.id<-oi$CFAR_PID unique.oi.id<-unique(oi.id) stuff<-NULL for (i in 1:length(unique.oi.id)) { stuff[i]<-NA if (sum(oi$INDICATOR_OI_0_1[oi.id==unique.oi.id[i]])>0){ stuff[i]<-min(oi$AGE_AT_OI[oi.id==unique.oi.id[i]]) } } ade<-ifelse(!is.na(stuff),1,0) age.ade<-stuff oi.d<-data.frame(id=unique.oi.id,ade,age.ade) ####### Merging d.ade1<-merge(d5,oi.d,by="id",all.x=TRUE) exclude<-ifelse(is.na(d.ade1$age.ade),0,ifelse(d.ade1$age.ade<=d.ade1$age.haart,1,0)) d.ade<-d.ade1[exclude==0,] d1<-d.ade d2<-merge(d1,prior.oi,by.x="id",by.y="CFAR_PID",all.x=TRUE) d3<-merge(d2,prior.art,by.x="id",by.y="CFAR_PID",all.x=TRUE) exclude1<-ifelse(d3$HAS_OI_BEFORE_FIRST_VISIT==1|d3$HAS_ART_BEFORE_FIRST_VISIT==1,1,0) d<-d3[exclude1==0,] id<-d$id male<-d$male black<-d$black white<-d$white idu<-d$idu route<-d$route year.haart<-d$year.haart ltfu<-d$ltfu age.last<-d$age.last deceased<-d$deceased age.death<-d$age.death cd4<-d$cd4 pcd4<-d$pcd4 vl<-d$vl reg<-d$reg age.haart<-d$age.haart time<-d$time ade<-d$ade age.ade<-d$age.ade logvl1<-log10(vl) logvl<-ifelse(is.na(logvl1)&!is.na(cd4),mean(logvl1,na.rm=TRUE),logvl1) age.first.visit<-d$age.first.visit time.in.care<-round(age.haart-age.first.visit,2) cd4.401t500<-ifelse(cd4>400&cd4<=500,1,0) cd4.301t400<-ifelse(cd4>300&cd4<=400,1,0) cd4.201t300<-ifelse(cd4>200&cd4<=300,1,0) cd4.101t200<-ifelse(cd4>100&cd4<=200,1,0) ############################################################################# ###################################################### ###### ADE/Mortality analysis library(survival) event<-ifelse(ade==1,1,deceased) time.to.ade<-age.ade-age.haart time.ade<-NULL for (i in 1:length(id)) { time.ade[i]<-min(time.to.ade[i],time[i],na.rm=TRUE) } S.ade<-Surv(time.ade,event) ######################## Now adding lead times ########## Simulating generalized gamma data qggamma <- function(beta, sigma, lambda, p) { s<-lambda^(-2) r<-s*exp((-beta*lambda)/sigma) b<-lambda/sigma quan<-(lambda>0)*((qgamma(p, shape=s, rate=r))^(1/b)) quan<- quan+(lambda<0)*((qgamma(1-p,shape=s, rate=r))^(1/b)) quan } ##### Imputing for lead times for AIDS/death par.ests<- matrix(c(-0.03244,1.4193,0.7969,0.4088,1.0706,1.2097,0.01610,450, 0.1132,1.4260,0.8284,0.5776,0.9573,1.4753,0.02273,425, 0.1211,1.4686,0.6908,0.5509,1.0282,1.3532,0.02691,400, 0.1938,1.4293,0.7916,0.1801,1.1994,0.8842,0.02830,375, 0.2147,1.3643,0.8099,0.2086,2.3632,0.8688,0.03197,350, 0.2830,1.2848,0.9275,0.002276,1.1950,0.8346,0.03249,325, 0.2590,1.2710,0.8762,0.07611,1.3065,0.6431,0.03828,300, 0.3509,1.2137,1.0532,0.3388,1.4096,0.7650,0.05289,275, 0.3211,1.2222,0.9638,0.5177,1.2067,1.1048,0.06128,250, 0.2951,1.2273,0.9177,0.4338,1.1252,1.2341,0.06974,225, 0.3231,1.1839,0.9771,0.3187,1.2826,1.0666,0.08119,200, 0.2906,1.1298,0.9100,0.2165,1.1235,1.1676,0.09472,175, 0.3235,1.0885,0.9236,0.1419,1.0550,1.2139,0.1079,150, 0.2989,1.0653,0.8941,0.1774,1.0114,1.2725,0.1309,125, 0.2417,1.0461,0.8189,0.1073,1.1676,1.0102,0.1759,100),ncol=8,byrow=TRUE) meds<-NULL for (i in 1:15) { meds[i]<-median(qggamma(par.ests[i,1],par.ests[i,2],par.ests[i,3],runif(10000,0,1))) } hr.compute.ade.cov.lt<-function(lower,upper,defer.cd4,n.impute) { select1<-ifelse(lower==1|upper==1,1,0) select<-ifelse(is.na(select1),0,select1) t<-time.ade[select==1] event.orig<-event[select==1] defer.orig<-lower[select==1] row<-which(par.ests[,8]==defer.cd4) mod.wgt<-glm(lower~idu+black+male+age.haart+logvl+sqrt(time.in.care),family="binomial",subset=select==1) pred.logit<-predict(mod.wgt) pred<-exp(pred.logit)/(1+exp(pred.logit)) p1<-pred*lower[select==1] + (1-pred)*(1-lower[select==1]) p<-ifelse(p1quantile(p1,.95),quantile(p1,.95),p1)) w.orig<-((1/p)/sum(1/p))*length(p) ### To get the weights right pi.est<-par.ests[row,7] n<-sum(defer.orig) beta<-se<-beta.w<-se.w<-beta.ps<-se.ps<-n.unseen<-NULL for (i in 1:n.impute) { event<-event.orig defer<-defer.orig w<-w.orig ps<-pred t.f<-ifelse(defer==0,0,qggamma(par.ests[row,1],par.ests[row,2],par.ests[row,3],runif(n,0,1))) t.star<-t+t.f u<-rbinom(1,n,pi.est/(1-pi.est)) if (u>0) { t.g<-qggamma(par.ests[row,4],par.ests[row,5],par.ests[row,6],runif(u,0,1)) t.star<-c(t.star,t.g) defer<-c(defer,rep(1,u)) event<-c(event,rep(1,u)) samp<-sample(c(1:length(w.orig[defer.orig==1])),u,replace=TRUE) w<-c(w,w.orig[defer.orig==1][samp]) ps<-c(ps,pred[defer.orig==1][samp]) } S<-Surv(t.star,event) mod<-coxph(S~defer) beta[i]<-summary(mod)$coeff[1] se[i]<-summary(mod)$coeff[3] n.unseen[i]<-u mod.w<-coxph(S~defer,weight=w) beta.w[i]<-summary(mod.w)$coeff[1] se.w[i]<-summary(mod.w)$coeff[3] mod.ps<-coxph(S~defer+ps) beta.ps[i]<-summary(mod.ps)$coeff[1] se.ps[i]<-summary(mod.ps)$coeff[1,3] } hr<-exp(mean(beta)) var<-mean(se^2)+var(beta)*(1+1/n.impute) ci<-exp(mean(beta)+1.96*sqrt(var)*c(-1,1)) hr.w<-exp(mean(beta.w)) var.w<-mean(se.w^2)+var(beta.w)*(1+1/n.impute) ci.w<-exp(mean(beta.w)+1.96*sqrt(var.w)*c(-1,1)) hr.ps<-exp(mean(beta.ps)) var.ps<-mean(se.ps^2)+var(beta.ps)*(1+1/n.impute) ci.ps<-exp(mean(beta.ps)+1.96*sqrt(var.ps)*c(-1,1)) u.mean<-mean(n.unseen) N<-length(t) return(c(hr,ci,hr.w,ci.w,hr.ps,ci.ps,u.mean,N,n)) } ### Large numbers of imputations set.seed(812) hr.compute.ade.cov.lt(cd4.301t400,cd4.401t500,400,1000) hr.compute.ade.cov.lt(cd4.201t300,cd4.301t400,300,1000) hr.compute.ade.cov.lt(cd4.101t200,cd4.201t300,200,1000)