rm(list=ls(all=TRUE)) cols=c("city","sex","doseg","agexg","calg","kerma","PY","adjPY","num.entering", "age","agex","tsx","cal","sv","gam","neut","lymphoma","NHL","leukemia","AML","ALL","CML","ATL","MML") d<-read.table("/data/abomb/HEMA87.DAT", header=F,col.names=cols); d=d[d$adjPY>0,] #remove any recs with zero py d$adjPY=d$adjPY*1e4 # convert to real person years d=d[d$kerma==1,] # take only kerma < 4 Gy CITY=1 if (CITY>0) d=d[d$city==CITY,]#hiroshima = 1 d=d[d$sex==1,]#males attach(d); #summary(d) svc=cut(sv,c(-1,.02,1,10)) head(d) #sort(unique(agex[CML>0])) (PYT=tapply(adjPY,list(svc,agexg),sum,na.rm=TRUE)) (AGET=tapply(agex,list(svc,agexg),mean,na.rm=TRUE)) (CMLT=tapply(CML,list(svc,agexg),sum,na.rm=TRUE)) (CMLt=cbind(apply(CMLT[,1:4],1,sum),apply(CMLT[,5:8],1,sum),apply(CMLT[,9:13],1,sum))) (PYt=cbind(apply(PYT[,1:4],1,sum),apply(PYT[,5:8],1,sum),apply(PYT[,9:13],1,sum))) (AGEt=cbind(apply(AGET[,1:4],1,mean),apply(AGET[,5:8],1,mean),apply(AGET[,9:13],1,mean))) (AGEw=cbind(apply(AGET[,1:4]*PYT[,1:4],1,sum),apply(AGET[,5:8]*PYT[,5:8],1,sum),apply(AGET[,9:13]*PYT[,9:13],1,sum))) (AGEtw=AGEw/PYt) (incid=CMLt/PYt) #detach(d) sum(PYt) brb=c("blue","red","black") pchs=c(2,1,22) #graphics.off() windows(width=10,height=5) par(mfrow=c(1,2),mar=c(4.7,0,0.5,0.9),lwd=3,cex.lab=1.8,cex.axis=2,cex.main=1.7,oma=c(0,5.5,1.4,0.7)) matplot(t(AGEtw),t(incid),log="y",type='b',ylab="", xlab="age at exposure", ylim=c(1e-6,15e-4),lty=1,col=brb,cex=2.5,pch=pchs,lwd=3,cex.lab=2,yaxt="n") axis(2,at=c(1e-6,1e-5,1e-4,1e-3),labels=c(0.1,1,10,100),las=1,cex.axis=1.7,outer=T) mtext(expression(paste("Cases per ",10^5," Person-Year-Sv")),side=2,line=3.3,cex=1.8,outer=T) mtext("Males",side=3,line=-2,cex=1.8,adj=0.05) text(30,25e-5,"mostly radiogenic",cex=1.5) legend(20,1.3e-5,c("High Dose","Medium Dose","Low Dose"),pch=pchs[3:1],col=brb[3:1],cex=1.5,pt.cex=2,pt.lwd=2,bty="n") title("Japanese A-bomb Survivors",outer=T) detach(d) d<-read.table("/data/abomb/HEMA87.DAT", header=F,col.names=cols); d=d[d$adjPY>0,] #remove any recs with zero py d$adjPY=d$adjPY*1e4 # convert to real person years d=d[d$kerma==1,] # take only kerma < 4 Gy if (CITY>0) d=d[d$city==CITY,]#hiroshima d=d[d$sex==2,]#2=females attach(d); #summary(d) svc=cut(sv,c(-1,.02,1,10)) head(d) #sort(unique(agex[CML>0])) (PYT=tapply(adjPY,list(svc,agexg),sum,na.rm=TRUE)) (AGET=tapply(agex,list(svc,agexg),mean,na.rm=TRUE)) (CMLT=tapply(CML,list(svc,agexg),sum,na.rm=TRUE)) (CMLt=cbind(apply(CMLT[,1:4],1,sum),apply(CMLT[,5:8],1,sum),apply(CMLT[,9:13],1,sum))) (PYt=cbind(apply(PYT[,1:4],1,sum),apply(PYT[,5:8],1,sum),apply(PYT[,9:13],1,sum))) (AGEt=cbind(apply(AGET[,1:4],1,mean),apply(AGET[,5:8],1,mean),apply(AGET[,9:13],1,mean))) (AGEw=cbind(apply(AGET[,1:4]*PYT[,1:4],1,sum),apply(AGET[,5:8]*PYT[,5:8],1,sum),apply(AGET[,9:13]*PYT[,9:13],1,sum))) (AGEtw=AGEw/PYt) (incid=CMLt/PYt) par(mar=c(4.7,1,0.5,0)) matplot(t(AGEtw),t(incid),log="y",type='b',ylab="", xlab="age at exposure", ylim=c(1e-6,15e-4),lty=1,col=brb,cex=2.5,pch=pchs,lwd=3,cex.lab=2,yaxt="n") text(30,25e-5,"mostly radiogenic",cex=1.5) mtext("Females",side=3,line=-2,cex=1.8,adj=0.05) detach(d)