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","MM") d<-read.table("c:/data/abomb/hema87.dat", header=F,col.names=cols); d=d[d$adjPY>0,] #remove two recs with zero py d=d[d$kerma==1,] # take only kerma < 4 Gy d=d[d$city==1,] # hiroshima only d$calg=as.integer(cut(d$calg,c(0,2,4,6,8,10))) head(d$calg,20) d$agex=as.integer(cut(d$agex,c(0,20,40,60,80))) agem=c(10,30,50,70) N=5 M=3 X0=c(c1=-13,k=0.05,rep(-10,N),rep(-5,M));X0 #am km L years=c(1951,1953,1956,1959,1963,1968,1973,1978,1983,1987) years=c(1952,1958,1966,1976,1985) par(mfrow=c(1,1),mar=c(4.7,4.1,1.3,1)) m=d[d$sex==1,] # males only f=d[d$sex==2,] # females only flin<-function(x) { c1=x[1];k=x[2];L=x[3:(2+N)];AX=c(0,x[(2+N+1):(2+N+3)]); mn = exp(c1+k*age)*py + sv*exp(L[calg]+AX[agex])*py; -sum(CML*log(mn) - mn) } getCI<-function(sol) { if (det(sol$hessian)>0) { sig=sqrt(diag(solve(sol$hessian))) print("Hessian OK") } else { sig=Inf print("Hessian Singular") } sig upper=signif(sol$par+1.96*sig,5) lower=signif(sol$par-1.96*sig,5) point=signif(sol$par,5) CI=cbind(lower,point,upper) CI=exp(CI) CI } attach(m) py=10^4*adjPY sol=optim(X0,flin,method="L-BFGS-B",hessian=TRUE,control=list(maxit=400));sol # the idea here is to let the data speak through L, a bit like a one way anova detach(m) waita=exp(sol$par[3:(2+N)]) Xa=exp(sol$par[(2+N+1):(2+N+3)]) getCI(sol) attach(f) py=10^4*adjPY sol=optim(X0,flin,method="L-BFGS-B",hessian=TRUE,control=list(maxit=400));sol # the idea here is to let the data speak through L, a bit like a one way anova detach(f) waitb=exp(sol$par[3:(2+N)]) Xb=exp(sol$par[(2+N+1):(2+N+3)]) getCI(sol) plot(years,waita,pch=1,xlab="Year",ylab="CML/(Sievert * Person-Year)", ylim=c(0,1e-3),col="blue") points(years,waitb,pch=2,col="red") legend(1967,1e-3,c("Males","Females"),col=c("blue","red"),pch=1:2) #mtext("A",side=3,line=.1,cex=1.5,adj=0)