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$sv<1,] # #d=d[(d$sv<.2)|(d$sv>=1),] # attach(d) summary(d) py=10^4*adjPY X0=c(c1=-13,k=0.05,rep(-10,10));X0 #am km L flin<-function(x) { c1=x[1];k=x[2];L=x[3:12]; mn = exp(c1+k*age)*py + sv*exp(L[calg])*py; -sum(CML*log(mn) - mn) } sol=optim(X0,flin,method="L-BFGS-B",control=list(maxit=400));sol # the idea here is to let the data speak through L, a bit like a one way anova wait=exp(sol$par[3:12]) X1=c(c1=-13,k=0.05,c2=-8,kt=0.4,c3=-9.4,mu=25,sd2=100) flingam2<-function(x) { c1=x[1];k=x[2];c2=x[3];kt=x[4];c3=x[5];mu=x[6];sd2=x[7]; mn = (exp(c1+k*age) + sv*(tsx^2)*exp(c2-kt*tsx)+ sv*exp(c3-(tsx-mu)^2/sd2))*py; -sum(CML*log(mn) - mn) } solg=optim(X1,flingam2,method="L-BFGS-B",control=list(maxit=400));solg t=seq(0,45,0.1) attach(as.list(solg$par)) #w=(t^2)*exp(solg$par[3]-solg$par[4]*t)+exp(solg$par[5]-(t-solg$par[6])^2/solg$par[7]) w=(t^2)*exp(c2-kt*t)+ exp(c3-(t-mu)^2/sd2) year1=t+1945 year2=c(1951,1953,1956,1959,1963,1968,1973,1978,1983,1987) plot(year1,w,type="l",xlab="Year",ylab="w(t)",main="IR-to-CML waiting time density shape") points(year2,wait,pch=5) #text(1970,2e-4,paste("kt=",format(c(kt, 0.123456789), digits=2)[1],"; mu=",format(c(mu, 0.123456789), digits=2)[1]))