R/Dose.Response.R

Defines functions DR.Lyman DR.Goitein DR.Niemierko DR.Munro DR.Okunieff DR.Warkentin DR.Bentzen outcome.def NTCP.curve probit.nLL.Lyman logit.nLL.Niemierko fit.NTCP fit.NTCP.CI print.NTCP summary.NTCP outcome.NTCP plot.NTCP plot3D.NTCP probit.nLL.Lyman.TD50 logit.nLL.Niemierko.TD50 probit.nLL.Lyman.gamma50 logit.nLL.Niemierko.gamma50 probit.nLL.Lyman.a logit.nLL.Niemierko.a CI bisect proflik calc.boot.el bootstrap.CI NTCP.kernel fit.DR ProfLik.DR

Documented in DR.Bentzen DR.Goitein DR.Lyman DR.Munro DR.Niemierko DR.Okunieff DR.Warkentin

###########################################################################
## script generating dose/response models over two  over two parameters  ##
## (TD50, gamma50), or three parameters (TD50, gamma50, a)               ##
###########################################################################

## Dose/Response according Lyman 1985, Kutcher and Burman 1989, Deasy 2000 (probit) ##
DR.Lyman <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series 
  if (!is.null(dose)) p <- pnorm(q=((dose - TD50)*gamma50*sqrt(2*pi))/TD50)
  if (!is.null(diffdvh)) p <- pnorm(q=((EUD(dvh.matrix=diffdvh, a=aa) - TD50)*gamma50*sqrt(2*pi))/TD50)
  return(p)
}

## Dose/Response according Goitein 1979 (logprobit) ##
DR.Goitein <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series 
  if (!is.null(dose)) p <- pnorm(q=(log(dose/TD50)*gamma50*sqrt(2*pi)))
  if (!is.null(diffdvh)) p <- pnorm(q=(log(EUD(dvh.matrix=diffdvh, a=aa)/TD50)*gamma50*sqrt(2*pi)))
  return(p)
}

## Dose/Response according Niemierko 1991 (loglogit) ##
DR.Niemierko <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series  
  if (is.vector(dose) && is.null(diffdvh)) p <- 1/(1+(TD50/dose)^(4*gamma50))
  if (is.null(dose) && is.matrix(diffdvh)) p <- 1/(1+(TD50/EUD(dvh.matrix=diffdvh, a=aa))^(4*gamma50))
  return(p)
}

## Dose/Response according Munro, Gilbert, Kallman 1992 (Poisson approximation) ##
DR.Munro <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series  
  if (is.vector(dose) && is.null(diffdvh)) p <- 2^(-(exp(exp(1)*gamma50*(1-dose/TD50))))
  if (is.null(dose) && is.matrix(diffdvh)) p <- 2^(-(exp(exp(1)*gamma50*(1-EUD(dvh.matrix=diffdvh, a=aa)/TD50))))
  return(p) 
}

## Dose/Response according Okunieff 1995 (logit) ##
DR.Okunieff <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series  
  if (is.vector(dose) && is.null(diffdvh)) p <- 1/(1+exp(4*gamma50*(1-(dose/TD50))))
  if (is.null(dose) && is.matrix(diffdvh)) p <- 1/(1+exp(4*gamma50*(1-(EUD(dvh.matrix=diffdvh, a=aa)/TD50))))
  return(p) 
}

## Dose/Response according Warkentin 2004 (Poisson) ##
DR.Warkentin <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series  
  if (is.vector(dose) && is.null(diffdvh)) p <- 0.5^(exp(2*gamma50/log(2)*(1-dose/TD50)))
  if (is.null(dose) && is.matrix(diffdvh)) p <- 0.5^(exp(2*gamma50/log(2)*(1-EUD(dvh.matrix=diffdvh, a=aa)/TD50)))
  return(p) 
}

## Dose/Response according Bentzen 1997 (log Poisson) ##
DR.Bentzen <- function (TD50=45, gamma50=1.5, aa=1, diffdvh=NULL, dose=NULL) {
  # check the single choice between dvh matrix or dose series
  if (!is.null(dose) && !is.null(diffdvh)) stop("Select either a DVH or a point dose to calculate NTCP")
  # check the single choice between dvh matrix or dose series  
  if (is.vector(dose) && is.null(diffdvh)) p <- 0.5^(TD50/dose)^(2*gamma50/log(2))
  if (is.null(dose) && is.matrix(diffdvh)) p <- 0.5^(TD50/EUD(dvh.matrix=diffdvh, a=aa))^(2*gamma50/log(2))
  return(p) 
}

## definition of outcome based by dose
outcome.def <- function (dvh.matrix=NULL, dose=NULL, TD50=45, gamma50=2, a=1, 
                         model=c("Lyman", "Niemierko", "Okunieff", "Munro")) {
  # checks for not empty both dose and dvh.matrix 
  if (is.matrix(dvh.matrix) && is.vector(dose)) 
    stop("You can choice to simulate the outcome EITHER for a dvh based patients series ", "\n", 
         "OR for a single doses vector patients series")
  # check matrix input in dvh.matrix
  if (!is.matrix(dvh.matrix) && is.null(dose))    
    stop("dvh.matrix MUST be a matrix containing relative differential DVHs")
  # check vector input in dose
  if (!is.vector(dose) && is.null(dvh.matrix))
    stop("dose MUST be a vector containing patients doses")
  # generating outcome simulation for a dvh matrix
  if (is.matrix(dvh.matrix)) { 
    if ((dvh.matrix[1,2]==1) && (dvh.matrix[2,2]>0)) 
      warning("dvh.matrix seems to be a cumulative DVH, ", "\n", 
              "  If so, and subsequent fit.NTCP converges,", "\n",
              "  results are not consisntent.", "\n",
              "  If dvh.matrix is a cumulative DVH matrix", "\n",
              "  convert into a differential matrix ", "\n" ,
              "  by using diff.dvh function.")
    dvh.size <- dim(dvh.matrix)                                       # size of NTCP series
    temp.outcome <- seq(1, dvh.size[2] - 1,1)                         # sequence of temporary outcomes
    eqdoses <- EUD(dvh.matrix=dvh.matrix, a=a)                        # sequence of EUD
    model<-match.arg(model)
    if (model=="Lyman") {
      NTCP.Series <- DR.Lyman(TD50=TD50, gamma50=gamma50, aa=a, diffdvh=dvh.matrix)
      for (m in 1:dvh.size[2]) {
        ifelse(runif(1,0,1) < NTCP.Series[m],temp.outcome[m] <- 1, temp.outcome[m] <-0) # outcome definition for Lyman model 
      }
    }
    if (model=="Niemierko") {
      NTCP.Series <- DR.Niemierko(TD50=TD50, gamma50=gamma50, aa=a, diffdvh=dvh.matrix)
      for (m in 1:dvh.size[2]) {
        ifelse(runif(1,0,1) < NTCP.Series[m],temp.outcome[m] <- 1, temp.outcome[m] <-0) # outcome definition for Niemierko 
      }
    }
    outcome <- matrix(nrow=length(temp.outcome), ncol=3)              # create the outcome matrix 
    outcome[,1] <- eqdoses                                            # 1st column: EUDs
    outcome[,2] <- NTCP.Series                                        # 2nd column: NTCP
    outcome[,3] <- temp.outcome                                       # 3rd column: Outcome
    colnames(outcome) <- c("Dose", "NTCP", "Outcome")
    return(outcome)
  }
  if (is.vector(dose)) {
    l <- length(x=dose)
    temp.outcome <- seq(from=1, to=l, 1)
    model<-match.arg(model)
    if (model=="Lyman") {
      NTCP.Series <- DR.Lyman(TD50=TD50, gamma50=gamma50, dose=dose)
      for (m in 1:l) {
        ifelse(runif(1,0,1) < NTCP.Series[m],temp.outcome[m] <- 1, temp.outcome[m] <-0) # outcome definition for Lyman model 
      }
    }
    if (model=="Niemierko") {
      NTCP.Series <- DR.Niemierko(TD50=TD50, gamma50=gamma50, dose=dose)
      for (m in 1:dvh.size[2]) {
        ifelse(runif(1,0,1) < NTCP.Series[m],temp.outcome[m] <- 1, temp.outcome[m] <-0) # outcome definition for Niemierko 
      }
    }
    outcome <- matrix(nrow=length(temp.outcome), ncol=3)              # create the outcome matrix 
    outcome[,1] <- dose                                               # 1st column: vector dose
    outcome[,2] <- NTCP.Series                                        # 2nd column: NTCP
    outcome[,3] <- temp.outcome                                       # 3rd column: Outcome
    colnames(outcome) <- c("Dose", "NTCP", "Outcome")
    return(outcome)
  }
}

## create NTCP curve
NTCP.curve <- function(TD50=45,gamma50=1.5,model=c("Lyman","Niemierko"),dmin=0,dmax=70,add=TRUE) {
  if (model=="Lyman") {
    curve(expr=pnorm(q=((x - TD50)*gamma50*sqrt(2*pi)/TD50)),from=dmin,to=dmax,ylim=c(0,1),
          add=add,xlab="Dose [Gy]",ylab="NTCP [*100%]")
  }
  if (model=="Niemierko") {
    curve(expr=1/(1+(TD50/EUD(dvh.matrix=diffdvh, a=aa))^(4*gamma50)),from=dmin,to=dmax,ylim=c(0,1),
          add=add,xlab="Dose [Gy]",ylab="NTCP [*100%]")
  }  
}



# Log-Likelihood of Lyman function
probit.nLL.Lyman <- function (z,dvh.matrix, outcome) {
  TD50 <- z[1]
  gamma50 <- z[2]
  aa <- z[3]
  p <- pnorm(gamma50*sqrt(2*pi)*(EUD(dvh.matrix, a=aa)/TD50-1))
  return(-sum(outcome*log(p)+(1-outcome)*log(1-p)))
}



# Log-Likelihood of Niemierko function
logit.nLL.Niemierko <- function (z,dvh.matrix, outcome) {
  TD50 <- z[1]
  gamma50 <- z[2]
  aa <- z[3]
  p <- 1/(1+(TD50/EUD(dvh.matrix, a=aa))^(4*gamma50))
  return(-sum(outcome*log(p)+(1-outcome)*log(1-p)))
}

# function for fitting DVH data
fit.NTCP <- function(model=c("Lyman","Niemierko"),dvh.matrix, outcome, 
                     calcCI=TRUE, verbose=FALSE, C.I.width=.95, n.boot=1000,
                     min.TD50=0, max.TD50=150, min.gamma50=0, max.gamma50=15,
                     min.a=-100, max.a=100, C.I.type=c("ProfLik", "Boot")) {
  model<-match.arg(model)
  C.I.type<-match.arg(C.I.type)
  ca1 <<- qchisq(C.I.width,1)/2  # set the bounds for C.I. calculation
  if (verbose==TRUE) message("Fitting ", model, " model")
  if (!is.vector(outcome)) stop("'outcome' MUST be a vector")
  if (verbose==TRUE) message("Outcome series contains ", length(outcome), " observations", "\n")
  if (model=="Lyman") {
    out <- nlminb(c(40,2,2), probit.nLL.Lyman,lower=c(10,0.01,-100),upper=c(180,15,100),
                  dvh.matrix=dvh.matrix,outcome=outcome)
  }
  if (model=="Niemierko") {
    out <- nlminb(c(40,2,2), logit.nLL.Niemierko,lower=c(10,0.01,-100),upper=c(180,15,100),
                  dvh.matrix=dvh.matrix,outcome=outcome)
  }
  # MLE calculation
  MLE <<- -out$objective
  # modelAIC<- 6 -2*out$objective # Akaike information criterion
  modelAIC <- -2*out$objective + 6 +(24 / (ncol(dvh.matrix)-5)) # Akaike information criterion with num cases correction
  # fitted parameters
  D50 <<- out$par[1]
  g50 <<- out$par[2]
  aaa <<- out$par[3]
  iter <- out$iterations
  if (verbose==TRUE) {
    message(model, " model successfully converged")
    message("      TD50 = ", format(D50, digits=3, nsmall=3, justify="right", width=8))
    message("   gamma50 = ", format(g50, digits=3, nsmall=3, justify="right", width=8))
    message("         a = ", format(aaa, digits=3, nsmall=3, justify="right", width=8))
    message("Iterations = ", format(iter, justify="right", width=8), "\n")
  }
  euddose<-EUD(dvh.matrix=dvh.matrix,a=aaa)     ## dose of patients
  # calculation of predicted probabilities
  # outcome MUST be a vector
  predicted<-matrix(nrow=length(outcome),ncol=2) 
  predicted[,2]<-outcome
  
  if (model=="Lyman") {
    predicted[,1]<-DR.Lyman(TD50=D50,gamma50=g50,aa=aaa,diffdvh=dvh.matrix)    
  }
  if (model=="Niemierko") {
    predicted[,1]<-DR.Niemierko(TD50=D50,gamma50=g50,aa=aaa,diffdvh=dvh.matrix)
  }

  # calculation of deviance http://http://stats.stackexchange.com/questions/6581/what-is-deviance-specifically-in-cart-rpart
  res<-outcome - predicted[,1]  # first calculates the residuals
  moddev<-t(res)%*%res          # calculates the deviance 
  # exclusion of values '1' that give NaN in LL function
  predicted[predicted[,1]!=1,]
  # initialize variables for C.I. calculation
  bmatrix<-NULL
  if (calcCI==TRUE) {
    if (C.I.type=="ProfLik") {
      # Confidence interval calculation with old procedure
      # outcome values rounded at 2nd decimal place
      d50L <- CI(MLE=MLE,arg="d50L",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
      d50U <- CI(MLE=MLE,arg="d50U",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
      g50L <- CI(MLE=MLE,arg="g50L",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
      g50U <- CI(MLE=MLE,arg="g50U",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
      aL   <- CI(MLE=MLE,arg="aL",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
      aU   <- CI(MLE=MLE,arg="aU",dvh.matrix=dvh.matrix,outcome=outcome,model=model,
                 min.TD50=min.TD50, max.TD50=max.TD50, min.gamma50=min.gamma50,
                 max.gamma50=max.gamma50, min.a=min.a, max.a=max.a)
    }
    if (C.I.type=="Boot") {
      # Confidence interval calculation using bootstrapping method
      temp.boot<-bootstrap.CI(model=model, TD50=D50, gamma50=g50, a=aaa, n=n.boot, 
                              dvh.matrix=dvh.matrix, C.I.width=C.I.width)
      d50L<-temp.boot$TD50L
      d50U<-temp.boot$TD50U
      g50L<-temp.boot$gamma50L
      g50U<-temp.boot$gamma50U
      aL  <-temp.boot$aL
      aU  <-temp.boot$aU
      # calculate max interation
      iter<-iter + temp.boot$TD50L[4] + temp.boot$TD50U[4] + 
        temp.boot$gamma50L[4] + temp.boot$gamma50U[4] +
        temp.boot$aL[4] + temp.boot$aU[4]
      # insert the bootstrap matrix
      bmatrix <- temp.boot$bootstrap.matrix
    }
  }
  
  # set NA values in CI if CI are not calculated
  if (calcCI==FALSE) {
    d50L <- c(NA,NA,NA,NA)
    d50U <- c(NA,NA,NA,NA)
    g50L <- c(NA,NA,NA,NA)
    g50U <- c(NA,NA,NA,NA)
    aL <- c(NA,NA,NA,NA)
    aU <- c(NA,NA,NA,NA)
  }
  
  output <- list("TD50"=D50,"TD50L"=d50L,"TD50U"=d50U,
                 "gamma50"=g50,"gamma50L"=g50L,"gamma50U"=g50U,
                 "a"=aaa,"aL"=aL,"aU"=aU,"model"=model,"Iterations"=iter,
                 "MLE"=MLE,"aic"=modelAIC,"outcome"=outcome,"dose"=euddose,
                 "deviance"=moddev,"predicted"=predicted[,1],
                 "DVH.matrix"=dvh.matrix,"bootstrap.matrix"=bmatrix)
  class(output)<-"NTCP"
  # remove global objects form .GlobalEnv
  if (!is.null(highD50)) rm(highD50, lowD50, highg50, lowg50, higha, lowa, envir=.GlobalEnv)
  rm(D50, aaa, g50, MLE, envir=.GlobalEnv)
  return(output)
}

# light version of fit.NTCP it's a function for fitting
# DVH data without C.I., used for bootstrap C.I. calculation
fit.NTCP.CI <- function(model=c("Lyman","Niemierko"),dvh.matrix, outcome, 
                     min.TD50=0, max.TD50=150, min.gamma50=0, max.gamma50=15,
                     min.a=-100, max.a=100) {
  model<-match.arg(model)
  if (model=="Lyman") {
    out <- nlminb(c(40,2,2), probit.nLL.Lyman, lower=c(min.TD50,min.gamma50,min.a),
                  upper=c(max.TD50,max.gamma50,max.a), dvh.matrix=dvh.matrix, outcome=outcome)
  }
  if (model=="Niemierko") {
    out <- nlminb(c(40,2,2), logit.nLL.Niemierko, lower=c(min.TD50,min.gamma50,min.a),
                  upper=c(max.TD50,max.gamma50,max.a), dvh.matrix=dvh.matrix, outcome=outcome)
  }
  # generate output
  output <- list("TD50"=out$par[1], "gamma50"=out$par[2], "a"=out$par[3],
                 "Iterations"=out$iterations, "MLE"=-out$objective)
  return(output)
}

# function for printing basic outcome of fitted functions
print.NTCP<-function(fittedmodel) {
  cat("                           -- NTCP model estimated parameters --","\n\n")
  cat("Model type:",fittedmodel$model,"\n\n")
  cat("Parameter:             Value:                         95% C.I.","\n\n")
  cat("TD50:   ",format(x=c(fittedmodel$TD50,fittedmodel$TD50L[1],fittedmodel$TD50U[1]),
                     justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("g50:    ",format(x=c(fittedmodel$gamma50,fittedmodel$gamma50L[1],fittedmodel$gamma50U[1]),
                     justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("a:      ",format(x=c(fittedmodel$a,fittedmodel$aL[1],fittedmodel$aU[1]),
                     justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("n (1/a):",format(x=c(1/fittedmodel$a,1/fittedmodel$aL[1],1/fittedmodel$aU[1]),
                     justify="right",digits=3,nsmall=3,width=20),"\n")
}

# function that shows model detailed features and alternate parameters
# found during profile likelihood procedures for confidence intervals calculation
summary.NTCP<-function(fittedmodel) {
  cat("                           -- NTCP model estimated parameters --","\n")
  cat("                           -- with alternate parameter values --","\n")
  cat("                           -- corresponding to C.I. limits    --","\n\n\n")
  cat("Model type:",fittedmodel$model,"\n")
  cat("MLE:",fittedmodel$MLE,"\n")
  cat("AIC:",fittedmodel$aic,"\n")
  cat("Deviance:",fittedmodel$deviance,"\n")
  cat("Model iterations number:",fittedmodel$Iterations,"\n\n\n")

  cat("Parameter:             Value:                         95% C.I.","\n\n")
  cat("TD50:   ",format(x=c(fittedmodel$TD50,fittedmodel$TD50L[1],fittedmodel$TD50U[1]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding g50:           ",format(x=c(fittedmodel$TD50L[2],fittedmodel$TD50U[2]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding a:             ",format(x=c(fittedmodel$TD50L[3],fittedmodel$TD50U[3]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("corresponding n=(1/a):       ",format(x=c(1/fittedmodel$TD50L[3],1/fittedmodel$TD50U[3]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  it<-fittedmodel$TD50L[4]+fittedmodel$TD50U[4]
  cat("C.I. Iterations number: ",it,"\n\n")

  cat("g50:    ",format(x=c(fittedmodel$gamma50,fittedmodel$gamma50L[1],fittedmodel$gamma50U[1]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding TD50:          ",format(x=c(fittedmodel$gamma50L[2],fittedmodel$gamma50U[2]),
                                             justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding a:             ",format(x=c(fittedmodel$gamma50L[3],fittedmodel$gamma50U[3]),
                                             justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("corresponding n=(1/a):       ",format(x=c(1/fittedmodel$gamma50L[3],1/fittedmodel$gamma50U[3]),
                                             justify="right",digits=3,nsmall=3,width=20),"\n")
  it<-fittedmodel$gamma50L[4]+fittedmodel$gamma50U[4]
  cat("C.I. Iterations number: ",it,"\n\n")
  
  cat("a:      ",format(x=c(fittedmodel$a,fittedmodel$aL[1],fittedmodel$aU[1]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("n=(1/a):",format(x=c(1/fittedmodel$a,1/fittedmodel$aL[1],1/fittedmodel$aU[1]),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding TD50:          ",format(x=c(fittedmodel$aL[2],fittedmodel$aU[2]),
                                             justify="right",digits=3,nsmall=3,width=20),"\n")
  cat("Corresponding gamma50:       ",format(x=c(fittedmodel$aL[3],fittedmodel$aU[3]),
                                             justify="right",digits=3,nsmall=3,width=20),"\n")

  it<-fittedmodel$aL[4]+fittedmodel$aU[4]
  cat("C.I. Iterations number: ",it,"\n\n")  
}

# function for predicting outcome in a given NTCP model
# with a specific level of dose
outcome.NTCP <- function(fittedmodel, dose) {
  x <- dose
  if (fittedmodel$model=="Lyman") {
    # y values of min C.I. curve
    ymin <- min(pnorm(q=((x - fittedmodel$TD50L[1])*fittedmodel$TD50L[2]*sqrt(2*pi))/fittedmodel$TD50L[1]),
                 pnorm(q=((x - fittedmodel$TD50U[1])*fittedmodel$TD50U[2]*sqrt(2*pi))/fittedmodel$TD50U[1]),
                 pnorm(q=((x - fittedmodel$gamma50L[2])*fittedmodel$gamma50L[1]*sqrt(2*pi))/fittedmodel$gamma50L[2]),
                 pnorm(q=((x - fittedmodel$gamma50U[2])*fittedmodel$gamma50U[1]*sqrt(2*pi))/fittedmodel$gamma50U[2]),
                 pnorm(q=((x - fittedmodel$aL[2])*fittedmodel$aL[3]*sqrt(2*pi))/fittedmodel$aL[2]),
                 pnorm(q=((x - fittedmodel$aU[2])*fittedmodel$aU[3]*sqrt(2*pi))/fittedmodel$aU[2]),
                 pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi))/fittedmodel$TD50))
    # y values of max C.I. curve
    ymax <- max(pnorm(q=((x - fittedmodel$TD50L[1])*fittedmodel$TD50L[2]*sqrt(2*pi))/fittedmodel$TD50L[1]),
                 pnorm(q=((x - fittedmodel$TD50U[1])*fittedmodel$TD50U[2]*sqrt(2*pi))/fittedmodel$TD50U[1]),
                 pnorm(q=((x - fittedmodel$gamma50L[2])*fittedmodel$gamma50L[1]*sqrt(2*pi))/fittedmodel$gamma50L[2]),
                 pnorm(q=((x - fittedmodel$gamma50U[2])*fittedmodel$gamma50U[1]*sqrt(2*pi))/fittedmodel$gamma50U[2]),
                 pnorm(q=((x - fittedmodel$aL[2])*fittedmodel$aL[3]*sqrt(2*pi))/fittedmodel$aL[2]),
                 pnorm(q=((x - fittedmodel$aU[2])*fittedmodel$aU[3]*sqrt(2*pi))/fittedmodel$aU[2]),
                 pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi))/fittedmodel$TD50))
    # predicted outcome with optimal parameters
    y <- pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi)/fittedmodel$TD50))
  }
  if (fittedmodel$model=="Niemierko") {
    # y values of min C.I. curve
    ymin <- min(1/(1+(fittedmodel$TD50L[1]/x)^(4*fittedmodel$TD50L[2])),
                 1/(1+(fittedmodel$TD50U[1]/x)^(4*fittedmodel$TD50U[2])),
                 1/(1+(fittedmodel$gamma50L[2]/x)^(4*fittedmodel$gamma50L[1])),
                 1/(1+(fittedmodel$gamma50U[2]/x)^(4*fittedmodel$gamma50U[1])),
                 1/(1+(fittedmodel$aL[2]/x)^(4*fittedmodel$aL[3])),
                 1/(1+(fittedmodel$aU[2]/x)^(4*fittedmodel$aU[3])),
                 1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)))
    # y values of max C.I. curve
    ymax <- max(1/(1+(fittedmodel$TD50L[1]/x)^(4*fittedmodel$TD50L[2])),
                 1/(1+(fittedmodel$TD50U[1]/x)^(4*fittedmodel$TD50U[2])),
                 1/(1+(fittedmodel$gamma50L[2]/x)^(4*fittedmodel$gamma50L[1])),
                 1/(1+(fittedmodel$gamma50U[2]/x)^(4*fittedmodel$gamma50U[1])),
                 1/(1+(fittedmodel$aL[2]/x)^(4*fittedmodel$aL[3])),
                 1/(1+(fittedmodel$aU[2]/x)^(4*fittedmodel$aU[3])),
                 1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)))
    #predicted outcome with optimal parameters
    y <- 1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50))
  }  
  cat("Dose:          Probability:                     95% C.I.","\n\n")
  cat(x, format(x=c(y, ymin, ymax),
                        justify="right",digits=3,nsmall=3,width=20),"\n")
}

# function for plotting NTCP curves:
# cases=TRUE: shows all cases on horizontal axes NTCP=0 and NTCP=1
# C.I.=TRUE: shows the confidence interval of the NTCP curve
# quantiles=FALSE: creates quantiles of cases for showing dots along the NTCP curve
# quantiles.no: quantiles number to be shown
# C.I.quantiles: shows the C.I. for NTCP of given quantiles
# C.I.width: width of confidence interval for quatiles
plot.NTCP <- function(fittedmodel,cases=TRUE,C.I.=TRUE,quantiles=FALSE,xmin=0,xmax=0,
                      quantiles.no=4,C.I.quantiles=FALSE,C.I.width=0.95,bounds=FALSE,...) {  
  # plot cases in the graph
  if (cases==TRUE){
    if (xmax==0) {
      plot(x=fittedmodel$dose,y=fittedmodel$outcome,xlim=c(xmin,max(fittedmodel$dose)+5),xlab="Dose [Gy]",ylab="NTCP [*100%]")
    } 
    else {
      plot(x=fittedmodel$dose,y=fittedmodel$outcome,xlim=c(xmin,xmax),xlab="Dose [Gy]",ylab="NTCP [*100%]")  
    }
  }
  # plot Lyman model
  if (fittedmodel$model=="Lyman") {
    if (xmax==0) {
      curve(expr=pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi)/fittedmodel$TD50)),from=xmin,to=max(fittedmodel$dose)+5,ylim=c(0,1),
            xlab="Dose [Gy]",ylab="NTCP [*100%]",add=cases,lwd=2)
    } else {
      curve(expr=pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi)/fittedmodel$TD50)),from=xmin,to=xmax,ylim=c(0,1),
            xlab="Dose [Gy]",ylab="NTCP [*100%]",add=cases,lwd=2)  
    }
  }
  # plot Niemierko model
  if (fittedmodel$model=="Niemierko") {
    if (xmax==0) {
      curve(expr=1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)),from=xmin,to=max(fittedmodel$dose)+5,ylim=c(0,1),
            xlab="Dose [Gy]",ylab="NTCP [*100%]",add=cases,lwd=2)
    } else {
      curve(expr=1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)),from=xmin,to=xmax,ylim=c(0,1),
            xlab="Dose [Gy]",ylab="NTCP [*100%]",add=cases,lwd=2)
    }
  }
  # plot confidence intervals of the models using the extreme values in the parameters matrix
  if (C.I.==TRUE){
    if (xmax==0) {
      x<-seq(0,max(fittedmodel$dose)+5,.25)
    } else {
      x<-seq(0,xmax,.25)
    }
    if (fittedmodel$model=="Lyman") {
      # y values of min C.I. curve
      ymin <- pmin(pnorm(q=((x - fittedmodel$TD50L[1])*fittedmodel$TD50L[2]*sqrt(2*pi))/fittedmodel$TD50L[1]),
                   pnorm(q=((x - fittedmodel$TD50U[1])*fittedmodel$TD50U[2]*sqrt(2*pi))/fittedmodel$TD50U[1]),
                   pnorm(q=((x - fittedmodel$gamma50L[2])*fittedmodel$gamma50L[1]*sqrt(2*pi))/fittedmodel$gamma50L[2]),
                   pnorm(q=((x - fittedmodel$gamma50U[2])*fittedmodel$gamma50U[1]*sqrt(2*pi))/fittedmodel$gamma50U[2]),
                   pnorm(q=((x - fittedmodel$aL[2])*fittedmodel$aL[3]*sqrt(2*pi))/fittedmodel$aL[2]),
                   pnorm(q=((x - fittedmodel$aU[2])*fittedmodel$aU[3]*sqrt(2*pi))/fittedmodel$aU[2]),
                   pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi))/fittedmodel$TD50))
      # y values of max C.I. curve
      ymax <- pmax(pnorm(q=((x - fittedmodel$TD50L[1])*fittedmodel$TD50L[2]*sqrt(2*pi))/fittedmodel$TD50L[1]),
                   pnorm(q=((x - fittedmodel$TD50U[1])*fittedmodel$TD50U[2]*sqrt(2*pi))/fittedmodel$TD50U[1]),
                   pnorm(q=((x - fittedmodel$gamma50L[2])*fittedmodel$gamma50L[1]*sqrt(2*pi))/fittedmodel$gamma50L[2]),
                   pnorm(q=((x - fittedmodel$gamma50U[2])*fittedmodel$gamma50U[1]*sqrt(2*pi))/fittedmodel$gamma50U[2]),
                   pnorm(q=((x - fittedmodel$aL[2])*fittedmodel$aL[3]*sqrt(2*pi))/fittedmodel$aL[2]),
                   pnorm(q=((x - fittedmodel$aU[2])*fittedmodel$aU[3]*sqrt(2*pi))/fittedmodel$aU[2]),
                   pnorm(q=((x - fittedmodel$TD50)*fittedmodel$gamma50*sqrt(2*pi))/fittedmodel$TD50))
    }
    if (fittedmodel$model=="Niemierko") {
      # y values of min C.I. curve
      ymin <- pmin(1/(1+(fittedmodel$TD50L[1]/x)^(4*fittedmodel$TD50L[2])),
                   1/(1+(fittedmodel$TD50U[1]/x)^(4*fittedmodel$TD50U[2])),
                   1/(1+(fittedmodel$gamma50L[2]/x)^(4*fittedmodel$gamma50L[1])),
                   1/(1+(fittedmodel$gamma50U[2]/x)^(4*fittedmodel$gamma50U[1])),
                   1/(1+(fittedmodel$aL[2]/x)^(4*fittedmodel$aL[3])),
                   1/(1+(fittedmodel$aU[2]/x)^(4*fittedmodel$aU[3])),
                   1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)))
      # y values of max C.I. curve
      ymax <- pmax(1/(1+(fittedmodel$TD50L[1]/x)^(4*fittedmodel$TD50L[2])),
                   1/(1+(fittedmodel$TD50U[1]/x)^(4*fittedmodel$TD50U[2])),
                   1/(1+(fittedmodel$gamma50L[2]/x)^(4*fittedmodel$gamma50L[1])),
                   1/(1+(fittedmodel$gamma50U[2]/x)^(4*fittedmodel$gamma50U[1])),
                   1/(1+(fittedmodel$aL[2]/x)^(4*fittedmodel$aL[3])),
                   1/(1+(fittedmodel$aU[2]/x)^(4*fittedmodel$aU[3])),
                   1/(1+(fittedmodel$TD50/x)^(4*fittedmodel$gamma50)))
    }    
    #lines(x,ymin,lty="dotted",lwd=2)
    #lines(x,ymax,lty="dotted",lwd=2)
    lines(x,ymin)
    lines(x,ymax)
  }
  # calculation of quantiles of cases in the plot for plotting quantiles
  if (quantiles==TRUE) {
    fq<-1/quantiles.no                    # fraction of each quantile
    meandoses<-c(1)                       # vector of mean doses
    length(meandoses)<-quantiles.no       # set length of vector of means
    meanoutcome<-c(1)                     # vector of mean outcomes
    length(meanoutcome)<-quantiles.no
    CImatrix<-matrix(nrow=quantiles.no,ncol=2)
    colnames(CImatrix)<-c(
      paste((1-C.I.width)/2,"%",sep=""),
      paste(C.I.width+(1-C.I.width)/2,"%",sep="")
    )
    series<-matrix(nrow=length(fittedmodel$dose), ncol=2)
    series[,1]<-fittedmodel$dose
    series[,2]<-fittedmodel$outcome
    for (m in 1:quantiles.no) {           # loop through quantiles
      doseV<-(1)                          # empty vector of mean doses
      outcomeV<-(1)                       # empty vector of mean outcome
      count<-0                            # counter of elements for each vector of means and results
      for  (n in 1:nrow(series)) {        # loop through all elements in the list        
        if (m==1){                        # 1st quantile
          if (series[n,1][[1]]<=quantile(x=fittedmodel$dose,fq)) {
            count<-count+1
            length(doseV)<-count
            length(outcomeV)<-count
            doseV[count]<-series[n,1][[1]]    # set the element of the dose vector
            outcomeV[count]<-series[n,2][[1]] # set the element of the outcome vector            
          }
        }
        if ((series[n,1][[1]]<=quantile(x=fittedmodel$dose,fq*m))&&(series[n,1][[1]]>quantile(x=fittedmodel$dose,fq*(m-1)))) {
          count<-count+1
          length(doseV)<-count
          length(outcomeV)<-count
          doseV[count]<-series[n,1][[1]]    # set the element of the dose vector
          outcomeV[count]<-series[n,2][[1]] # set the element of the outcome vector
        }          
      }
      meandoses[m]<-mean(doseV)        # set the element of vector of mean doses
      meanoutcome[m]<-mean(outcomeV)   # set the element of vector of mean outcomes
      bstrap <- c(1)                   # vector of bootstrap series for outcome C.I. calculation
      for (i in 1:1000) {
        bsample <- sample(outcomeV,length(outcomeV),replace=T) # sample generation
        bestimate <- mean(bsample)                             # bootstrapped sample mean
        bstrap <- c(bstrap,bestimate)
      }
      CImatrix[m,1]<-quantile(bstrap,(1-C.I.width)/2)          # lower bound of C.I. for mean
      CImatrix[m,2]<-quantile(bstrap,1-(1-C.I.width)/2)        # higher bound of C.I. for mean
    }
    for (p in 1:quantiles.no) {        # plot patients quantiles
      points(x=meandoses,y=meanoutcome,pch=15,cex=1.5)
    }
    if (C.I.quantiles==TRUE) {         # plot outcome C.I. for quantiles
      for (p in 1:quantiles.no) {
        lines(x=c(meandoses[p],meandoses[p]),y=c(CImatrix[p,1],CImatrix[p,2]))      # lines of C.I.
        lines(x=c(meandoses[p]-1,meandoses[p]+1),y=c(CImatrix[p,1],CImatrix[p,1]))  # lower bracket of C.I.
        lines(x=c(meandoses[p]-1,meandoses[p]+1),y=c(CImatrix[p,2],CImatrix[p,2]))  # higher bracket of C.I.
      }
      list("quantiles number"=quantiles.no,"quantiles mean doses"=meandoses,
           "quantiles mean outcomes"=meanoutcome,"quantiles confidence intervals"=CImatrix)
    }
  }  
}

# function for plotting the 3D density kernel of NTCP values of a single differential DVH
plot3D.NTCP <- function(dvh=NULL, model=NULL) {
  label<-c("\n")
  err.state<-0  # no error in getting varibles
  if (sum(dvh[,2])!=1) {
    err.state<-1
    label<-c(label,"dvh doesn't seem to be a relative differential dvh","\n")
  }
  if ((is.null(dvh)) || (is.matrix(dvh)==FALSE)) {
    err.state<-1
    label<-c(label,"You MUST provide a valid differential dvh for calculating 3D plot","\n")
  }
  if (ncol(dvh)>2) {
    err.state<-1
    label<-c(label,"dvh must be a matrix with two columns, 1st as dose bins, 2nd as volume bins", "\n")
  }
  if (is.null(model)) {
    err.state<-1
    label<-c(label, "You MUST provide a valid NTCP model for calculating 3D plot", "\n")
  }
  if (err.state==1) stop(label)
  # calculate the column of EUDs and NTCPs
  doses<-c()
  NTCPs<-c()
  if (model$model=="Lyman")
    for (n in 1:nrow(model$bootstrap.matrix)) {
      doses<-c(doses, EUD(dvh.matrix=dvh, a=model$bootstrap.matrix$a[n]))
      NTCPs<-c(NTCPs, DR.Lyman(TD50=model$bootstrap.matrix$TD50[n], 
                                 gamma50=model$bootstrap.matrix$gamma50[n], 
                                 dose=doses[n]))
    }
  if (model$model=="Niemierko") 
    for (n in 1:nrow(model$bootstrap.matrix)) {
      doses<-c(doses, EUD(dvh.matrix=dvh, a=model$bootstrap.matrix$a[n]))
      NTCPs<-c(NTCPs, DR.Niemierko(TD50=model$bootstrap.matrix$TD50[n], 
                                     gamma50=model$bootstrap.matrix$gamma50[n], 
                                     dose=doses[n]))
    }
  # creates 3D density kernel
  require(MASS)
  # density estimate
  k<-kde2d(x=doses, y=NTCPs, n=100)
  # plot perspective graph
  persp(k, box=TRUE, xlab="EUD (Gy)", ylab="NTCP", r=30, theta=135, phi=45, ticktype="detailed")  
  # calculate final values for EUD and NTCP for optimized parameters
  final.EUD<-EUD(dvh.matrix=dvh, a=model$a)
  if (model$model=="Lyman") final.NTCP<-DR.Lyman(TD50=model$TD50, gamma50=model$gamma50, aa=model$a, diffdvh=dvh)
  if (model$model=="Niemierko") final.NTCP<-DR.Niemierko(TD50=model$TD50, gamma50=model$gamma50, aa=model$a, diffdvh=dvh)
  points(final.EUD, y=final.NTCP, pch=20, col="red")
  plot(doses, NTCPs, pch="+")    
  points(final.EUD, y=final.NTCP, pch=20, col="red")
  abline(v=EUD(dvh.matrix=dvh, a=model$a), col="red", lty=2)
  abline(h=DR.Lyman(TD50=model$TD50, gamma50=model$gamma50, aa=model$a, diffdvh=dvh), col="red", lty=2)
  # calculates the filled contour of the kerneldensity  matrix
  # first calculate the volume of the 3D surface and normalizes to 1
  deltaX<-k$x[2]-k$x[1]
  deltaY<-k$y[2]-k$y[1]
  z.mat<-deltaX*deltaY*k$z       # matrix of volumes of bins of dose/NTCP
  total.abs.volume<-sum(z.mat)   # absolute total volume under the surface (3D integral)
  filled.contour(k)
  return(list("kernel.matrix"=k,"dose"=doses,"NTCP"=NTCPs))
}


# Log-Likelihood for Lyman TD50
probit.nLL.Lyman.TD50 <- function(z,C.I.TD50,dvh.matrix,outcome) { 
  gamma50 <- z[1]
  a <- z[2]
  pp <- pnorm(gamma50*sqrt(2*pi)*(EUD(dvh.matrix, a)/C.I.TD50-1))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# Log-Likelihood for Niemierko TD50
logit.nLL.Niemierko.TD50 <- function(z,C.I.TD50,dvh.matrix,outcome) { 
  gamma50 <- z[1]
  a <- z[2]
  pp <- 1/(1+(C.I.TD50/EUD(dvh.matrix, a))^(4*gamma50))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# Log-Likelihood for Lyman gamma50
probit.nLL.Lyman.gamma50 <- function(z,C.I.gamma50,dvh.matrix,outcome) { 
  TD50 <- z[1]
  a <- z[2]
  pp <- pnorm(C.I.gamma50*sqrt(2*pi)*(EUD(dvh.matrix, a)/TD50-1))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# Log-Likelihood for Niemierko gamma50
logit.nLL.Niemierko.gamma50 <- function(z,C.I.gamma50,dvh.matrix,outcome) { 
  TD50 <- z[1]
  a <- z[2]
  pp <- 1/(1+(TD50/EUD(dvh.matrix, a))^(4*C.I.gamma50))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# Log-Likelihood for Lyman a
probit.nLL.Lyman.a <- function(z,C.I.a,dvh.matrix,outcome) {
  TD50 <- z[1]
  gamma50 <- z[2]
  pp <- pnorm(gamma50*sqrt(2*pi)*(EUD(dvh.matrix, C.I.a)/TD50 - 1))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# Log-Likelihood for Niemierko a
logit.nLL.Niemierko.a <- function(z,C.I.a,dvh.matrix,outcome) {
  TD50 <- z[1]
  gamma50 <- z[2]
  pp <- 1/(1+(TD50/EUD(dvh.matrix, C.I.a))^(4*gamma50))
  p <- c(1:length(pp))
  m<-length(pp)
  for (n in 1:m) {
    if (pp[n]==1) p[n]<-0 else p[n]<-outcome[n]*log(pp[n])+(1-outcome[n])*log(1-pp[n])
  }
  return(-sum(p))
}

# building bounds of 95% Confidence Interval: profile likelihood
CI <- function(MLE, arg, dvh.matrix, outcome, model,
               min.TD50=0, max.TD50=150, min.gamma50=0, 
               max.gamma50=10, min.a=-100, max.a=100, verbose=TRUE) {
  if (arg=="d50L") {
    lab <-"lower TD50"
    par0<-"   Parameter Lower TD50 = "
    par1<-"  Corresponding gamma50 = "
    par2<-"        Corresponding a = "
  }
  if (arg=="d50U") {
    lab <-"upper TD50"
    par0<-"   Parameter Upper TD50 = "
    par1<-"  Corresponding gamma50 = "
    par2<-"        Corresponding a = "
  }
  if (arg=="g50L") {
    lab<- "lower gamma50"
    par0<-"Parameter Lower gamma50 = "
    par1<-"     Corresponding TD50 = "
    par2<-"        Corresponding a = "
  }
  if (arg=="g50U") {
    lab<- "upper gamma50"
    par0<-"Parameter Upper gamma50 = "
    par1<-"     Corresponding TD50 = "
    par2<-"        Corresponding a = "
  }
  if (arg=="aL") {
    lab<- "lower a"
    par0<-"      Parameter Lower a = "
    par1<-"     Corresponding TD50 = "
    par2<-"  Corresponding gamma50 = "
    
  }
  if (arg=="aU") {
    lab<- "upper a"
    par0<-"      Parameter Upper a = "
    par1<-"     Corresponding TD50 = "
    par2<-"  Corresponding gamma50 = "
  }
  if (verbose==TRUE) message("Calculating boundaries of Confidence Intervals: ", lab)
  it<-0 # number of iterations performed
  # set the model type
  if (model=="Lyman") {
    nLL.TD50 <- probit.nLL.Lyman.TD50
    nLL.gamma50 <- probit.nLL.Lyman.gamma50
    nLL.a <- probit.nLL.Lyman.a
  }
  if (model=="Niemierko") {
    nLL.TD50 <- logit.nLL.Niemierko.TD50
    nLL.gamma50 <- logit.nLL.Niemierko.gamma50
    nLL.a <- logit.nLL.Niemierko.a
  }
  # lower bound of TD50 C.I.
  if (arg=="d50L") {
    highD50 <- D50                 # starting point for interval search
    lowD50 <- D50 - 10              # ending point for interval search
    lowout <- nlminb(start=c(g50,aaa),objective=nLL.TD50,lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                     C.I.TD50=lowD50,dvh.matrix=dvh.matrix,outcome=outcome)
    it <- lowout$iterations
    while((-lowout$objective-MLE+ca1)>0) { # function to be repeated until LL-MLE+Chi2/2>0
      highD50 <- lowD50
      lowD50 <- lowD50 - 10
      lowout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                       C.I.TD50=lowD50,dvh.matrix=dvh.matrix,outcome=outcome)
      it <- it + lowout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    highD50 <<- highD50
    lowD50 <<- lowD50
    # start bisection search
    newD50<-(highD50+lowD50)/2
    newout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                     C.I.TD50=newD50,dvh.matrix=dvh.matrix,outcome=outcome)
    it <- it + newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      if (highD50-lowD50 < 1e-8) break
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        highD50 <- newD50
      } else lowD50 <- newD50
      newD50<-(highD50+lowD50)/2
      newout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                       C.I.TD50=newD50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations  
    }
    result <- (highD50+lowD50)/2 # final value of lower bound of D50 C.I. and iterations number
  }
  
  # higher bound of TD50 C.I.
  else if (arg=="d50U") {
    lowD50 <- D50                   # starting point for interval search
    highD50 <- D50 + 10             # ending point for interval search    
    highout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                      C.I.TD50=highD50,dvh.matrix=dvh.matrix,outcome=outcome) 
    it<-highout$iterations
    while((-highout$objective-MLE+ca1)>0) { # routine to be repeated until LL-MLE+Chi2/2>0
      lowD50 <- highD50
      highD50 <- highD50 + 10
      highout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                        C.I.TD50=highD50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+highout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    lowD50 <<- lowD50
    highD50 <<- highD50
    # start bisection search
    newD50<-(highD50+lowD50)/2
    newout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                     C.I.TD50=newD50,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-it+newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      if (highD50-lowD50 < 1e-8) break
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        lowD50 <- newD50
      } else highD50 <- newD50
      newD50<-(highD50+lowD50)/2
      newout <- nlminb(start=c(g50,aaa),objective=nLL.TD50, lower=c(min.gamma50,min.a), upper=c(max.gamma50,max.a),
                       C.I.TD50=newD50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations      
    }
    result <- (highD50+lowD50)/2 # final value of higher bound of D50 C.I.    
  }
  
  # lower bound of gamma50 C.I.
  else if (arg=="g50L") {
    highg50 <- g50                 # starting point for interval search
    lowg50 <- g50 - 0.5            # ending point for interval search
    lowout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                     C.I.gamma50=lowg50,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-lowout$iterations
    while((-lowout$objective-MLE+ca1)>0) { # routine to be repeated until LL-MLE+Chi2/2>0
      highg50 <- lowg50
      lowg50 <- lowg50 - 0.5
      lowout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                       C.I.gamma50=lowg50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+lowout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    highg50 <<- highg50
    lowg50 <<- lowg50
    # start bisection search
    newg50 <- (highg50+lowg50)/2
    newout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                     C.I.gamma50=newg50,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-it+newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        highg50 <- newg50
      } else lowg50 <- newg50
      newg50<-(highg50+lowg50)/2
      newout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                       C.I.gamma50=newg50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations      
    }
    result <- (highg50+lowg50)/2 # final value of lower bound of gamma50 C.I.
  }
  
  # higher bound of gamma50 C.I.
  else if (arg=="g50U") {
    lowg50 <- g50                   # starting point for interval search
    highg50 <- g50 + 0.5            # ending point for interval search
    highout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                      C.I.gamma50=highg50,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-highout$iterations
    while((-highout$objective-MLE+ca1)>0) { # routine to be repeated until LL-MLE+Chi2/2>0
      lowg50 <- highg50
      highg50 <- highg50 + 0.5
      highout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                        C.I.gamma50=highg50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+highout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    lowg50 <<- lowg50
    highg50 <<- highg50
    # start bisection search
    newg50<-(highg50+lowg50)/2
    newout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                     C.I.gamma50=newg50,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-it+newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        lowg50 <- newg50
      } else highg50 <- newg50
      newg50<-(highg50+lowg50)/2
      newout <- nlminb(start=c(D50,aaa),objective=nLL.gamma50, lower=c(min.TD50,min.a), upper=c(max.TD50,max.a),
                       C.I.gamma50=newg50,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations      
    }
    result <- (highg50+lowg50)/2 # final value of higher bound of D50 C.I.
  }
  
  # lower bound of a C.I.
  else if (arg=="aL") {
    if (aaa>5)  inc <- 1 else inc <- 0.25
    higha <- aaa                 # starting point for interval search
    lowa <- aaa - inc            # ending point for interval search
    lowout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                     C.I.a=lowa,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-lowout$iterations
    while((-lowout$objective-MLE+ca1)>0) { # routine to be repeated until LL-MLE+Chi2/2>0
      higha <- lowa
      lowa <- lowa - inc
      lowout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                       C.I.a=lowa,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+lowout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    higha <<- higha
    lowa <<- lowa
    # start bisection search
    newa <- (higha+lowa)/2
    newout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                     C.I.a=newa,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-it+newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        higha <- newa
      } else lowa <- newa
      newa<-(higha+lowa)/2
      newout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                       C.I.a=newa,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations      
    }
    result <- (higha+lowa)/2 # final value of lower bound of a C.I.
  }
  
  # higher bound of a C.I.
  else if (arg=="aU") {
    if (aaa>5)  inc <- 1 else inc <- 0.25
    lowa <- aaa                   # starting point for interval search
    higha <- aaa + inc            # ending point for interval search
    highout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                      C.I.a=higha,dvh.matrix=dvh.matrix,outcome=outcome) 
    it<-highout$iterations
    while((-highout$objective-MLE+ca1)>0) { # routine to be repeated until LL-MLE+Chi2/2>0
      lowa <- higha
      higha <- higha + inc
      highout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                        C.I.a=higha,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+highout$iterations
      if (it > 10000) stop("Failed convergence for C.I.")
    }
    # starting point for bisection
    lowa <<- lowa
    higha <<- higha
    # start bisection search
    newa<-(higha+lowa)/2
    newout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                     C.I.a=newa,dvh.matrix=dvh.matrix,outcome=outcome)
    it<-it+newout$iterations
    while (abs(-newout$objective-MLE+ca1)>1e-5) {
      # grid search of D50 value      
      if ((-newout$objective-MLE+ca1)>0) {   # moving around the zero
        lowa <- newa
      } else higha <- newa
      newa<-(higha+lowa)/2
      newout <- nlminb(start=c(D50,g50),objective=nLL.a, lower=c(min.TD50,min.gamma50), upper=c(max.TD50,max.gamma50),
                       C.I.a=newa,dvh.matrix=dvh.matrix,outcome=outcome)
      it<-it+newout$iterations      
    }
    result <- (higha+lowa)/2 # final value of higher bound of D50 C.I.
  }
  # vector of result for one side C.I. and corresponding other parameters and iteraions number
  bound <- c(result, newout$par[1], newout$par[2], it)
  # return one value of C.I. for each time and corresponding other parameters optimized
  # displays messages during C.I. calculation
  if (verbose==TRUE) {
    message(par0, format(result, digits=3, nsmall=3, justify="right", width=8))
    message(par1, format(newout$par[1], digits=3, nsmall=3, justify="right", width=8))
    message(par2, format(newout$par[2], digits=3, nsmall=3, justify="right", width=8))
    message("             Iterations = ", format(it, justify="right", width=8), "\n")
  }
  return(bound) 
}

# function for finding root of not linear equations using bisection method
bisect <- function(f, low, high, delta=1e-8) {
  a<-low
  b<-high
  d<-(a + b)/2
  it<-0
  errorbound <- (abs(b - a))/2
  if (f(a) == 0) return(list("root"=a, "iterations"=it))    
  if (f(b) == 0) return(list("root"=b, "iterations"=it))    
  while (errorbound > delta) {
    it<-it+1    
    if (f(d) == 0) return(list("root"=d, "iterations"=it))    
    if ((sign(f(a))*sign(f(d))) < 0) b<-d else a<-d
    d<-(a+b)/2
    errorbound<-errorbound/2    
  }
  if ((abs(d-low)<delta) ||(abs(b-high)<delta)) return(list("root"=NA, "iterations"=it)) else
    return(list("root"=d, "iterations"=it))
}


# function for calculating profile likelihood of DVH based model
proflik <- function(fitted.model, dvh.matrix, min.TD50=40, max.TD50=50,
                    min.gamma50=0.1, max.gamma50=5, min.a=-100, max.a=100, nbins=10) {
  if (fitted.model$model=="Lyman") {
    nLL.TD50 <- probit.nLL.Lyman.TD50
    nLL.gamma50 <- probit.nLL.Lyman.gamma50
    nLL.a <- probit.nLL.Lyman.a
  }
  if (fitted.model$model=="Niemierko") {
    nLL.TD50 <- logit.nLL.Niemierko.TD50
    nLL.gamma50 <- logit.nLL.Niemierko.gamma50
    nLL.a <- logit.nLL.Niemierko.a
  }
  message("Generating profile likelihood for TD50")  
  # create the matrix for TD50 profile likelihood
  proflik.TD50.matrix <- matrix(nrow=(nbins + 2), ncol=2)
  # create the vector of the TD50 and put into the matrix
  proflik.TD50.matrix[,1] <- sort(c(seq(from=min.TD50, to=max.TD50, by=((max.TD50 - min.TD50)/nbins)), fitted.model$TD50))
  # creates the progress bar
  pb <- txtProgressBar(min = 0, max = nbins + 2, style = 3)
  # creating the values of MLE for fixed parameter
  for (m in 1:nrow(proflik.TD50.matrix)) {
    setTxtProgressBar(pb, m)
    out<-nlminb(start=c(fitted.model$gamma50,fitted.model$a),objective=nLL.TD50, lower=c(0,0.01), upper=c(5,1000),
                C.I.TD50=proflik.TD50.matrix[m, 1],dvh.matrix=dvh.matrix,outcome=fitted.model$outcome)
    if (out$objective==0) proflik.TD50.matrix[m,2]<-NA else proflik.TD50.matrix[m,2]<- -out$objective
  }
  close(pb)  # closes the progress bar
  message("Generating profile likelihood for gamma50")  
  # create the matrix for gamma50 profile likelihood
  proflik.gamma50.matrix <- matrix(nrow=(nbins + 2), ncol=2)
  # create the vector of the TD50 and put into the matrix
  proflik.gamma50.matrix[,1] <- sort(c(seq(from=min.gamma50, to=max.gamma50, by=((max.gamma50 - min.gamma50)/nbins)), fitted.model$gamma50))
  # creates the progress bar
  pb <- txtProgressBar(min = 0, max = nbins + 2, style = 3)
  # creating the values of MLE for fixed parameter
  for (m in 1:nrow(proflik.gamma50.matrix)) {
    setTxtProgressBar(pb, m)
    out<-nlminb(start=c(fitted.model$TD50,fitted.model$a),objective=nLL.gamma50, lower=c(0,0.01), upper=c(100,1000),
                C.I.gamma50=proflik.gamma50.matrix[m, 1],dvh.matrix=dvh.matrix,outcome=fitted.model$outcome)
    if (out$objective==0) proflik.gamma50.matrix[m,2]<-NA else proflik.gamma50.matrix[m,2]<- -out$objective
  }
  close(pb)  # closes the progress bar
  message("Generating profile likelihood for a")
  # create the matrix for a profile likelihood
  proflik.a.matrix <- matrix(nrow=(nbins + 2), ncol=2)
  # create the vector of the TD50 and put into the matrix
  proflik.a.matrix[,1] <- sort(c(seq(from=min.a, to=max.a, by=((max.a - min.a)/nbins)), fitted.model$a))
  # creates the progress bar
  pb <- txtProgressBar(min = 0, max = nbins + 2, style = 3)
  # creating the values of MLE for fixed parameter
  for (m in 1:nrow(proflik.a.matrix)) {
    setTxtProgressBar(pb, m)
    out<-nlminb(start=c(fitted.model$TD50,fitted.model$gamma50),objective=nLL.a, lower=c(0,0), upper=c(100,5),
                C.I.a=proflik.a.matrix[m, 1],dvh.matrix=dvh.matrix,outcome=fitted.model$outcome)
    if (out$objective==0) proflik.a.matrix[m,2]<- NA else proflik.a.matrix[m,2]<- -out$objective
  }
  close(pb)  # closes the progress bar
  # output results
  return(list("TD50"=proflik.TD50.matrix, "gamma50"=proflik.gamma50.matrix, "a"=proflik.a.matrix))
}

# function for calculation of each element simulated in bootstrap during multicore calculation
calc.boot.el <- function(model, dvhnumber, maxdose, dosebin, TD50, gamma50, a) {
  boot.dvh.matrix <- gen.dvh(dvhnumber=dvhnumber, type="random", mindose=0, 
                             maxdose=maxdose, dosebin=dosebin, relative=TRUE)
  boot.outcome <- outcome.def(dvh.matrix=boot.dvh.matrix, TD50=TD50, gamma50=gamma50, 
                              a=a, model=model)
  boot.NTCP <- fit.NTCP.CI(model=model, dvh.matrix=boot.dvh.matrix, outcome=boot.outcome[,3])
  return(list("TD50"=boot.NTCP$TD50, "gamma50"=boot.NTCP$gamma50, "a"=boot.NTCP$a, "Iterations"=boot.NTCP$Iterations))
}

# function for calculation of C.I. of model using bootstrapping method
bootstrap.CI <- function(model, TD50, gamma50, a, n=1000, dvh.matrix, C.I.width=.95) {
  dvh.number<-ncol(dvh.matrix) - 1                          # find the number of DVH in the input DVH matrix 
  result.matrix<-matrix(nrow = n, ncol=4)                   # create the result matrix
  colnames(result.matrix)<-c("TD50", "gamma50", "a", "Iterations") # set the names of result matrix 
  max.dose<-max(dvh.matrix[,1])                             # set the maximum dose for bootstrapped DVHs
  dbin<-dvh.matrix[2,1]-dvh.matrix[1,1]                     # set the dose.bin for bootstrapped DVHs  
  if (Sys.info()["sysname"]=="Linux") {                     # check OS for using parallel computing under Linux
    require(foreach)
    require(doMC)
    require(parallel)
    registerDoMC(cores=detectCores())                       # compute the numbers of cores in CPU and register them
    boot.list <- foreach(m = 1:n) %dopar% calc.boot.el(model=model, dvhnumber=dvh.number, maxdose=max.dose, 
                                                       dosebin=dbin, TD50=TD50, gamma50=gamma50, a=a)
    for (p in 1:n) {
      result.matrix[p,1] <- boot.list[[p]]$TD50
      result.matrix[p,2] <- boot.list[[p]]$gamma50
      result.matrix[p,3] <- boot.list[[p]]$a
      result.matrix[p,4] <- boot.list[[p]]$Iterations
    }  
  } else {
    pb <- txtProgressBar(min = 0, max = n, style = 3)         # creates progressbar
    for (m in 1:n) {
      boot.dvh.matrix <- gen.dvh(dvhnumber=dvh.number, type="random", mindose=0, 
                                 maxdose=max.dose, dosebin=dbin, relative=TRUE)
      boot.outcome <- outcome.def(dvh.matrix=boot.dvh.matrix, TD50=TD50, gamma50=gamma50, 
                                  a=a, model=model)
      boot.NTCP <- fit.NTCP.CI(model=model, dvh.matrix=boot.dvh.matrix, outcome=boot.outcome[,3])
      # fills the value in the result matrix
      result.matrix[m,1] <- boot.NTCP$TD50
      result.matrix[m,2] <- boot.NTCP$gamma50
      result.matrix[m,3] <- boot.NTCP$a
      result.matrix[m,4] <- boot.NTCP$Iterations
      setTxtProgressBar(pb, m)
    }
    close(pb)
  }
  # output of confidence intervals by accessing into the result.matrix
  # find the bootstrapped series that is closest to the given C.I. values  
  TD50L<-which.min(abs(result.matrix[,1] - quantile(x=result.matrix[,1], probs=(1-C.I.width)/2)))
  TD50L<-c(as.numeric(result.matrix[TD50L,1]), as.numeric(result.matrix[TD50L,2]), 
                      as.numeric(result.matrix[TD50L,3]), as.numeric(result.matrix[TD50L,4]))
  TD50U<-which.min(abs(result.matrix[,1] - quantile(x=result.matrix[,1], probs=1-(1-C.I.width)/2)))
  TD50U<-c(as.numeric(result.matrix[TD50U,1]), as.numeric(result.matrix[TD50U,2]),
           as.numeric(result.matrix[TD50U,3]), as.numeric(result.matrix[TD50U,4]))
  gamma50L<-which.min(abs(result.matrix[,2] - quantile(x=result.matrix[,2], probs=(1-C.I.width)/2)))
  gamma50L<-c(as.numeric(result.matrix[gamma50L,2]), as.numeric(result.matrix[gamma50L,1]),
             as.numeric(result.matrix[gamma50L,3]), as.numeric(result.matrix[gamma50L,4]))
  gamma50U<-which.min(abs(result.matrix[,2] - quantile(x=result.matrix[,2], probs=1-(1-C.I.width)/2)))
  gamma50U<-c(as.numeric(result.matrix[gamma50U,2]), as.numeric(result.matrix[gamma50U,1]),
              as.numeric(result.matrix[gamma50U,3]), as.numeric(result.matrix[gamma50U,4]))
  aL<-which.min(abs(result.matrix[,3] - quantile(x=result.matrix[,3], probs=(1-C.I.width)/2)))
  aL<-c(as.numeric(result.matrix[aL,3]), as.numeric(result.matrix[aL,1]), 
                 as.numeric(result.matrix[aL,2]), as.numeric(result.matrix[aL,4]))
  aU<-which.min(abs(result.matrix[,3] - quantile(x=result.matrix[,3], probs=1-(1-C.I.width)/2)))
  aU<-c(as.numeric(result.matrix[aU,3]), as.numeric(result.matrix[aU,1]),
                 as.numeric(result.matrix[aU,2]), as.numeric(result.matrix[aU,4]))
  # generate output list
  output<-list("TD50L"=TD50L, "TD50U"=TD50U, "gamma50L"=gamma50L, "gamma50U"=gamma50U, 
               "aL"=aL, "aU"=aU, "bootstrap.matrix"=as.data.frame(result.matrix))
  return(output)
}


# function for calculating NTCP kernel from bootstrap matrix for a single DVH
NTCP.kernel <- function (model, dvh) {  
  doses <- c()   # empty vector of EUDs
  NTCP  <- c()   # empty vector of NTCP
  for (n in 1:nrow(model$bootstrap.matrix)) {
    doseV <- dvh[,1]^(model$bootstrap.matrix$a[n])  # calculates the doses
    addenda.EUD<-dvh[,2]*doseV
    doses<- c(doses, (sum(addenda.EUD))^(1/model$bootstrap.matrix$a[n])) # (apply(X=addenda.EUD,MARGIN=2,FUN=sum))^(1/a)
  }
  if (model$model == "Lyman") {
    for (n in 1:nrow(model$bootstrap.matrix))
    NTCP <- c(NTCP, DR.Lyman(TD50=model$bootstrap.matrix$TD50[n], 
                               gamma50=model$bootstrap.matrix$gamma50[n], 
                               dose=doses[n])) 
  }
  if (model$model == "Niemierko") {
    for (n in 1:nrow(model$bootstrap.matrix))
      NTCP <- c(NTCP, DR.Niemierko(TD50=model$bootstrap.matrix$TD50[n], 
                                 gamma50=model$bootstrap.matrix$gamma50[n], 
                                 dose=doses[n])) 
  }
  return(as.matrix(cbind(doses, NTCP)))
}

# function for fitting generic dose response function
fit.DR <- function(dvh, dose, outcome, start.TD50=50, start.gamma50=1.5, start.a=2,
                   DR.FUN=c("Lyman","Goitein","Niemierko","Munro","Okunieff","Warkentin","Bentzen"),
                   calcCI=TRUE, verbose=TRUE, C.I.width=.95, n.bins=200) {
  require(bbmle)
  # entry checks
  if (missingArg(dvh) && missingArg(dose)) stop("Please enter EITHER a dvhmatrix OR a dose vector")
  if (missingArg(outcome)) stop("You MUST provide an outcome vector for dose-response model fitting")
  if (!missingArg(dvh)) {
    if (class(dvh)!="dvhmatrix") warning("dvh isn't a dvhmatrix class object. If dvh isn't a relative differential dvh\n", 
                                         "  achieved results are wrong")
    if (class(dvh)=="dvhmatrix") dvh<-dvh@dvh
    if (ncol(dvh)!=(length(outcome)+1)) stop("Number of observation in dvhmatrix differs from outcome vector length")
  }
  fname<-DR.FUN
  DR.FUN=match.arg(arg=DR.FUN)
  DR.FUN=match.fun(FUN=paste("DR.", DR.FUN, sep=""))   # match the dose-response function
  # fits three parameters model
  if (missingArg(dose)) {
    if (verbose==TRUE) message("fitting 3 parameters ", fname, " model...")
    negLogLik<- function(TD50, gamma50, aa)      
      -sum(outcome*log(DR.FUN(diffdvh=dvh, TD50=TD50, gamma50=gamma50, aa=aa)) +
             (1-outcome)*log(1-DR.FUN(diffdvh=dvh, TD50=TD50, gamma50=gamma50, aa=aa)))
    model<-mle2(minuslogl=negLogLik, start=list(TD50=start.TD50, gamma50=start.gamma50, aa=start.a))
    predicted<-DR.FUN(diffdvh=dvh, TD50=model@coef[1], gamma50=model@coef[2], aa=model@coef[3])
    if (verbose==TRUE) message(fname, " model fitting successful")
    # calculate confidence interval of model
    if (calcCI==TRUE) {
      # calculates "fast" C.I. using confint function
      CI.model<-confint(model, level=C.I.width)
      # profiling Likelihood for a wider range of parameters
      # TD50 profiling
      low.TD50<-CI.model[1,1] - (CI.model[1,2]-CI.model[1,1])/20
      high.TD50<-CI.model[1,2] + (CI.model[1,2]-CI.model[1,1])/20
      TD50.seq<-seq(from=low.TD50, to=high.TD50, length.out=n.bins)    
      nLL.TD50.seq<-c()    
      negLogLikTD50<- function(gamma50, aa)      
        -sum(outcome*log(DR.FUN(diffdvh=model$dvh, TD50=TD50.seq[n], gamma50=gamma50, aa=aa)) +
               (1-outcome)*log(1-DR.FUN(diffdvh=model$dvh, TD50=TD50.seq[n], gamma50=gamma50, aa=aa)))
      cat("Profiling TD50...\n")
      pb <- txtProgressBar(min = 0, max = length(TD50.seq), style = 3)
      for (n in 1:length(TD50.seq)) {
        setTxtProgressBar(pb, n)
        nLL.TD50.seq<-c(nLL.TD50.seq, mle2(minuslogl=negLogLikTD50, start=list(gamma50=model$model@coef[2], aa=model$model@coef[3]))@min)
      }
      close(pb)
      # exploring right side of profile for TD50
      n<-which.min(nLL.TD50.seq)
      while (n<length(nLL.TD50.seq)) {
        if (((min(nLL.TD50.seq) + qchisq(p = .95, df = 1)/2 - nLL.TD50.seq[n])>0) &&  
               ((min(nLL.TD50.seq) + qchisq(p = .95, df = 1)/2 - nLL.TD50.seq[n + 1])<0)) {
          low<-n
          high<-n+1
          break
        }
      }
    }
  }
  # fits two parameters model
  if (!missingArg(dose)) {
    if (verbose==TRUE) message("fitting 2 parameters ", fname, " model...")
    negLogLik<- function(TD50, gamma50) 
      -sum(outcome*log(DR.FUN(dose=dose, TD50=TD50, gamma50=gamma50)) +
             (1-outcome)*log(1-DR.FUN(dose=dose, TD50=TD50, gamma50=gamma50)))
    model<-mle2(minuslogl=negLogLik, start=list(TD50=start.TD50, gamma50=start.gamma50))
    predicted<-DR.FUN(dose=dose, TD50=model@coef[1], gamma50=model@coef[2])
    if (verbose==TRUE) message(fname, " model fitting successful")
    # calculate confidence interval of model
    if (calcCI==TRUE) {
      # calculates "fast" C.I. using confint function
      CI.model<-confint(model, level=C.I.width)
      # start manual profiling
    }
  }
  # fills empty arguments
  if (missingArg(dose)) dose<-NULL
  if (missingArg(dvh))  dvh<-NULL
  # creates the output
  output <- list("model"=model, "CI.model"=CI.model, "Predicted"=predicted, 
                 "dvh.matrix"=dvh, "dose"=dose, "outcome"=outcome, "DR.FUN"=DR.FUN)
  # S3 class for DoseRespone
  class(output)<-"DoseResponse"
  return(output)  
}

# function for creating profile likelihood of a DoseResponse object
ProfLik.DR<-function(model, n.bins=200) {
  if (class(model)!="DoseResponse") stop("model must be a DoseResponse object")
  # profile likelihood for 3 parameters model
  DR.FUN<-model$DR.FUN
  outcome<-model$outcome
  if (is.null(model$dose)) {
    # profile likelihood of TD50
    low.TD50<-model$CI.model[1,1] - (model$CI.model[1,2]-model$CI.model[1,1])/20
    high.TD50<-model$CI.model[1,2] + (model$CI.model[1,2]-model$CI.model[1,1])/20
    TD50.seq<-seq(from=low.TD50, to=high.TD50, length.out=n.bins)    
    nLL.TD50.seq<-c()    
    negLogLikTD50<- function(gamma50, aa)      
      -sum(outcome*log(DR.FUN(diffdvh=model$dvh, TD50=TD50.seq[n], gamma50=gamma50, aa=aa)) +
             (1-outcome)*log(1-DR.FUN(diffdvh=model$dvh, TD50=TD50.seq[n], gamma50=gamma50, aa=aa)))
    cat("Profiling TD50...\n")
    pb <- txtProgressBar(min = 0, max = length(TD50.seq), style = 3)
    for (n in 1:length(TD50.seq)) {
      setTxtProgressBar(pb, n)
      nLL.TD50.seq<-c(nLL.TD50.seq, mle2(minuslogl=negLogLikTD50, start=list(gamma50=model$model@coef[2], aa=model$model@coef[3]))@min)
    }
    close(pb)
    # profile likelihood of gamma50
    low.gamma50<-model$CI.model[2,1] - (model$CI.model[2,2]-model$CI.model[2,1])/20
    high.gamma50<-model$CI.model[2,2] + (model$CI.model[2,2]-model$CI.model[2,1])/20
    gamma50.seq<-seq(from=low.gamma50, to=high.gamma50, length.out=n.bins)    
    nLL.gamma50.seq<-c()    
    negLogLikgamma50<- function(TD50, aa)      
      -sum(outcome*log(DR.FUN(diffdvh=model$dvh, gamma50=gamma50.seq[n], TD50=TD50, aa=aa)) +
             (1-outcome)*log(1-DR.FUN(diffdvh=model$dvh, gamma50=gamma50.seq[n], TD50=TD50, aa=aa)))
    cat("Profiling gamma50...\n")
    pb <- txtProgressBar(min = 0, max = length(gamma50.seq), style = 3)
    for (n in 1:length(gamma50.seq)) {
      setTxtProgressBar(pb, n)
      nLL.gamma50.seq<-c(nLL.gamma50.seq, mle2(minuslogl=negLogLikgamma50, start=list(TD50=model$model@coef[1], aa=model$model@coef[3]))@min)
    }
    close(pb)
    # profile likelihood of a
    low.a<-model$CI.model[3,1] - (model$CI.model[3,2]-model$CI.model[3,1])/20
    high.a<-model$CI.model[3,2] + (model$CI.model[3,2]-model$CI.model[3,1])/20
    a.seq<-seq(from=low.a, to=high.a, length.out=n.bins)
    nLL.a.seq<-c()    
    negLogLika<- function(TD50, gamma50)      
      -sum(outcome*log(DR.FUN(diffdvh=model$dvh, aa=a.seq[n], TD50=TD50, gamma50=gamma50)) +
             (1-outcome)*log(1-DR.FUN(diffdvh=model$dvh, aa=a.seq[n], TD50=TD50, gamma50=gamma50)))
    cat("Profiling a...\n")
    pb <- txtProgressBar(min = 0, max = length(a.seq), style = 3)
    for (n in 1:length(a.seq)) {
      setTxtProgressBar(pb, n)
      nLL.a.seq<-c(nLL.a.seq, mle2(minuslogl=negLogLika, start=list(TD50=model$model@coef[1], gamma50=model$model@coef[2]))@min)
    }
    close(pb)
  }
  return(list("TD50.ProfLik"=cbind(TD50.seq, nLL.TD50.seq), "gamma50.ProfLik"=cbind(gamma50.seq, nLL.gamma50.seq),
              "a.ProfLik"=cbind(a.seq, nLL.a.seq)))
}
robertogattabs/RadAgent documentation built on June 30, 2018, 12:02 a.m.