filecanc=c('breast','digothr','malegen','femgen','other','respir','colrect','lymyleuk','urinary') load("/data/SEER/pops.RData") # this loads in pops head(pops) pops1=subset(pops,popyear<1985) pops2=subset(pops,(popyear>=1985)&(popyear<1997)) pops3=subset(pops,popyear>=1997) getSexPY<-function(pops) { attach(pops) pym=NULL;pyf=NULL for (i in 0:18) { pym[i+1]=sum(population[(popsex==1)&(popage==i)]) pyf[i+1]=sum(population[(popsex==2)&(popage==i)]) } detach(pops) list(m=pym,f=pyf) } (py1=getSexPY(pops1)) (py2=getSexPY(pops2)) (py3=getSexPY(pops3)) M1=F1=numeric(19) M2=F2=numeric(19) M3=F3=numeric(19) graphics.off() par(mfcol=c(3,3),mar=c(5,4,1,1)+.1) for (k in c(1:9)) { load(paste("/data/SEER/",filecanc[k],".RData",sep="")) # this loads in DF ndat1=subset(DF,yDiag<1985) ndat2=subset(DF,(yDiag>=1985)&(yDiag<1997)) ndat3=subset(DF,yDiag>=1997) if ((k!=1)&(k!=4)) # if not breast and not fem gen { m1=with(ndat1,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) m2=with(ndat2,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) m3=with(ndat3,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) M1=M1+m1 M2=M3+m2 M3=M3+m3 } if (k!=3) # if not male gen { f1=with(ndat1,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) f2=with(ndat2,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) f3=with(ndat3,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)$counts) F1=F1+f1 F2=F2+f2 F3=F3+f3 } } # k loop on files age=c(0.5,3,seq(7.5,87.5,5)) graphics.off() windows(height=7,width=12) par(mfrow=c(1,2),mar=c(4.5,5.2,0.8,1),lwd=3,cex.lab=1.5,cex=1.3,cex.main=1.3,oma=c(0,0,2,0)) plot(age,M1/py1$m,log="y",type='b',xlab="Age",pch=1,col='blue',ylab="Cases per Person-Year",yaxt="n",ylim=c(1e-4,5e-2),cex=1.3) lines(age,M2/py2$m,type='b',col='red',pch=2) lines(age,M3/py3$m,type='b',col='black',pch=3) axis(side=2,las=1, at=c(1e-4,1e-3,1e-2),labels=expression(10^-4,10^-3,10^-2)) mtext("Males",side=3,line=-1.8,cex=2,adj=0.05) legend("bottomright",c("1973-1984","1985-1996","1997-2008"), col=c("blue","red","black"),pch=1:3,bty="n",cex=1.5) plot(age,F1/py1$f,log="y",type='b',xlab="Age",pch=1,col='blue', ylab="Cases per Person-Year",yaxt="n",ylim=c(1e-4,5e-2) ) lines(age,F2/py2$f,type='b',col='red',pch=2) lines(age,F3/py3$f,type='b',col='black',pch=3) mtext("Females",side=3,line=-1.8,cex=2,adj=0.05) axis(side=2,las=1, at=c(1e-4,1e-3,1e-2),labels=expression(10^-4,10^-3,10^-2)) title(main="Incidence of All Cancers",outer=T)