#2103.10.30

setwd("/Users/jameshanley/Dropbox/Courses/bios601/Epidemiology")

# usa data on

# deaths

d=read.table("usa-d.txt",head=T)
head(d) ; tail(d)

# population 

p=read.table("usa-p.txt",head=T)
head(p) ; tail(p)


D=(aggregate(d$d,by=list(d$a),FUN="sum")$x)[1:105]
P=(aggregate(p$p,by=list(p$a),FUN="sum")$x)[1:105]

death.rate = D/P ; summary(death.rate) ; 

plot(1:105,death.rate)
plot(1:105,log(appendicitis.rate))

# Michelle T. Buckius, 2010 recent usa rates of appendectomy 2005-2008
# rate per 10,000 person years, by decade of age

appendicitis.rate = c(5.9,15.3,13.0,10.8,7.8,7.3,6.7,5.5,4.0,2.2,0.7)/10000
appendicitis.rate= rep(appendicitis.rate,each=10)[1:105]
plot(1:105,appendicitis.rate)


Alive.Appendix.Intact  = 100*exp(- 
                cumsum(appendicitis.rate*1 + death.rate*1 ) )


YLIM=c(95,100); XLIM=c(0,25)
#YLIM=c(90,100); XLIM=c(0,100)


par(mfrow=c(1,1),mar = c(4.5,4,0,0))

plot(0:100,c(100,Alive.Appendix.Intact[1:100]),
    ylim=YLIM,xlim=XLIM+c(0,1),cex=0.5,pch=19,col="white",xlab="Age",ylab="Number")

n.appendices.removed = Alive.Appendix.Intact *
                       appendicitis.rate
                       
cum.n.removed = cumsum(n.appendices.removed)
            
n.deaths =             Alive.Appendix.Intact * 
                       death.rate

sum(n.appendices.removed)

segments(0,100,0,0,col="green", lwd=3 )

text(-0.7,100,"Alive, Intact Appendix",col="green", lwd=3, srt=90, adj=c(1,1) )
 
text(mean(XLIM), 99.9,"Appendices removed", col="blue", adj=c(0.5,1))


for(a in 1:XLIM[2]) {
	
 segments(
   a-0.5,Alive.Appendix.Intact[a],
   a-0.5,Alive.Appendix.Intact[a]+n.deaths[a],
   col="red", lwd=3, lend=1
 )
 if(a==1) text(0.6, Alive.Appendix.Intact[a]+n.deaths[a]/2,
     "Deaths", srt=42, col="red", adj=c(0,0.5))
     
 if(a==XLIM[2]) {
 	
 	text(a+1.5, Alive.Appendix.Intact[a]+n.deaths[a]/2+0.25,
     toString(round(cumsum(n.appendices.removed)[a],1)), col="blue",
     adj=c(1,-0.25))
     
 	text(a+1.5, Alive.Appendix.Intact[a]+n.deaths[a]/2,
     toString(round(cumsum(n.deaths)[a],1)), col="red",
     adj=c(1,-0.25))
     
     text(a+1.5, Alive.Appendix.Intact[a],
     toString(round(Alive.Appendix.Intact[a],1)), col="green",
     adj=c(1,1.25))


 	
 }
   
 segments(
   a-0.5,Alive.Appendix.Intact[a]+n.deaths[a],
   a-0.5,Alive.Appendix.Intact[a]+n.deaths[a]+n.appendices.removed[a],
   col="blue", lwd=3, lend=1
 )
 text(a-0.5,
      Alive.Appendix.Intact[a]+n.deaths[a]+n.appendices.removed[a],
       paste(toString(round(n.appendices.removed[a],2) ),
             " cum:",
             toString(round(cum.n.removed[a],1) )
            ),
       srt=45,cex=0.75,col="blue", adj=c(0,0)
     )
      
 segments(
   a,Alive.Appendix.Intact[a],
   a,0,
   col="green", lwd=3, lend=1
 )

 
}