#
# R CODE FOR REPRODUCING CONTOUR FIGURE 6.10 IN JW'S 
#       "BAYESIAN AND FREQUENTIST REGRESSION ANALYSIS" CHAPTER 6.
# CODE WRITTEN BY RUTH SALWAY
##############################################
#
Plot.PKlik.comp<-function(data,range,fixed,n=10, log=F, lev=NULL,
    nlev=10,xlab=NULL,ylab=NULL,llmax=0)
{
# data is a data.frame with columns conc, dose and time
# range is  the range for each of the two parameters on the x and y axes
# fixed lists the values at which the other params are kept fixed
# both are lists with elements names V, ka and ke as required
#
# log=T plots on the log(V) scale for normal lik - this is the lik 
# actually maximised. Check all parameters are given
#
data$logconc=log(data$conc)
    pnames=c("V","ka","ke","sigma")
    fnames=names(fixed)
    rnames=names(range)
    if (any((c(fnames,rnames) %in% pnames)==F)) 
                stop(paste("fixed and range must contain ",
                paste(pnames,collapse=", ",sep="")))
    if (any(fnames %in% rnames)==T) 
       stop("Each parameter can appear in either fixed or range, not both.")
    if (length(fnames)!=2) stop("fixed must be a list with two parameters.")
    if (length(rnames)!=2) stop("range must be a list with two parameters.")
#
# Set up parameters
#
    if ("V" %in% fnames) {
        V=rep(fixed$V,n*n)
        if (log) V=log(V)
    }
    else {
        if (log) V=seq(log(range$V[1]), log(range$V[2]), len = n)
        else V=seq(range$V[1], range$V[2], len = n)
    }
    if ("ka" %in% fnames) {
        ka=rep(fixed$ka,n*n)
        if (log) ka=log(ka)
    }
    else {
        if (log) ka=seq(log(range$ka[1]), log(range$ka[2]), len = n)
        else ka=seq(range$ka[1], range$ka[2], len = n)
    }
    if ("ke" %in% fnames) {
        ke=rep(fixed$ke,n*n)
        if (log) ke=log(ke)
    }
    else {
        if (log) ke=seq(log(range$ke[1]), log(range$ke[2]), len = n)
        else ke=seq(range$ke[1], range$ke[2], len = n)
    }
    if ("sigma" %in% fnames) sigma=rep(fixed$sigma,n*n)
    else sigma=seq(range$sigma[1], range$sigma[2], len = n)
#
    eval(parse(text=paste(rnames[1],"=rep(",rnames[1],",n)",sep="")))
    eval(parse(text=paste(rnames[2],"=rep(",rnames[2],",rep(n,n))",sep="")))
#
# parameters should now be n*n vectors
#
    ll=rep(0,n*n)
    for (i in 1:(n*n)) {
        if (!log) {
            fit=log( data$dose*ka[i]/(V[i]*(ka[i]-ke[i])) * 
                  (exp(-ke[i]*data$time)-exp(-ka[i]*data$time)) )
            ll[i] = loglik.lnorm(data$conc,fit,sigma[i])
        }
        else {
            fit = log(data$dose) + ka[i] - V[i] - log(exp(ka[i])-exp(ke[i])) +
                    log( exp(-exp(ke[i])*data$time)-exp(-exp(ka[i])*data$time))
            ll[i] = loglik.norm(data$logconc,fit,sigma[i])
#if(ll[i] >llmax)cat("llmax lli V ka ke: ",llmax,ll[i],V[i],ka[i],ke[i],"\n")
        }
    }
ll <- ll - llmax
cat("max ll = ",max(ll),"\n")
    if (is.null(lev))  lev=pretty(range(ll, finite = TRUE),nlev)

    if (!log) {
if (is.null(xlab)) xlab=rnames[1]
if (is.null(ylab)) ylab=rnames[2]
eval(parse(text=paste("contour(seq(range$",rnames[1],"[1],range$",
      rnames[1],"[2],len=n) ,seq(range$",rnames[2],"[1],range$",
      rnames[2],"[2],len=n), matrix(ll,n,n), lev = lev,xlab='",xlab,"',
      ylab='",ylab,"')",sep="")))
#        title("Compartment")
    }
    else {
if (is.null(xlab)) xlab=paste("log(",rnames[1],")",sep="")
if (is.null(ylab)) ylab=paste("log(",rnames[2],")",sep="")
eval(parse(text=paste("contour(seq(log(range$",rnames[1],"[1]),
     log(range$",rnames[1],"[2]),len=n) ,seq(log(range$",rnames[2],"[1]),
     log(range$",rnames[2],"[2]),len=n), matrix(ll,n,n), lev = lev,xlab='",
     xlab,"',ylab='",ylab,"')",sep="")))
#        title("Compartment")
    }
}
#+++++++++++++++++++++++++++++++++++++
loglik.norm <- function(y,mu,sigma){
    n <- length(y)
    -n*log(sigma) - 0.5*n*log(2*pi) - 0.5 * sum((y-mu)^2)/(sigma^2)
}
#+++++++++++++++++++++++++++++++++++++
loglik.lnorm <- function(y,mu,sigma){
         n <- length(y)
         -n*log(sigma) - sum(log(y)) - 0.5*n*log(2*pi) - 
         0.5 * sum((log(y)-mu)^2)/(sigma^2)
    }
#+++++++++++++++++++++++++++++++++++++
loglik.gam <- function(y,alpha,beta){
        n <- length(y)
        # Y ~ Ga(a,b) with mean=a/b var=a/b^2
        -n*lgamma(alpha) + sum(alpha*log(beta) +(alpha-1)*log(y)-beta*y)
    } 
##############################################
# Code to produce plots
##############################################
#
# Plot likelihood contours for compartment model
#
y <- Theoph$conc[24:33];x <- Theoph$Time[24:33];dose <- Theoph$Dose[24] #Ind 3
z <- log(y)
pkdata <- data.frame(z=z,t=x,dose=dose)
modpk <- nls(z~log(dose)+theta1-theta0-log(exp(theta1)-exp(theta2))+
         log( exp(-exp(theta2)*t)-exp(-exp(theta1)*t) ),data=pkdata,
         start=list(theta0=log(1.14),theta1=log(0.565),theta2=log(0.109)))
#
# range over which to plot: nsd gives number of sds in each direction
#
nsd=3 # Change to 30 for wider range in book (Fig 6.10)
Vr=exp(coef(modpk)[1]+c(-1,1)*nsd*sqrt(vcov(modpk)[1,1]))
kar=exp(coef(modpk)[2]+c(-1,1)*nsd*sqrt(vcov(modpk)[2,2]))
ker=exp(coef(modpk)[3]+c(-1,1)*nsd*sqrt(vcov(modpk)[3,3]))
conc <- y
time <- x
dose <- dose
pkdata2 <- data.frame(conc,time,dose)
#These are the MLEs - your code will presumably vary
V.est=exp(coef(modpk)[1])
ka.est=exp(coef(modpk)[2])
ke.est=exp(coef(modpk)[3])
sig.est=0.06322
# fit at the max
fitmax = log(pkdata2$dose) + log(ka.est) - log(V.est) - log(ka.est-ke.est) +
         log( exp(-ke.est*pkdata2$time)-exp(-ka.est*pkdata2$time))
llmax = loglik.norm(log(pkdata2$conc),fitmax,sig.est)
#
# fix ke and sigma at estimated values, and plot V and ka on axes
#
levels30 <- c(-10000,-8000,-6000,-4000,-2000,-1000,-500,-250,-100,-50,-20)
pdf("logliklogVlogka3.pdf",height=4,width=4)
par(mfrow=c(1,1),pty="s")
#
# uncomment out the following if 30 sds wanted 
#
#Plot.PKlik.comp(data=pkdata2,range=list(V=Vr,ka=kar),fixed=list(ke=ke.est,sigma=sig.est),lev=levels30,n=50,log=T,xlab="",ylab="",llmax=llmax)
Plot.PKlik.comp(data=pkdata2,range=list(V=Vr,ka=kar),fixed=list(ke=ke.est,sigma=sig.est),n=50,log=T,xlab="",ylab="",llmax=llmax)
mtext(expression(paste("log ",V,,sep="")),side=1,line=3)
mtext(expression(paste("log ",k[a],sep="")),side=2,line=3)
points(coef(modpk)[1],coef(modpk)[2],pch=19)
dev.off()
pdf("logliklogVlogke3.pdf",height=4,width=4)
par(mfrow=c(1,1),pty="s")
#Plot.PKlik.comp(data=pkdata2,range=list(V=Vr,ke=ker),fixed=list(ka=ka.est,sigma=sig.est),lev=levels30,n=50,log=T,xlab="",ylab="",llmax=llmax)
Plot.PKlik.comp(data=pkdata2,range=list(V=Vr,ke=ker),fixed=list(ka=ka.est,sigma=sig.est),n=50,log=T,xlab="",ylab="",llmax=llmax)
mtext(expression(paste("log ",V,sep="")),side=1,line=3)
mtext(expression(paste("log ",k[a],sep="")),side=2,line=3)
points(coef(modpk)[1],coef(modpk)[3],pch=19)
dev.off()
pdf("logliklogkalogke3.pdf",height=4,width=4)
par(mfrow=c(1,1),pty="s")
#Plot.PKlik.comp(data=pkdata2,range=list(ka=kar,ke=ker),fixed=list(V=V.est,sigma=sig.est),lev=levels30,n=50,log=T,xlab="",ylab="",llmax=llmax)
Plot.PKlik.comp(data=pkdata2,range=list(ka=kar,ke=ker),fixed=list(V=V.est,sigma=sig.est),n=50,log=T,xlab="",ylab="",llmax=llmax)
mtext(expression(paste("log ",k[a],sep="")),side=1,line=3)
mtext(expression(paste("log ",k[e],sep="")),side=2,line=3)
points(coef(modpk)[2],coef(modpk)[3],pch=19)
dev.off()
#abline(h=coef(modpk)[3],lwd=2,lty=3)
stop()
#
# Plot likelihood contours for GLM model: these are not included in the
# book chapter
#
#range over which to plot: nsd gives number of sds in each direction
b0r=glm$parameters[1]+c(-1,1)*nsd*sqrt(glm$var[1,1])
b1r=glm$parameters[2]+c(-1,1)*nsd*sqrt(glm$var[2,2])
b1r[b1r>0]=0
b2r=glm$parameters[3]+c(-1,1)*nsd*sqrt(glm$var[3,3])
b2r[b2r>0]=0
#
# MLEs from fitted model
#
b2.est=glm$parameters[3]
phi.est=glm$phi
#
# fix beta2 and plot beta0 and beta1 on axes
#
Plot.PKlik.glm(data=pkdata,range=list(beta0=b0r,beta1=b1r),fixed=list(beta2=b2.est,phi=phi.est),lev=c(-10000,-8000,-6000,-4000,-2000,-1000,-500,-250,-100,-50),n=50,xlab="",ylab="")
mtext(expression(beta[0]),side=1,line=3)
mtext(expression(beta[1]),side=2,line=3)
points(glm$parameters[1],glm$parameters[2],pch=19)
#
# plot on beta0 scale
#
Plot.PKlik.glm<-function(data,range,fixed, n=10, lev=NULL,nlev=10,xlab=NULL,ylab=NULL)
{
#
# data is a data.frame with columns conc, dose and time
# range is  the range for each of the two parameters on the x and y axes
# fixed lists the values at which the other params are kept fixed
# both are lists with elements names beta0, beta1 and beta2 as required
#
data$invtime=1/data$time
#
# Check all parameters are given
#
    pnames=c("beta0","beta1","beta2","phi")
    fnames=names(fixed)
    rnames=names(range)
    if (any((c(fnames,rnames) %in% pnames)==F)) stop(paste("fixed and range must contain ",paste(pnames,collapse=", ",sep="")))
    if (any(fnames %in% rnames)==T) stop("Each parameter can appear in either fixed or range, not both.")
    if (length(fnames)!=2) stop("fixed must be a list with two parameters.")
    if (length(rnames)!=2) stop("range must be a list with two parameters.")
#
# Set up parameters
#
    if ("beta0" %in% fnames) beta0=rep(fixed$beta0,n*n)
        else beta0=seq(range$beta0[1], range$beta0[2], len = n)
    if ("beta1" %in% fnames) beta1=rep(fixed$beta1,n*n)
        else beta1=seq(range$beta1[1], range$beta1[2], len = n)
    if ("beta2" %in% fnames) beta2=rep(fixed$beta2,n*n)
        else beta2=seq(range$beta2[1], range$beta2[2], len = n)
    if ("phi" %in% fnames) phi=rep(fixed$phi,n*n)
        else phi=seq(range$phi[1], range$phi[2], len = n)
    eval(parse(text=paste(rnames[1],"=rep(",rnames[1],",n)",sep="")))
    eval(parse(text=paste(rnames[2],"=rep(",rnames[2],",rep(n,n))",sep="")))
#
# parameters should now be n*n vectors
#
    ll=rep(0,n*n)
    for (i in 1:(n*n)) {
        fit=data$dose*exp(beta0[i]+beta1[i]*data$time+beta2[i]*data$invtime)
        ll[i] <- loglik.gam(data$conc,1/phi[i],1/(fit*phi[i]))
    }
#    
    if (is.null(lev))  lev=pretty(range(ll, finite = TRUE),nlev) 
    if (is.null(xlab)) xlab=rnames[1]
    if (is.null(ylab)) ylab=rnames[2]
    eval(parse(text=paste("contour(seq(range$",rnames[1],"[1],range$",
         rnames[1],"[2],len=n) ,seq(range$",rnames[2],"[1],range$",
         rnames[2],"[2],len=n), matrix(ll,n,n), lev = lev,xlab='",
         xlab,"',ylab='",ylab,"')",sep="")))
#
    title("GLM")
}