rm(list=ls(all=TRUE)) filecanc=c('breast','digothr','malegen','femgen','other','respir','colrect','lymyleuk','urinary') load("/data/SEER/pops.RData") # this loads in 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)) age=c(0.5,3,seq(7.5,87.5,5)) graphics.off() windows(height=5,width=12) par(mfrow=c(1,2),mar=c(4.7,6.2,0.8,1),lwd=3,cex.lab=2,cex.axis=2,cex.main=1.3) for (k in c(8:8)) { load(paste("/data/SEER/",filecanc[k],".RData",sep="")) # this loads in DF DF=DF[!is.na(DF$ICD9),] # get rid of missing ICD9 entries if (k==1) { canc=c('Breast'); code=list(c(1741,1748)) } if (k==2) {canc=c('Esophogeal','Stomach','Liver','Pancreatic'); code=list(c(1504,1505),c(1510,1519),c(1550,1550),c(1570,1579)) } if (k==3) {canc=c("Prostate"); code=list(c(185,185))} if (k==4) {canc=c("Cervical","Uterine","Vaginal"); code=list(c(1809,1809),c(1820,1820),c(1830,1830))} if (k==5) {canc=c("Thyroid","Brain"); code=list(c(193,193),c(1911,1919))} if (k==6){canc=c("Laryngeal","Lung");code=list(c(1610,1610),c(1622,1629))} if (k==7){canc=c("Colorectal"); code=list(c(1530,1539)) } if (k==8) { canc=c("NHL", "Hodgkins", "MML" ,"ALL"," CLL", "AML", "CML") code=list(c(2000,2008),c(2015,2016),c(2030,2030),c(2040,2040),c(2041,2041),c(2050,2050),c(2051,2051))} if (k==9) { canc=c("Bladder","Kidney"); code=list(c(1880,1889),c(1890,1891))} for (i in 1:length(code)) { cat("k =",k," i =",i,"\n") Indx=((DF$ICD9>=code[[i]][1])&(DF$ICD9<=code[[i]][2])&(DF$numPrim==1)) # &(DF$numPrim==1) ) ndat=DF[Indx,] ndat1=subset(ndat,yDiag<1985) ndat2=subset(ndat,(yDiag>=1985)&(yDiag<1997)) ndat3=subset(ndat,yDiag>=1997) if ((k!=1)&(k!=4)) # if not breast and not fem gen i.e. if males { m1=with(ndat1,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) m2=with(ndat2,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) m3=with(ndat3,hist(ageRC[sex==1],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) plot(age,m1$counts/py1$m,log="y",type='b',xlab="Age",pch=1,col='red',ylab="",yaxt="n", ylim=c(4.9e-7,6e-2),cex=1.3) lines(age,m2$counts/py2$m,type='b',col='blue',pch=2,cex=1.3) lines(age,m3$counts/py3$m,type='b',col='black',pch=3,cex=1.3) # axis(side=2,las=1, at=c(1e-6,1e-5,1e-4),labels=expression(1,10,10^2)) axis(side=2,las=1, at=c(1e-7,1e-6,1e-5,1e-4,1e-3,1e-2),labels=expression(.1,1,10,10^2,10^3,10^4)) mtext(expression(paste("Cases per ",10^6," PY")),side=2,line=3.7,cex=2) mtext("Males",side=3,line=-1.8,cex=2,adj=0.05) mtext(canc[i],side=3,line=-1.8,cex=2,adj=0.95) legend("bottomright",c("1973-1984","1985-1996","1997-2009"), col=c("red","blue","black"),pch=1:3,bty="n",cex=1.5) } else plot.new() if (k!=3) # not males only i.e. if females { f1=with(ndat1,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) f2=with(ndat2,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) f3=with(ndat3,hist(ageRC[sex==2],breaks=c(seq(-.5,17.5,1),100),plot=FALSE)) plot(age,f1$counts/py1$f,log="y",type='b',xlab="Age",pch=1,cex=1.3,col='red',ylab="",yaxt="n", ylim=c(4.9e-7,6e-2) ) lines(age,f2$counts/py2$f,type='b',col='blue',pch=2,cex=1.3) lines(age,f3$counts/py3$f,type='b',col='black',pch=3,cex=1.3) mtext("Females",side=3,line=-1.8,cex=2,adj=0.05) mtext(canc[i],side=3,line=-1.8,cex=2,adj=0.95) # axis(side=2,las=1, at=c(1e-6,1e-5,1e-4),labels=expression(1,10,10^2)) axis(side=2,las=1, at=c(1e-7,1e-6,1e-5,1e-4,1e-3,1e-2),labels=expression(.1,1,10,10^2,10^3,10^4)) mtext(expression(paste("Cases per ",10^6," PY")),side=2,line=3.5,cex=2) legend("bottomright",c("1973-1984","1985-1996","1997-2009"), col=c("red","blue","black"),pch=1:3,bty="n",cex=1.5) } else plot.new() locator(n=1) } # the i for loop over within file cancer types } # k loop on files