### This file creates figures and includes computation used in the paper.  The HR that are hard-coded (unnecessarily!) 
#   at the beginning of the file are taken from computations done in the other analysis file.

### Taken from file iwhod-figure.R

rm(list=ls())



NAME.ade<-c(1.21,0.96,1.53)
NAME.ade.death<-c(1.09, 0.88, 1.34)
NAME.ltfu<-c(0.71,0.56,0.89)
NAME.ade.ntb<-c(1.11,0.86,1.44)
NAME.ade.death.ntb<-c(1.00, 0.80, 1.26)
NAME.death<-c(0.79, 0.52, 1.20)

SSA.ade<-c(1.25, 1.10, 1.42)
SSA.ade.death<-c(1.17, 1.05, 1.31)
SSA.ltfu<-c(1.08, 1.00, 1.18)
SSA.ade.ntb<-c(1.12, 0.97, 1.28)
SSA.ade.death.ntb<-c(1.06, 0.94, 1.20)
SSA.death<-c(0.91, 0.72, 1.15)

ASIA.ade<-c(0.99,0.72, 1.37)
ASIA.ade.death<-c(0.94, 0.70,1.26)
ASIA.ltfu<-c(0.98, 0.75, 1.29)
ASIA.ade.ntb<-c(1.00,0.72, 1.40)
ASIA.ade.death.ntb<-c(0.94, 0.69,1.27)
ASIA.death<-c(0.93, 0.53, 1.61)

LA.ade<-c(1.17, 0.98, 1.39)
LA.ade.death<-c(1.12, 0.96, 1.31)
LA.ltfu<-c(1.27,1.12,1.43)
LA.ade.ntb<-c(1.12, 0.92, 1.35)
LA.ade.death.ntb<-c(1.08, 0.92, 1.27)
LA.death<-c(0.98, 0.73, 1.32)

nn.ade<-c(1.21, 1.09, 1.34)
nn.ade.death<-c(1.13, 1.03, 1.24)
nn.ltfu<-c(1.08, 1.00, 1.16)
nn.ade.ntb<-c(1.11, 0.99, 1.24)
nn.ade.death.ntb<-c(1.05, 0.95, 1.16)
nn.death<-c(0.91, 0.76, 1.09)

nn.tb<-c(1.94, 1.53, 2.46)
NAME.tb<-c(2.09, 1.26, 3.47)
SSA.tb<-c(2.17, 1.65, 2.86)
ASIA.tb<-c(0.97, 0.36, 2.64)
LA.tb<-c(1.64, 1.10, 2.46)



#### After Jonathan's suggestions

#pdf("iwhod-figure.pdf")
#pdf("figure2-jonathan.pdf")
postscript("figure2-paper.eps",width=8,height=8)

k<- -.6
k1<-.8
k2<-.1
a<-.8
lw<-1.5


#par(mar=c(3,.5,1,1))
plot(c(k-1.4,.6+k1),c(17,37),type="n",xlab="",ylab="",bty="n",axes=FALSE)
abline(v=0,col=gray(.8))


text(k-1.5,37,"ADE",pos=4,cex=a)
text(k-1.4,36,"Migrants",pos=4,cex=a)
text(k-1.3,35,"Sub-Saharan Africa",pos=4,cex=a)
text(k-1.3,34,"Latin America",pos=4,cex=a)
text(k-1.3,33,"North Africa/Middle East",pos=4,cex=a)
text(k-1.3,32,"Asia",pos=4,cex=a)

lines(log(c(nn.ade[2],nn.ade[3])),c(36,36),lwd=lw)
lines(log(c(SSA.ade[2],SSA.ade[3])),c(35,35),lwd=lw)
lines(log(c(LA.ade[2],LA.ade[3])),c(34,34),lwd=lw)
lines(log(c(NAME.ade[2],NAME.ade[3])),c(33,33),lwd=lw)
lines(log(c(ASIA.ade[2],ASIA.ade[3])),c(32,32),lwd=lw)

points(log(nn.ade[1]),36,pch=20,cex=1.5)
points(log(SSA.ade[1]),35,pch=20,cex=1.5)
points(log(LA.ade[1]),34,pch=20,cex=1.5)
points(log(NAME.ade[1]),33,pch=20,cex=1.5)
points(log(ASIA.ade[1]),32,pch=20,cex=1.5)



text(k-1.5,30,"ADE without TB",pos=4,cex=a)
text(k-1.4,29,"Migrants",pos=4,cex=a)
text(k-1.3,28,"Sub-Saharan Africa",pos=4,cex=a)
text(k-1.3,27,"Latin America",pos=4,cex=a)
text(k-1.3,26,"North Africa/Middle East",pos=4,cex=a)
text(k-1.3,25,"Asia",pos=4,cex=a)


lines(log(c(nn.ade.ntb[2],nn.ade.ntb[3])),c(29,29),lwd=lw)
lines(log(c(SSA.ade.ntb[2],SSA.ade.ntb[3])),c(28,28),lwd=lw)
lines(log(c(LA.ade.ntb[2],LA.ade.ntb[3])),c(27,27),lwd=lw)
lines(log(c(NAME.ade.ntb[2],NAME.ade.ntb[3])),c(26,26),lwd=lw)
lines(log(c(ASIA.ade.ntb[2],ASIA.ade.ntb[3])),c(25,25),lwd=lw)

points(log(nn.ade.ntb[1]),29,pch=20,cex=1.5)
points(log(SSA.ade.ntb[1]),28,pch=20,cex=1.5)
points(log(LA.ade.ntb[1]),27,pch=20,cex=1.5)
points(log(NAME.ade.ntb[1]),26,pch=20,cex=1.5)
points(log(ASIA.ade.ntb[1]),25,pch=20,cex=1.5)



text(k-1.5,23,"Death",pos=4,cex=a)
text(k-1.4,22,"Migrants",pos=4,cex=a)
text(k-1.3,21,"Sub-Saharan Africa",pos=4,cex=a)
text(k-1.3,20,"Latin America",pos=4,cex=a)
text(k-1.3,19,"North Africa/Middle East",pos=4,cex=a)
text(k-1.3,18,"Asia",pos=4,cex=a)

lines(log(c(nn.death[2],nn.death[3])),c(22,22),lwd=lw)
lines(log(c(SSA.death[2],SSA.death[3])),c(21,21),lwd=lw)
lines(log(c(LA.death[2],LA.death[3])),c(20,20),lwd=lw)
lines(log(c(NAME.death[2],NAME.death[3])),c(19,19),lwd=lw)
lines(log(c(ASIA.death[2],ASIA.death[3])),c(18,18),lwd=lw)

points(log(nn.death[1]),22,pch=20,cex=1.5)
points(log(SSA.death[1]),21,pch=20,cex=1.5)
points(log(LA.death[1]),20,pch=20,cex=1.5)
points(log(NAME.death[1]),19,pch=20,cex=1.5)
points(log(ASIA.death[1]),18,pch=20,cex=1.5)




text(k1-k2,36,paste(nn.ade[1]," (",paste(nn.ade[2],nn.ade[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,35,paste(SSA.ade[1]," (",paste(SSA.ade[2],SSA.ade[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,34,paste(LA.ade[1]," (",paste(LA.ade[2],LA.ade[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,33,paste(NAME.ade[1]," (",paste(NAME.ade[2],NAME.ade[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,32,paste(ASIA.ade[1]," (",paste(ASIA.ade[2],ASIA.ade[3],sep=", "),")",sep=""),pos=4,cex=a)


text(k1-k2,29,paste(nn.ade.ntb[1]," (",paste(nn.ade.ntb[2],nn.ade.ntb[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,28,paste(SSA.ade.ntb[1]," (",paste(SSA.ade.ntb[2],SSA.ade.ntb[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,27,paste(LA.ade.ntb[1]," (",paste(LA.ade.ntb[2],LA.ade.ntb[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,26,paste(NAME.ade.ntb[1]," (",paste(NAME.ade.ntb[2],NAME.ade.ntb[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,25,paste(round(ASIA.ade.ntb[1],2)," (",paste(round(ASIA.ade.ntb[2],2),round(ASIA.ade.ntb[3],2),sep=", "),")",sep=""),pos=4,cex=a)


text(k1-k2,22,paste(nn.death[1]," (",paste(nn.death[2],nn.death[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,21,paste(SSA.death[1]," (",paste(SSA.death[2],SSA.death[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,20,paste(LA.death[1]," (",paste(LA.death[2],LA.death[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,19,paste(NAME.death[1]," (",paste(NAME.death[2],NAME.death[3],sep=", "),")",sep=""),pos=4,cex=a)
text(k1,18,paste(ASIA.death[1]," (",paste(ASIA.death[2],ASIA.death[3],sep=", "),")",sep=""),pos=4,cex=a)



mtext("Adjusted Hazard Ratio",side=1,at=0,line=3,cex=a)
#mtext("Adjusted Hazard Ratio",side=3,at=0+k1,line=1,cex=a)
#mtext("(without TB)",side=3,at=0+k1,line=0,cex=a)

axis(1,at=log(c(0.5,.67,1,1.5,2)),labels=c(0.5,0.67,1,1.5,2.0),cex.axis=a)


dev.off()



#### Taken from file incidence-ci.R

rm(list=ls())

load("sub_baseline_for_ci.rda")


#####  Comparing those with and without NB status
d2<-sub.baseline[sub.baseline$cohort!="EuroSIDA"&sub.baseline$cohort!="HAVACS"&sub.baseline$cohort!="ICONA"&sub.baseline$cohort!="Royal Free"&
                   sub.baseline$cohort!="UAB"&sub.baseline$cohort!="VACS",]

summary(d2$haart_d)
unique.cohort<-unique(d2$cohort)
earliest.sd<-NULL
for (i in 1:length(unique.cohort)) {
  earliest.sd[i]<-as.character(min(as.Date(d2$haart_d[d2$cohort==unique.cohort[i]])))
}
data.frame(unique.cohort,earliest.sd)

d2$missing.nb<-ifelse(is.na(d2$native.born.factor),1,0)

with(d2, table(missing.nb,gender.factor))
with(d2, table(missing.nb,gender.factor)[,1]/apply(table(missing.nb,gender.factor),1,sum))
with(d2, chisq.test(table(missing.nb,gender.factor)))

with(d2, summary(cd4_t0_corrected[missing.nb==1]))
with(d2, summary(cd4_t0_corrected[missing.nb==0]))
with(d2, wilcox.test(cd4_t0_corrected[missing.nb==0],cd4_t0_corrected[missing.nb==1]))

with(d2, summary(lrna_t0_corrected[missing.nb==1]))
with(d2, summary(lrna_t0_corrected[missing.nb==0]))
with(d2, wilcox.test(lrna_t0_corrected[missing.nb==0],lrna_t0_corrected[missing.nb==1]))

with(d2, summary(age[missing.nb==1]))
with(d2, summary(age[missing.nb==0]))
with(d2, wilcox.test(age[missing.nb==0],age[missing.nb==1]))

with(d2, summary(haartyr[missing.nb==1]))
with(d2, summary(haartyr[missing.nb==0]))
with(d2, wilcox.test(haartyr[missing.nb==0],haartyr[missing.nb==1]))

with(d2, table(missing.nb,stagec.factor))
1-with(d2, table(missing.nb,stagec.factor)[,1]/apply(table(missing.nb,stagec.factor),1,sum))
with(d2, chisq.test(table(missing.nb,stagec.factor)))

with(d2, table(missing.nb,pi_based.factor))
with(d2, table(missing.nb,pi_based.factor)[,1]/apply(table(missing.nb,pi_based.factor),1,sum))
with(d2, chisq.test(table(missing.nb,pi_based.factor)))

with(d2, table(missing.nb,nnrti_based.factor))
with(d2, table(missing.nb,nnrti_based.factor)[,1]/apply(table(missing.nb,nnrti_based.factor),1,sum))
with(d2, chisq.test(table(missing.nb,nnrti_based.factor)))

with(d2, table(missing.nb,risk.factor))
with(d2, table(missing.nb,risk.factor)[,1]/apply(table(missing.nb,risk.factor)[,-4],1,sum))
with(d2, table(missing.nb,risk.factor)[,2]/apply(table(missing.nb,risk.factor)[,-4],1,sum))
with(d2, table(missing.nb,risk.factor)[,3]/apply(table(missing.nb,risk.factor)[,-4],1,sum))
with(d2, table(missing.nb,risk.factor)[,2]/apply(table(missing.nb,risk.factor),1,sum))
with(d2, chisq.test(table(missing.nb,risk.factor)))



d1<-sub.baseline[!is.na(sub.baseline$native.born.factor),]

library(foreign)
d.cod<-read.dta("cause_of_death100811.dta")

d<-merge(d1,d.cod,by=c("patient","cohort"),all.x=TRUE)




fup<-d$incidence_fup


table(d$event1_name)
table(d$event2_name)
table(d$event3_name)
table(d$event4_name)
table(d$event8_name)

time1<-with(d, as.numeric(event_d1)-as.numeric(haart_d))
time2<-with(d, as.numeric(event_d2)-as.numeric(haart_d))
time3<-with(d, as.numeric(event_d3)-as.numeric(haart_d))
time4<-with(d, as.numeric(event_d4)-as.numeric(haart_d))
time5<-with(d, as.numeric(event_d5)-as.numeric(haart_d))
time6<-with(d, as.numeric(event_d6)-as.numeric(haart_d))
time7<-with(d, as.numeric(event_d7)-as.numeric(haart_d))

clin.stop.time<-as.numeric(d$clinfup_d)-as.numeric(d$haart_d)

num.unique<-function(x) {
  length(unique(x[x!=9999]))
}

make.ade.var<-function(x) {
  stuff1<-d$event1==x&!is.na(time1)&time1<=365&time1<=clin.stop.time
  stuff2<-d$event2==x&!is.na(time2)&time2<=365&time2<=clin.stop.time
  stuff3<-d$event3==x&!is.na(time3)&time3<=365&time3<=clin.stop.time
  stuff4<-d$event4==x&!is.na(time4)&time4<=365&time4<=clin.stop.time
  stuff5<-d$event5==x&!is.na(time5)&time5<=365&time5<=clin.stop.time
  stuff6<-d$event6==x&!is.na(time6)&time6<=365&time6<=clin.stop.time
  stuff7<-d$event7==x&!is.na(time7)&time7<=365&time7<=clin.stop.time
  
  times1<-ifelse(stuff1,time1,9999)
  times2<-ifelse(stuff2,time2,9999)
  times3<-ifelse(stuff3,time3,9999)
  times4<-ifelse(stuff4,time4,9999)
  times5<-ifelse(stuff5,time5,9999) 
  times6<-ifelse(stuff6,time6,9999)
  times7<-ifelse(stuff7,time7,9999)
  
  n.events<-apply(cbind(times1,times2,times3,times4,times5,times6,times7),1,num.unique)
  #  n.events<-apply(cbind(stuff1,stuff2,stuff3,stuff4,stuff5,stuff6,stuff7),1,sum)
  t.event<-apply(cbind(times1,times2,times3,times4,times5,times6,times7),1,min)
  
  list(n.events,t.event)
}

a.var<-make.ade.var("BPN")
a.bpn<-a.var[[1]]
t.bpn<-a.var[[2]]
a.var<-make.ade.var("CAN")
a.can<-a.var[[1]]
t.can<-a.var[[2]]
a.var<-make.ade.var("CCM")
a.ccm<-a.var[[1]]
t.ccm<-a.var[[2]]
a.var<-make.ade.var("CMV")
a.cmv<-a.var[[1]]
t.cmv<-a.var[[2]]
a.var<-make.ade.var("COC")
a.coc<-a.var[[1]]
t.coc<-a.var[[2]]
a.var<-make.ade.var("DEM")
a.dem<-a.var[[1]]
t.dem<-a.var[[2]]
a.var<-make.ade.var("HGL")
a.hgl<-a.var[[1]]
t.hgl<-a.var[[2]]
a.var<-make.ade.var("HIS")
a.his<-a.var[[1]]
t.his<-a.var[[2]]
a.var<-make.ade.var("HSD")
a.hsd<-a.var[[1]]
t.hsd<-a.var[[2]]
a.var<-make.ade.var("ICC")
a.icc<-a.var[[1]]
t.icc<-a.var[[2]]
a.var<-make.ade.var("ILE")
a.ile<-a.var[[1]]
t.ile<-a.var[[2]]
a.var<-make.ade.var("ISO")
a.iso<-a.var[[1]]
t.iso<-a.var[[2]]
a.var<-make.ade.var("KSA")
a.ksa<-a.var[[1]]
t.ksa<-a.var[[2]]
a.var<-make.ade.var("LEI")
a.lei<-a.var[[1]]
t.lei<-a.var[[2]]
a.var<-make.ade.var("LOB")
a.lob<-a.var[[1]]
t.lob<-a.var[[2]]
a.var<-make.ade.var("MAC")
a.mac<-a.var[[1]]
t.mac<-a.var[[2]]
a.var<-make.ade.var("MYC")
a.myc<-a.var[[1]]
t.myc<-a.var[[2]]
a.var<-make.ade.var("NHL")
a.nhl<-a.var[[1]]
t.nhl<-a.var[[2]]
a.var<-make.ade.var("OTH")
a.oth<-a.var[[1]]
t.oth<-a.var[[2]]
a.var<-make.ade.var("PCP")
a.pcp<-a.var[[1]]
t.pcp<-a.var[[2]]
a.var<-make.ade.var("PML")
a.pml<-a.var[[1]]
t.pml<-a.var[[2]]
a.var<-make.ade.var("RET")
a.ret<-a.var[[1]]
t.ret<-a.var[[2]]
a.var<-make.ade.var("SAL")
a.sal<-a.var[[1]]
t.sal<-a.var[[2]]
a.var<-make.ade.var("SPO")
a.spo<-a.var[[1]]
t.spo<-a.var[[2]]
a.var<-make.ade.var("TOX")
a.tox<-a.var[[1]]
t.tox<-a.var[[2]]
a.var<-make.ade.var("TUB")
a.tub<-a.var[[1]]
t.tub<-a.var[[2]]
a.var<-make.ade.var("UNK")
a.unk<-a.var[[1]]
t.unk<-a.var[[2]]
a.var<-make.ade.var("WAS")
a.was<-a.var[[1]]
t.was<-a.var[[2]]


a.ade<-apply(cbind(a.bpn,a.can,a.ccm,a.cmv,a.coc,a.dem,a.hgl,a.his,a.hsd,a.icc,a.ile,a.iso,a.ksa,a.lei,a.lob,a.mac,a.myc,
                   a.nhl,a.oth,a.pcp,a.pml,a.ret,a.sal,a.spo,a.tox,a.tub,a.unk,a.was),
             1,sum)
t.ade<-apply(cbind(t.bpn,t.can,t.ccm,t.cmv,t.coc,t.dem,t.hgl,t.his,t.hsd,t.icc,t.ile,t.iso,t.ksa,t.lei,t.lob,t.mac,t.myc,
                   t.nhl,t.oth,t.pcp,t.pml,t.ret,t.sal,t.spo,t.tox,t.tub,t.unk,t.was),
             1,min)


####  Combining any event that doesn't have a total number of 100 in 'other'
a.other<-apply(cbind(a.ccm,a.coc,a.hgl,a.his,a.icc,a.ile,a.iso,a.lei,a.lob,
                     a.oth,a.ret,a.sal,a.spo,a.unk,a.was),
               1,sum)
t.other<-apply(cbind(t.ccm,t.coc,t.hgl,t.his,t.icc,t.ile,t.iso,t.lei,t.lob,
                     t.oth,t.ret,t.sal,t.spo,t.unk,t.was),
               1,min)


####  Combining any event that doesn't have a total number of 60 in 'other' (this was chosen because there is a big gap from 63 to 35)
a.other1<-apply(cbind(a.ccm,a.hgl,a.his,a.icc,a.ile,a.iso,a.lei,a.lob,
                      a.oth,a.sal,a.spo,a.unk),
                1,sum)
t.other1<-apply(cbind(t.ccm,t.hgl,t.his,t.icc,t.ile,t.iso,t.lei,t.lob,
                      t.oth,t.sal,t.spo,t.unk),
                1,min)



nb<-d$native.born.factor

t<-as.numeric(d$incidence_end_of_fup)-as.numeric(d$haart_d)

sum(a.ade[nb=="Yes"])
sum(a.ade[nb=="No"])


library(survival)
library(rms)

inc.calc<-function(event,time,fup,nb,t) {
  inc.nb1<-sum(event[nb=="Yes"])/sum(fup[nb=="Yes"])*1000
  inc.nb0<-sum(event[nb=="No"])/sum(fup[nb=="No"])*1000
  
  mod0<-glm(event~1+offset(log(t/365.25)),family=poisson,subset=nb=="No")
  mod1<-glm(event~1+offset(log(t/365.25)),family=poisson,subset=nb=="Yes")
  
  t1<-ifelse(time<t,time,t)  # for Cox  regression
  e1<-ifelse(event>1,1,event)
  
  mod.cph<-coxph(Surv(t1,e1)~nb+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                   +rcs(d$age,3)+d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$risk.factor+d$stagec.factor+strata(d$cohort))
  
  
  
  list(total=sum(event),
       inc.natives=c(inc.nb1, exp(mod1$coeff[1])*1000,
                     exp(mod1$coeff[1]-1.96*summary(mod1)$coeff[1,2])*1000, exp(mod1$coeff[1]+1.96*summary(mod1)$coeff[1,2])*1000),
       inc.nonnatives=c(inc.nb0, exp(mod0$coeff[1])*1000,
                        exp(mod0$coeff[1]-1.96*summary(mod0)$coeff[1,2])*1000, exp(mod0$coeff[1]+1.96*summary(mod0)$coeff[1,2])*1000),
       hr.ref.natives=c(summary(mod.cph)$coeff[1,2],summary(mod.cph)$conf.int[1,3],summary(mod.cph)$conf.int[1,4])
  )
  
  
}


i.bpn<-inc.calc(a.bpn,t.bpn,fup,nb,t)
i.can<-inc.calc(a.can,t.can,fup,nb,t)
i.ccm<-inc.calc(a.ccm,t.ccm,fup,nb,t)
i.cmv<-inc.calc(a.cmv,t.cmv,fup,nb,t)
i.coc<-inc.calc(a.coc,t.coc,fup,nb,t)
i.dem<-inc.calc(a.dem,t.dem,fup,nb,t)
i.hgl<-inc.calc(a.hgl,t.hgl,fup,nb,t)
i.his<-inc.calc(a.his,t.his,fup,nb,t)
i.hsd<-inc.calc(a.hsd,t.hsd,fup,nb,t)
i.icc<-inc.calc(a.icc,t.icc,fup,nb,t)
i.ile<-inc.calc(a.ile,t.ile,fup,nb,t)
i.iso<-inc.calc(a.iso,t.iso,fup,nb,t)
i.ksa<-inc.calc(a.ksa,t.ksa,fup,nb,t)
i.lei<-inc.calc(a.lei,t.lei,fup,nb,t)
i.lob<-inc.calc(a.lob,t.lob,fup,nb,t)
i.mac<-inc.calc(a.mac,t.mac,fup,nb,t)
i.myc<-inc.calc(a.myc,t.myc,fup,nb,t)
i.nhl<-inc.calc(a.nhl,t.nhl,fup,nb,t)
i.oth<-inc.calc(a.oth,t.oth,fup,nb,t)
i.pcp<-inc.calc(a.pcp,t.pcp,fup,nb,t)
i.pml<-inc.calc(a.pml,t.pml,fup,nb,t)
i.ret<-inc.calc(a.ret,t.ret,fup,nb,t)
i.sal<-inc.calc(a.sal,t.sal,fup,nb,t)
i.spo<-inc.calc(a.spo,t.spo,fup,nb,t)
i.tox<-inc.calc(a.tox,t.tox,fup,nb,t)
i.tub<-inc.calc(a.tub,t.tub,fup,nb,t)
i.unk<-inc.calc(a.unk,t.unk,fup,nb,t)
i.was<-inc.calc(a.was,t.was,fup,nb,t)

i.ade<-inc.calc(a.ade,t.ade,fup,nb,t)
i.other<-inc.calc(a.other,t.other,fup,nb,t)
i.other1<-inc.calc(a.other1,t.other1,fup,nb,t)


sum(a.ade)/sum(fup)*1000   ## overall incidence of ADE  (used in paper)
mod.overall<-glm(a.ade~1+offset(log(t/365.25)),family=poisson)
c(exp(mod.overall$coeff[1])*1000,
  exp(mod.overall$coeff[1]-1.96*summary(mod.overall)$coeff[1,2])*1000,
  exp(mod.overall$coeff[1]+1.96*summary(mod.overall)$coeff[1,2])*1000)


i.ade
i.tub
i.ksa
i.can
i.mac
i.cmv
i.nhl
i.tox
i.pcp
i.bpn
i.hsd
i.pml
i.dem
i.myc
i.ret
i.coc
i.was
i.other1

i.spo
i.hgl
i.lei
i.icc #warning
i.iso #warning
i.his #warning
i.sal #warning
i.unk
i.ile #warning, no event
i.ccm #warning, no event
i.oth #warning, no event
i.ccm # #warning, no event

# Decided to use other 'other' (other1 not other)
i.other


table(d$artcc_origin_g.factor)

## incidence of ADE by region  (used in paper)
mod.ssa<-glm(a.ade~1+offset(log(t/365.25)),family=poisson,subset=d$artcc_origin_g.factor=="SSA")
c(exp(mod.ssa$coeff[1])*1000,
  exp(mod.ssa$coeff[1]-1.96*summary(mod.ssa)$coeff[1,2])*1000,
  exp(mod.ssa$coeff[1]+1.96*summary(mod.ssa)$coeff[1,2])*1000)

mod.la<-glm(a.ade~1+offset(log(t/365.25)),family=poisson,subset=d$artcc_origin_g.factor=="LA")
c(exp(mod.la$coeff[1])*1000,
  exp(mod.la$coeff[1]-1.96*summary(mod.la)$coeff[1,2])*1000,
  exp(mod.la$coeff[1]+1.96*summary(mod.la)$coeff[1,2])*1000)

mod.name<-glm(a.ade~1+offset(log(t/365.25)),family=poisson,subset=d$artcc_origin_g.factor=="NAME")
c(exp(mod.name$coeff[1])*1000,
  exp(mod.name$coeff[1]-1.96*summary(mod.name)$coeff[1,2])*1000,
  exp(mod.name$coeff[1]+1.96*summary(mod.name)$coeff[1,2])*1000)

mod.asia<-glm(a.ade~1+offset(log(t/365.25)),family=poisson,subset=d$artcc_origin_g.factor=="ASIA")
c(exp(mod.asia$coeff[1])*1000,
  exp(mod.asia$coeff[1]-1.96*summary(mod.asia)$coeff[1,2])*1000,
  exp(mod.asia$coeff[1]+1.96*summary(mod.asia)$coeff[1,2])*1000)


#####  Looking closer at the mortality analysis -- results are largely driven by age

mod.death<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                   +rcs(d$age,3)+d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$risk.factor+d$stagec.factor+strata(d$cohort))

mod.death.no.age<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                          +d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$risk.factor+d$stagec.factor+strata(d$cohort))

mod.death.only.age<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+rcs(d$age,3)+strata(d$cohort))


mod.death.no.risk.factor<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                                  +d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$stagec.factor+strata(d$cohort))

mod.death.hetero<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                          +rcs(d$age,3)+d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$stagec.factor+strata(d$cohort),subset=d$risk.factor=="Heterosexual")

mod.death.no.idu<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                          +rcs(d$age,3)+d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$stagec.factor+strata(d$cohort),subset=d$risk.factor!="IDU")

mod.death.hetero.unadj<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+strata(d$cohort),subset=d$risk.factor=="Heterosexual")

mod.death.no.idu.unadj<-coxph(Surv(d$followup_time_first_year,d$death_first_year_of_haart)~d$native.born.factor+strata(d$cohort),subset=d$risk.factor!="IDU")



boxplot(d$age~d$native.born.factor)
wilcox.test(d$age~d$native.born.factor)

boxplot(sqrt(d$cd4_t0_corrected)~d$native.born.factor)
wilcox.test(d$cd4_t0_corrected~d$native.born.factor)




##### Cox regression for TB and specific regions of origin
t1<-ifelse(t.tub<t,t.tub,t) 
e1<-ifelse(a.tub>1,1,a.tub)
mod.go.tb<-coxph(Surv(t,e1)~d$artcc_origin_g.factor+d$gender.factor+rcs(d$cd4_t0_corrected,4)+rcs(d$lrna_t0_corrected,3)+
                   +rcs(d$age,3)+d$pi_based.factor+d$nnrti_based.factor+rcs(d$haartyr,3)+d$risk.factor+d$stagec.factor+strata(d$cohort))


round(c(summary(mod.go.tb)$conf.int[1,1],summary(mod.go.tb)$conf.int[1,3],summary(mod.go.tb)$conf.int[1,4]),2)  # NAME
round(c(summary(mod.go.tb)$conf.int[2,1],summary(mod.go.tb)$conf.int[2,3],summary(mod.go.tb)$conf.int[2,4]),2)  # SSA
round(c(summary(mod.go.tb)$conf.int[3,1],summary(mod.go.tb)$conf.int[3,3],summary(mod.go.tb)$conf.int[3,4]),2)  # ASIA
round(c(summary(mod.go.tb)$conf.int[4,1],summary(mod.go.tb)$conf.int[4,3],summary(mod.go.tb)$conf.int[4,4]),2)  # LA







#pdf("ade-HR.pdf",height=7,width=8)
postscript("figure3-paper.eps",height=7, width=8)
par(mar=c(4.2,.2,2,.2))
plot.lines<-function(x,loc){
  lines(log(c(x[2],x[3])),c(loc,loc),lwd=1.5)
  points(log(x[1]),loc,pch=20,cex=1.5)
  text(k1-k2,loc,paste(round(x[1],2)," (",paste(round(x[2],2),round(x[3],2),sep=", "),")",sep=""),pos=4,cex=a)
}


k<- -2.2
k1<-1.8
k2<-.1
a<-.8

plot(c(k-1.3, 1+k1),c(12,31),type="n",xlab="",ylab="",bty="n",axes=FALSE)
abline(v=0,col=gray(.8))


k3<-30

plot.lines(i.ade$hr.ref.natives,k3+1)
plot.lines(i.tub$hr.ref.natives,k3-1)
plot.lines(i.ksa$hr.ref.natives,k3-2)
plot.lines(i.can$hr.ref.natives,k3-3)
plot.lines(i.mac$hr.ref.natives,k3-4)
plot.lines(i.cmv$hr.ref.natives,k3-5)
plot.lines(i.nhl$hr.ref.natives,k3-6)
plot.lines(i.tox$hr.ref.natives,k3-7)
plot.lines(i.pcp$hr.ref.natives,k3-8)
plot.lines(i.bpn$hr.ref.natives,k3-9)
plot.lines(i.hsd$hr.ref.natives,k3-10)
plot.lines(i.pml$hr.ref.natives,k3-11)
plot.lines(i.dem$hr.ref.natives,k3-12)
plot.lines(i.myc$hr.ref.natives,k3-13)
plot.lines(i.ret$hr.ref.natives,k3-14)
plot.lines(i.coc$hr.ref.natives,k3-15)
plot.lines(i.was$hr.ref.natives,k3-16)
plot.lines(i.other1$hr.ref.natives,k3-17)

text(k-1.4,k3+1,pos=4,cex=a,paste("Any ADE  (n=", i.ade$total, ")",sep=""))
text(k-1.4,k3-1,pos=4,cex=a,paste("Tuberculosis  (", i.tub$total, ")",sep=""))
text(k-1.4,k3-2,pos=4,cex=a,paste("Kaposi's Sarcoma  (", i.ksa$total, ")",sep=""))
text(k-1.4,k3-3,pos=4,cex=a,paste("Candidiasis  (", i.can$total, ")",sep=""))
text(k-1.4,k3-4,pos=4,cex=a,paste("MAC or Kansasii  (", i.mac$total, ")",sep=""))
text(k-1.4,k3-5,pos=4,cex=a,paste("CMV  (", i.cmv$total, ")",sep=""))
text(k-1.4,k3-6,pos=4,cex=a,paste("Non-Hodgkin  (", i.nhl$total, ")",sep=""))
text(k-1.4,k3-7,pos=4,cex=a,paste("Toxoplasmosis  (", i.tox$total, ")",sep=""))
text(k-1.4,k3-8,pos=4,cex=a,paste("PCP  (", i.pcp$total, ")",sep=""))
text(k-1.4,k3-9,pos=4,cex=a,paste("Bacterial pneumonia  (", i.bpn$total, ")",sep=""))
text(k-1.4,k3-10,pos=4,cex=a,paste("Herpes simplex virus ulcers  (", i.hsd$total, ")",sep=""))
text(k-1.4,k3-11,pos=4,cex=a,paste("Leucoencephalopathy  (", i.pml$total, ")",sep=""))
text(k-1.4,k3-12,pos=4,cex=a,paste("AIDS dementia  (", i.dem$total, ")",sep=""))
text(k-1.4,k3-13,pos=4,cex=a,paste("Mycobacterial extrapulm  (", i.myc$total, ")",sep=""))
text(k-1.4,k3-14,pos=4,cex=a,paste("CMV chorioretinitis  (", i.ret$total, ")",sep=""))
text(k-1.4,k3-15,pos=4,cex=a,paste("Cryptococcosis  (", i.coc$total, ")",sep=""))
text(k-1.4,k3-16,pos=4,cex=a,paste("HIV wasting syndrome  (", i.was$total, ")",sep=""))
text(k-1.4,k3-17,pos=4,cex=a,paste("Other  (", i.other$total, ")",sep=""))


axis(1,at=log(c(.25,.33,.5,.67,1,1.5,2,3,4)),labels=c(.25,.33,.5,.67,1,1.5,2,3,4),cex.axis=.7)

mtext("Adjusted Hazard Ratio", side=1, at=0, line=3, cex=a)

dev.off()



#####  p-value


chisq.test(table(nb,d$european.clinic))

with(d, chisq.test(table(as.character(origin_g_combined.factor),european.clinic)))



###### COD stuff

death1<-d$death_first_year_of_haart
ade1<-d$ade_first_year_of_haart
cod<-d$code
table(death1,ade1,nb)

table(cod[death1==1&nb=="No"])
table(cod[death1==1&nb=="Yes"])

sum(table(cod[death1==1&nb=="No"]))- sum(cod[death1==1&nb=="No"]=="920. Unknown")
sum(table(cod[death1==1&nb=="Yes"]))- sum(cod[death1==1&nb=="Yes"]=="920. Unknown",na.rm=TRUE)

sum(death1[nb=="No"])
sum(death1[nb=="Yes"])


table(cod[death1==1&nb=="No"])/(sum(table(cod[death1==1&nb=="No"]))- sum(cod[death1==1&nb=="No"]=="920. Unknown"))
table(cod[death1==1&nb=="Yes"])/(sum(table(cod[death1==1&nb=="Yes"]))- sum(cod[death1==1&nb=="Yes"]=="920. Unknown",na.rm=TRUE))


## number with an AIDS-related COD

sum(cod[death1==1&nb=="No"]=="10. AIDS"|cod[death1==1&nb=="No"]=="11. AIDS infection"|cod[death1==1&nb=="No"]=="12. AIDS malignancy")
sum(cod[death1==1&nb=="Yes"]=="10. AIDS"|cod[death1==1&nb=="Yes"]=="11. AIDS infection"|cod[death1==1&nb=="Yes"]=="12. AIDS malignancy",na.rm=TRUE)


## Of those who died in the first year with known cause of death, the percentage that was AIDS-related

sum(cod[death1==1&nb=="No"]=="10. AIDS"|cod[death1==1&nb=="No"]=="11. AIDS infection"|cod[death1==1&nb=="No"]=="12. AIDS malignancy")/
  (sum(table(cod[death1==1&nb=="No"]))- sum(cod[death1==1&nb=="No"]=="920. Unknown"))
sum(cod[death1==1&nb=="Yes"]=="10. AIDS"|cod[death1==1&nb=="Yes"]=="11. AIDS infection"|cod[death1==1&nb=="Yes"]=="12. AIDS malignancy",na.rm=TRUE)/
  (sum(table(cod[death1==1&nb=="Yes"]))- sum(cod[death1==1&nb=="Yes"]=="920. Unknown",na.rm=TRUE))

sum(cod[death1==1]=="10. AIDS"|cod[death1==1]=="11. AIDS infection"|cod[death1==1]=="12. AIDS malignancy",na.rm=TRUE)/
  (sum(table(cod[death1==1]))- sum(cod[death1==1]=="920. Unknown",na.rm=TRUE))

sum(death1)-(sum(table(cod[death1==1]))- sum(cod[death1==1]=="920. Unknown",na.rm=TRUE))

##

sum(death1[ade1==1&nb=="No"])
sum(ade1==1&nb=="No")
mean(death1[ade1==1&nb=="No"])

sum(death1[ade1==1&nb=="Yes"])
sum(ade1==1&nb=="Yes")
mean(death1[ade1==1&nb=="Yes"])

sum(death1[ade1==1])
sum(ade1)

sum(ade1[death1==1])
sum(death1)

sum(ade1[death1==1&nb=="Yes"])
sum(death1&nb=="Yes")

sum(ade1[death1==1&nb=="No"])
sum(death1&nb=="No")


table(cod[death1==1&ade1==1&nb=="Yes"])
table(cod[death1==1&ade1==1&nb=="No"])





sum(ade1[death1==1]==1|cod[death1==1]=="10. AIDS"|cod[death1==1]=="11. AIDS infection"|cod[death1==1]=="12. AIDS malignancy",na.rm=TRUE)
sum(death1)

sum(ade1[death1==1&nb=="Yes"]==1|cod[death1==1&nb=="Yes"]=="10. AIDS"|cod[death1==1&nb=="Yes"]=="11. AIDS infection"|cod[death1==1&nb=="Yes"]=="12. AIDS malignancy",na.rm=TRUE)
sum(death1&nb=="Yes")

sum(ade1[death1==1&nb=="No"]==1|cod[death1==1&nb=="No"]=="10. AIDS"|cod[death1==1&nb=="No"]=="11. AIDS infection"|cod[death1==1&nb=="No"]=="12. AIDS malignancy",na.rm=TRUE)
sum(death1&nb=="No")



sum(cod[death1==1&nb=="No"&ade1==1]=="10. AIDS"|cod[death1==1&nb=="No"&ade1==1]=="11. AIDS infection"|cod[death1==1&nb=="No"&ade1==1]=="12. AIDS malignancy")/
  (sum(table(cod[death1==1&nb=="No"&ade1==1]))- sum(cod[death1==1&nb=="No"&ade1==1]=="920. Unknown"))

sum(cod[death1==1&nb=="Yes"&ade1==1]=="10. AIDS"|cod[death1==1&nb=="Yes"&ade1==1]=="11. AIDS infection"|cod[death1==1&nb=="Yes"&ade1==1]=="12. AIDS malignancy")/
  (sum(table(cod[death1==1&nb=="Yes"&ade1==1]))- sum(cod[death1==1&nb=="Yes"&ade1==1]=="920. Unknown"))

sum(cod[death1==1&ade1==1]=="10. AIDS"|cod[death1==1&ade1==1]=="11. AIDS infection"|cod[death1==1&ade1==1]=="12. AIDS malignancy")/
  (sum(table(cod[death1==1&ade1==1]))- sum(cod[death1==1&ade1==1]=="920. Unknown"))




#######

table(a.tub>0,death1,nb)