#pre='yr1973_2002.seer9/' # this one is a bit slow on breast and prostate pre='yr1992_2002.sj_la_rg_ak/' # comment if you want all years filecanc=c('breast','digothr','malegen','femgen','other','respir','colrect','lymyleuk','urinary') load(paste("c:/SEER02/",pre,"pops.RData",sep="")) # this loads in 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)]) } p0=c(c0=-10,k=.04) fexp<-function(p,dat) { attach(dat) c0=p[1];k=p[2]; mn=exp(c0+k*age)*py ncases=cases detach(dat) -sum(ncases*log(mn) - mn) } windows(width=12,height=7) par(mfcol=c(2,4)) for (k in c(2,5,6,8)) { load(paste("c:/SEER02/",pre,filecanc[k],".RData",sep="")) # this loads in DF DF=DF[!is.na(DF$ICD9),] # get rid of missing ICD9 entries if (k==2){canc=c('stomach'); code=c(1510,1519); agLim=c(7,15)} if (k==5){canc=c("Brain"); code=c(1911,1919); agLim=c(5,15)} if (k==6){canc=c("Lung"); code=c(1622,1629); agLim=c(5,13)} if (k==8){canc=c("NHL"); code=c(2000,2008); agLim=c(5,16)} Indx=((DF$ICD9>=code[1])&(DF$ICD9<=code[2]))&(DF$numPrim==1)&((DF$ageRC>=agLim[1])&(DF$ageRC<=agLim[2])) ndat=DF[Indx,] attach(ndat) m=hist(ageRC[sex==1],breaks=seq(-.5,18.5,1),plot=FALSE) f=hist(ageRC[sex==2],breaks=seq(-.5,18.5,1),plot=FALSE) int=agLim[1]:agLim[2] age=c(0.5,3,seq(7.5,87.5,5))[int] ssolm=optim(p0,fexp,method="L-BFGS-B",control=list(maxit=400),dat=list(age=age,cases=m$counts[int],py=pym[int])) ssolf=optim(p0,fexp,method="L-BFGS-B",control=list(maxit=400),dat=list(age=age,cases=f$counts[int],py=pyf[int])) ym=exp(ssolm$par["c0"]+ssolm$par["k"]*age) yf=exp(ssolf$par["c0"]+ssolf$par["k"]*age) plot(age,m$counts[int]/pym[int],log="y",type='b',col='blue',ylab="Cases per Person-Years",main=paste(canc," Males ; k = ",format(c(ssolm$par["k"], 0.123456789), digits=2)[1])) lines(age,ym) plot(age,f$counts[int]/pyf[int],log="y",type='b',col='red' ,ylab="Cases per Person-Years",main=paste(canc," Females; k = ",format(c(ssolf$par["k"], 0.123456789), digits=2)[1])) lines(age,yf) print(ssolm$par) print(ssolf$par) detach(ndat) } # k loop on files detach(pops)