# Exam 1. 
#
# 1 (40pts). Compute the net lifetime low dose risk of CML in the absence of competing risks assuming 
#        mn = (exp(c1+k*age) + sv*(tsx^2)*exp(c2-kt*tsx) + sv*exp(c3-(tsx-mu)^2/sd2))*py
#        first for hiroshima males, then hiroshima females, then for both sexes, then both sexes and both cities. 
#        How did the Nagasaki data alter the risk estimates?
#
#  2 (30pts). Repeat the calculations using the same models but in the presence of (US 2001) competing risks. Assume an exposure at the age of 40 years. 
#        Which of the two risks estimates decreased by a greater percentage, HM or HF? Why?
#
#  3 (30pts). a) Since early diagnosis and early treatment improve treatment prognosis, one can conjecture that if one could survive the acute effects 
#        of a very high dose, the same dose which caused the cancer could cure it too. To explore this, remove the kerma constraint to obtain an
#        additional dose group at 5 Sv. Now let the dose-response speak for HM CML as in block 2 of CMLriskHM.r. What happens in the highest dose group?
#        Now do the same for AML and ALL with sexes and cities combined; assume a pure exponential waiting time function exp(-kt*tsx). Are the three 
#        leukemias self-consistent in the highest dose group? Describe the observed phenomenon.
#        b) The low dose regions of the leukemia dose-response curves differ slightly. Treating the dose-response parameters of 3a as "data", fit the 
#        seven data points to the model exp(c)*(sv+boa*sv^2)^k*exp(-ak*sv) where k is the smallest integer that leads to an acceptable fit. 
#        For chromosomal translocations and gamma rays, beta/alpha (of alphaD+betaD^2) is .06/(.01 to.02) or roughly 3-6. Approximating A-bomb 
#   	 doses in sieverts as gamma ray doses in gray, how many chromosomal tranlocation hits k are needed to cause CML?, AML?, and ALL? 
#  

# Please perform all of your work using R. Provide answers, supported by numbers and your codes.



# ************************** Solutions  *****************************************************

setwd("C:/cwru/classes/epbi473_05/lec4(riskModels)")
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:/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

# 1

par(mfrow=c(2,2))
run=c("HM","HF","HMF","HNMF")
risk=rep(0,4)
params=NULL
for (i in 1:4)
{
nd=d  # keep d as the full data frame
if (i!=4) nd=nd[nd$city==1,] # comment to get both cities
if (i<3) nd=nd[nd$sex==i,] # 1 for HM's, change to 2 for HF or comment for both sexes 
attach(nd)
py=10^4*adjPY
X0=c(c1=-13,k=0.05,rep(-10,10));X0 #am km L
fL<-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)
}
solL=optim(X0,fL,method="L-BFGS-B",control=list(maxit=400));solL  # the idea here is to let the fit speak through L, a bit like a one way anova
#
fixed=c(F,F,F,F,F,F,F)
X0=c(c1=-13,k=0.05,c2=-8,kt=0.4,c3=-9.4,mu=25,sd2=100)
fgauss<-function(x)
{
x[fixed]=X0[fixed]
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(X0,fgauss,method="L-BFGS-B",control=list(maxit=400));solg
#
t=seq(0,45,0.1)
attach(as.list(solg$par))
w=(t^2)*exp(c2-kt*t)+ exp(c3-(t-mu)^2/sd2)
detach(as.list(solg$par))
params=rbind(params,solg$par)
wait=exp(solL$par[3:12])
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)
cmlRiskD=sum(w)*.1       # see plots to see that integrals are OK w.r.t the function vanishing in time
cmlRiskD  # this is the lifetime excess risk of CML after a 1 Sv exposure.  
risk[i]=cmlRiskD
#mtext(paste(run[i],"Lifetime CML risk w/o competing risks =",format(cmlRiskD,digits=2),"per Sv"),line=-1)
mtext(paste(run[i]),line=-1)
paste(run[i],"lifetime CML risk w/o competing risks =",format(cmlRiskD,digits=2),"per Sv")
detach(nd)
} # for loop on conditions
par(mfrow=c(1,1))
rownames(params)<-run
params
risk

#
# For HM, second component is essentially zero, with c3 as it is
#          c1            k           c2           kt           c3           mu          sd2 
#-12.83788493   0.04137065  -8.44643679   0.39821282 -13.12321370  24.88421485  99.99366265 
# "Lifetime CML risk w/o competing risks = 0.0068 per Sv"

#
# For HF, second hump is strong, first is present but weak
#          c1            k           c2           kt           c3           mu          sd2 
#-15.04721267   0.06404139  -7.21642398   0.74604806  -8.89047517  25.20515380  99.89571947 
# Lifetime CML risk w/o competing risks = 0.006 per Sv"

#
# HF + HM combined (not asked for)
#          c1            k           c2           kt           c3           mu          sd2 
#-13.97109664   0.05269812  -8.09333465   0.51871975  -9.19365849  24.93846147  99.99505584 
#"Lifetime CML risk w/o competing risks = 0.0062 per Sv"

#
# both cities and sexes
#          c1            k           c2           kt           c3           mu          sd2 
#-13.96272303   0.04950036  -8.47953603   0.48822456  -9.48395618  24.93359908  99.99452218 
#"Lifetime CML risk w/o competing risks = 0.0049 per Sv"



# 2  **************************************************

ageX=40  # do everything for an exposed 40 year old
agef=seq(ageX,100,1)

# first find survival probability at age a>x given exposed (and alive) at age x 
lt<-read.table("lifeUSA2001.txt",header=T)
m=lt[lt$Sex=="males",];m
mage=c(0.5,3,seq(7.5,97.5,5),105)
m=m[mage>40,];m
mu=m$nqx # mortality rate
mu=rep(mu,each=5)
mu=mu[1:61]/5
Mu=cumsum(mu)
S=exp(-Mu);S
plot(agef,S,type="l",xlab="Age",ylab="Survival Probability",main="US Population in 2001")
#locator(n=1)

Rccr=NULL
for (i in 1:4)
{
attach(as.list(params[i,]))
Rccr[i]=sum( ((agef-ageX)^2*exp(c2-kt*(agef-ageX))+exp(c3-(agef-ageX-mu)^2/sd2)) *S  )  # the sum here is across the fine age grid
detach(as.list(params[i,]))
}

Rccr   # CML risk, controlled for competing other risks
# 0.006589233 0.004735048 0.005202627 0.004172087
# male risk only slightly reduced, female risk is more reduced
# The reason is that the female waiting times are much larger, so there is more time for competing CODs to take their lives



# 3  **************************************************

#            HM CML
#
# read in again to try without kerma constraint
setwd("C:/cwru/classes/epbi473_05/lec4(riskModels)")
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:/abomb/hema87.dat", header=F,col.names=cols);
d=d[d$adjPY>0,] #remove two recs with zero py
#d=d[d$kerma==1,] # commented to take all kermas and thus higher doses
d=d[d$city==1,] # hiro
d=d[d$sex==1,] # 1 for males
attach(d)
py=10^4*adjPY
svc=cut(sv,c(-1,.02,.2,.5,1,2,4,10))
table(svc[CML>0])   # a decent spread of cases
sum(table(svc[CML>0])) 
mdose=c(0,.1,.35,.75,1.5,3,5)
D0=c(-40,rep(-5,6));names(D0)<-paste("Ds",1:7,sep="");D0
X0=c(c1=-13.6,k=0.047,kt=0.5,D0)
fixed=c(F,F,F,T,rep(F,6))
fD<-function(x){
x[fixed]=X0[fixed]
 c1=x[1];k=x[2];kt=x[3];Ds=x[4:10]  		
 mn =( exp(c1+k*age) + exp(Ds[svc])*tsx^2*exp(-kt*tsx))*py;
 -sum(CML*log(mn) - mn)
 }
solDCML=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))
solDCML
ampsD=exp(solDCML$par[4:10])
plot(mdose,ampsD,pch=5,xlab="Sv",ylab="amplitude",main="Excess CML risk dose-response shape")
X0=c(c2=-8,boa=3,alfk=-.1)
fixed=c(F,F,F)
fD<-function(x){
x[fixed]=X0[fixed]
 c2=x[1];boa=x[2];alfk=x[3]  		
 E =exp(c2)*(mdose+boa*mdose^2)*exp(-alfk*mdose)
 sum((ampsD-E)^2)
 }
sol=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))
sol
fD=seq(0,5,.1)
attach(as.list(sol$par))
y=exp(c2)*(fD+boa*fD^2)*exp(-alfk*fD)
lines(fD,y)
detach(as.list(sol$par))




#                    AML  
#
setwd("C:/cwru/classes/epbi473_05/lec4(riskModels)")
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:/abomb/hema87.dat", header=F,col.names=cols);
d=d[d$adjPY>0,] #remove two recs with zero py
attach(d)

py=10^4*adjPY
svc=cut(sv,c(-1,.02,.2,.5,1,2,4,10))

table(svc[AML>0])   # a decent spread of cases
sum(table(svc[AML>0])) 
fixed=c(F,F,F,T,rep(F,6))
D0=c(-40,rep(-10,6));names(D0)<-paste("Ds",1:7,sep="");D0
X0=c(c1=-12.6,k=0.042,kt=0.411,D0)
fD<-function(x){
 x[fixed]=X0[fixed]
 c1=x[1];k=x[2];kt=x[3];Ds=x[4:10]   		
 mn =( exp(c1+k*age) + exp(Ds[svc])*exp(-kt*tsx) )*py;
 -sum(AML*log(mn) - mn)
 }
solDAML=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))

solDAML
ampsD=exp(solDAML$par[4:10])
plot(mdose,ampsD,pch=5,xlab="Sv",ylab="amplitude",main="Excess AML risk dose-response shape")
X0=c(c2=-8,boa=3,alfk=-.1)
fixed=c(F,F,F)
fD<-function(x){
x[fixed]=X0[fixed]
 c2=x[1];boa=x[2];alfk=x[3]  		
 E =exp(c2)*(mdose+boa*mdose^2)*exp(-alfk*mdose)
 sum((ampsD-E)^2)
 }
sol=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))
sol
fD=seq(0,5,.1)
attach(as.list(sol$par))
y=exp(c2)*(fD+boa*fD^2)*exp(-alfk*fD)
lines(fD,y)
detach(as.list(sol$par))



#                    ALL  
#
setwd("C:/cwru/classes/epbi473_05/lec4(riskModels)")
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:/abomb/hema87.dat", header=F,col.names=cols);
d=d[d$adjPY>0,] #remove two recs with zero py
attach(d)

py=10^4*adjPY
svc=cut(sv,c(-1,.02,.2,.5,1,2,4,10))

table(svc[ALL>0])   # a decent spread of cases
sum(table(svc[ALL>0])) 
fixed=c(F,F,F,T,rep(F,6))
D0=c(-40,rep(-10,6));names(D0)<-paste("Ds",1:7,sep="");D0
X0=c(c1=-12.6,k=0.042,kt=0.411,D0)
fD<-function(x){
 x[fixed]=X0[fixed]
 c1=x[1];k=x[2];kt=x[3];Ds=x[4:10]   		
 mn =( exp(c1+k*age) + exp(Ds[svc])*exp(-kt*tsx) )*py;
 -sum(ALL*log(mn) - mn)
 }
solDALL=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))
solDALL
ampsD=exp(solDALL$par[4:10])
plot(mdose,ampsD,pch=5,xlab="Sv",ylab="amplitude",main="Excess ALL risk dose-response shape")
X0=c(c2=-8,boa=3,alfk=1)
fixed=c(F,F,F)
fD<-function(x){
x[fixed]=X0[fixed]
 c2=x[1];boa=x[2];alfk=x[3]  		
 E =exp(c2)*(mdose+boa*mdose^2)*exp(-alfk*mdose)
 sum((ampsD-E)^2)
 }
sol=optim(X0,fD,method="L-BFGS-B",control=list(maxit=400))
sol

fD=seq(0,5,.1)
attach(as.list(sol$par))
y=exp(c2)*(fD+boa*fD^2)*exp(-alfk*fD)
lines(fD,y)
detach(as.list(sol$par))