#' Function that creates a \code{dvhmatrix} class object
#'
#' @param dvh.number The number of DVHs to be generated
#' @param type The type of DVH to be generated: \code{random} generates all possible variation in dose distribution,
#' \code{convex} for DVH similar to PTV structures with good dose homogeneity to high levels,
#' \code{concave} for DVH with low level of dose delivered to the volume as a spared structure,
#' \code{mix} for DVH where dose distribution is averaged in the whole volume spectrum.
#' @param dvh.type The type of volume distribution to be created: \code{differential} or \code{cumulative}.
#' @param vol.distr Defines if the volume bins have to be divided by the total volume of the structure (\code{relative}) or not (\code{absolute}).
#' @param max.dose The upper bound of the dose distribution to be simulated.
#' @param dose.bin The dose bin in Gy to compute the value of volume parts in the final DVH.
#' @param volbin.side The value of the side of each cube (in mm) that builds the final volume of the simulated structure.
#' @param min.vol The minimum volume of the range to be simulated.
#' @param max.vol The maximum volume of the range to be simulated.
#' @description Creates an object of class \code{dvhmatrix}.
#' @details Dose Volume Histograms (DVH) are the basic representation of how the radiation doses distribute inside structures.
#' Usually they are shown in two forms: \code{differential} or \code{cumulative} being related each other because
#' \code{cumultive} DVHs are the integral form of \code{"differential"} ones. This function provides a siumlated series
#' of structures with the features defined by function parameters. The simulated series are useful for producing simulated
#' Dose/Response models using the modeling functions defined in this package.
#' class objects.
#' @references Van den Heuvela F. \emph{Decomposition analysis of differential dose volume histograms.} Med Phys. 2006 Feb;33(2):297-307. PubMed PMID: 16532934.
#' @return An object of \code{dvhmatrix} class.
#' @examples ## creates a dvhmatrix object with 200 differential - relative histograms
#' a<-DVH.generate(dvh.number=200, dvh.type="differential", vol.distr="relative")
#'
#' ## creates a dvhmatrix object containing 150 cumulative - absolute histograms
#' ## maximum dose in simulated series is 60 Gy
#' b<-DVH.generate(dvh.number=150, dvh.type="cumulative", vol.distr="absolute", max.dose=60)
#' @export
DVH.generate<-function(dvh.number, type=c("random","convex","concave","mix"),
dvh.type=c("differential", "cumulative"), vol.distr=c("relative", "absolute"),
max.dose = 75, dose.bin = 0.5, volbin.side = 2.5, min.vol=180, max.vol=220) {
if (min.vol<2) {
warning("Minimum volumes under 2 cc aren't allowed, min.vol set to 2.")
min.vol <- 2
}
if (max.dose <= 10) stop("max.dose must by higher than 10 Gy")
# create the vector of volumes
volumes <- runif(n = dvh.number, min = min.vol, max = max.vol)
volbin.num <- round(volumes/((volbin.side/10)^3)) # number of bins for each volume
# function for generating convex DVHs voxels series
convex.dvh <- function(n) {
mean.dose <- runif(n = 1, min = max.dose - (max.dose/3), max = max.dose - max.dose/6)
sd.dose <- (max.dose - mean.dose)/2.5
return(rnorm(n = n, mean = mean.dose, sd = sd.dose))
}
# function for generating concave DVHs voxels series
concave.dvh <- function(n) {
mean.dose <- runif(n = 1, min = 2, max = max(10, max.dose/2.5))
sd.dose <- mean.dose/3
result<-rnorm(n = n, mean = mean.dose, sd = sd.dose)
# takes only the voxels with dose>=0
result<-result[which(result>=0)]
# adds more voxels to reach the expected number with random uniform distribution
return(c(result, runif(n = n - length(result), min = 0, max = 2*mean.dose)))
}
# function for creating mix DVHs voxels series
mix.dvh <- function(n) {
contrib<-c(runif(n = 1, min = .05, max = .3), runif(n = 1, min = .05, max = .3))
# proportions in contributions to final DVH
contrib<-c(contrib[1], 1 - sum(contrib), contrib[2])
part1<-concave.dvh(n = round(contrib[1] * n))
part3<-convex.dvh(n = round(contrib[3] * n))
part2.mean<-runif(n = 1, min = max.dose/20, max = max.dose - max.dose/5)
part2.sd <-(max.dose - part2.mean)/(8 * 1/runif(n = 1, min = 1.5, max = 3.5))
part2<-rnorm(n = n - (length(part1) + length(part3)),
mean = part2.mean,
sd = part2.sd)
result<-c(part1, part2, part3)
result<-result[which(result>=0)]
# compensate the negative values when available
return(c(result, runif(n = n - length(result), min = 0, max = max.dose)))
}
# function that creates random DVHs according the previous three given functions
random.dvh <- function(n) {
FUN <- sample(x = c(convex.dvh, concave.dvh, mix.dvh), size = 1, replace = T, prob = c(.25,.25,.5))
return(FUN[[1]](n))
}
# voxels creation
type <- match.arg(type)
if (type=="convex") dose.voxels<-sapply(X = volbin.num, FUN = convex.dvh)
if (type=="concave") dose.voxels<-sapply(X = volbin.num, FUN = concave.dvh)
if (type=="mix") dose.voxels<-sapply(X = volbin.num, FUN = mix.dvh)
if (type=="random") dose.voxels<-sapply(X = volbin.num, FUN = random.dvh)
VolBin<-volbin.side^3/1000 # Volume Bin in cc
# delete dose.voxels with exceeding value over 15% max
dose.voxels<-sapply(X = dose.voxels, FUN = function(x) return(x[which(x<(max.dose+max.dose*.15))]), USE.NAMES = FALSE)
# creates the vector of structures volumes
if (dvh.number > 1) volume<-sapply(X = dose.voxels, FUN = function(x) VolBin * length(x), simplify = TRUE, USE.NAMES = FALSE) else volume <- VolBin * length(dose.voxels)
result<-new("dvhmatrix")
# creates the dvhmatrix object
dvh.type<-match.arg(arg = dvh.type)
vol.distr<-match.arg(arg = vol.distr)
result@dvh.type<-"differential" # default value, corrected by DVH.diff.to.cum if dvh.type = "cumulative"
result@vol.distr<-"absolute" # default DVH vol.distr type
result@volume<-volume
# creates the list of differential histograms
if (dvh.number > 1) {
hlist<-lapply(X = dose.voxels, FUN = hist,
breaks = seq(from = 0, to = max(c(unlist(lapply(X = dose.voxels, FUN = max))) + dose.bin * 4, max.dose), by = dose.bin), plot = FALSE)
if (dvh.type=="differential") {
result@dvh<-cbind(hlist[[1]]$mids, sapply(X = hlist, FUN = function(x) cbind(x$counts * VolBin), simplify = TRUE, USE.NAMES = FALSE))
#if (vol.distr=="relative") for (n in 2:ncol(result@dvh)) result@dvh[,n]<-result@dvh[,n]/result@volume[n-1]
}
if (dvh.type=="cumulative") {
result@dvh<-cbind(hlist[[1]]$mids, sapply(X = hlist, FUN = function(x) cbind(x$counts * VolBin), simplify = TRUE, USE.NAMES = FALSE))
result<-DVH.diff.to.cum(dvh = result)
#if (vol.distr=="relative") result<-DVH.diff.to.cum(dvh = result, relative = TRUE) else result<-DVH.diff.to.cum(dvh = result, relative = FALSE)
}
} else {
hlist<-hist(x = dose.voxels, breaks = seq(from = 0, to = max(max(dose.voxels) + dose.bin * 4, max.dose), by = dose.bin), plot = FALSE)
if (dvh.type=="differential") {
result@dvh <- cbind(hlist$mids, hlist$counts * VolBin)
#if (vol.distr=="relative") result@dvh[,2] <- result@dvh[,2]/result@volume
}
if (dvh.type=="cumulative") {
result@dvh <- cbind(hlist$mids, hlist$counts * VolBin)
result<-DVH.diff.to.cum(dvh = result)
#if (vol.distr=="relative") result<-DVH.diff.to.cum(dvh = result, relative = TRUE) else result<-DVH.diff.to.cum(dvh = result, relative = FALSE)
}
}
if (vol.distr == "relative") result<-DVH.relative(dvh = result)
return(result)
}
#' Function that converts differential DVHs into cumulative ones
#'
#' @param dvh Either an object of class \code{dvhmatrix} or \code{matrix} type.
#' @description Function that converts an object of class \code{dvhmatrix} or a simple \code{matrix} that
#' represents a differential DVH into a cumulative one.
#' @return Either an object of class \code{dvhmatrix} or a \code{matrix} according the argument \code{dvh}.
#' @export
DVH.diff.to.cum <- function(dvh) {
if ((!is.matrix(dvh))&&(class(dvh)!="dvhmatrix")) stop("dvh MUST be either an object of class dvhmatrix or a matrix")
if (class(dvh)=="dvhmatrix") dvh.matrix<-dvh@dvh else dvh.matrix<-dvh
if (class(dvh)=="dvhmatrix") if (dvh@dvh.type=="cumulative") return(dvh)
dvh.size <- dim(dvh.matrix)
DVHList <- matrix(nrow=dvh.size[1] + 1, ncol=dvh.size[2]) # create the matrix of cumulative DVHs
for (m in 2:dvh.size[2]) {
total.volume <- sum(dvh.matrix[,m])
for (n in 1:dvh.size[1]) # loop for rows
DVHList[n+1, m] <- total.volume - sum(dvh.matrix[c(1:n),m]) # elements of the matrix as relative volume
DVHList[1,m] <- total.volume # first element is total volume by default
}
DVHList[1,1]<-0
for (n in 2:nrow(DVHList)) DVHList[n,1]<-2*dvh.matrix[n-1,1]-DVHList[n-1,1] # set doses
if (class(dvh)=="dvhmatrix") {
dvh@dvh<-DVHList
dvh@dvh.type<-"cumulative"
return(dvh)
} else return(DVHList)
}
#' Function that converts cumulative DVHs into differential ones
#'
#' @param dvh Either an object of class \code{dvhmatrix} or \code{matrix} type.
#' @description Function that converts an object of class \code{dvhmatrix} or a simple \code{matrix} that
#' represents a cumulative DVH into a differential one.
#' @return Either an object of class \code{dvhmatrix} or a \code{matrix} according the argument \code{dvh}.
#' @export
DVH.cum.to.diff <- function(dvh) {
if ((!is.matrix(dvh))&&(class(dvh)!="dvhmatrix")) stop("dvh MUST be either an object of class dvhmatrix or a matrix")
if (class(dvh)=="dvhmatrix") dvh.matrix<-dvh@dvh else dvh.matrix<-dvh
if (class(dvh)=="dvhmatrix") if (dvh@dvh.type=="differential") return(dvh)
dvh.size <- dim(dvh.matrix)
DVHList <- matrix(nrow=dvh.size[1] - 1, ncol=dvh.size[2]) # create the matrix of differential DVHs
for (m in 2:dvh.size[2]) { # loop for columns (volumes)
total.volume<-dvh.matrix[1,m] # set total volume to 1st element of the cumulative DVH
for (n in 2:dvh.size[1]) # loop for rows
DVHList[n-1,m] <- (dvh.matrix[n-1,m]-dvh.matrix[n,m]) # differential volume for absolute distribution
}
for (m in 2:dvh.size[1]) { # loop for column of dose values
DVHList[m-1,1]<-dvh.matrix[m-1,1]+(dvh.matrix[m,1]-dvh.matrix[m-1,1])/2 # sets the dose values
}
if (class(dvh)=="dvhmatrix") {
dvh@dvh<-DVHList
dvh@dvh.type<-"differential"
return(dvh)
} else return(DVHList)
}
#' Extracts a DVH from a vector of dose bins
#' @param x vector of dose bins.
#' @param max.dose Upper dose bound for limiting the DVH computation.
#' @param dose.bin The dose bin for DVH computation in Gy.
#' @param dvh.type The type of volume distribution to be created: \code{differential} or \code{cumulative}.
#' @param vol.distr Defines if the volume bins have to be divided by the total volume of the structure (\code{relative}) or not (\code{absolute}).
#' @param createObj if \code{TRUE} returns a \code{dvhmatrix} class object.
#' @param volbin.side The value of the side of each cube (in cm) that builds the final volume of the structure.
#' @param voxel.volume The value of the volume of a single voxel (in cc) to compute the ROI total volume
#' @param total.volume The value of the total volume of the ROI (in cc)
#' @description Function that given a vector of dose bins (either a vector got from sampling a 3D mesh or a vector of
#' simple values of dose) extracts the DVH that summarizes that vector. Either \code{volbin.side}, \code{voxel.volume} or \code{total.volume} have to be
#' not \code{NULL} values. The priority list for computing the volume (if more than one argument is not \code{NULL}) is total.volume \code{->} voxel.volume \code{->} volbin.side.
#' @return Either an object of class \code{dvhmatrix} or a \code{matrix} according the argument \code{dvh} type.
#' @examples ## simulate a vector of dose bins
#' doses <- c(rnorm(n = 10000, mean = 45, sd = 3), rnorm(n = 7000, mean = 65, sd = 2.5))
#'
#' ## creates a dvhmatrix class object
#' DVH<-DVH.extract(x = doses)
#' @export
DVH.extract<-function(x, max.dose=NULL, dose.bin=.25, dvh.type=c("differential","cumulative"),vol.distr=c("relative","absolute"), createObj=TRUE, volbin.side=NULL, voxel.volume = NULL, total.volume = NULL) {
# default VolBin is given in cm3
if (!is.null(total.volume)) VolBin <- total.volume / length(x)
else if (!is.null(voxel.volume)) VolBin <- voxel.volume
else if (is.null(volbin.side)) stop('Either total.volume, voxel.volume or volbin.side MUST be a not NULL value')
else VolBin <- volbin.side^3
TotalVol<-length(x)*VolBin
dvh.type=match.arg(dvh.type)
vol.distr=match.arg(vol.distr)
if (is.null(x=max.dose)) max.dose<-max(x)+dose.bin*4
h<-hist(x=x, breaks=seq(from=0, to=max.dose,by=dose.bin), plot=FALSE)
diff<-cbind(h$mids, h$density/sum(h$density)*TotalVol)
# function for converting differential DVHs into
# relative differential DVHS
rel.diff.dvh <- function(dvh.matrix) {
dvh.size <- dim(dvh.matrix)
DVHList <- matrix(nrow=dvh.size[1], ncol=dvh.size[2]) # create the matrix of DVHs
for (m in 2:dvh.size[2]) { # loop for columns (volumes)
tot.volume <- sum(dvh.matrix[,m]) # calculate the total volume
for (n in 1:dvh.size[1]) { # loop for rows
DVHList[n, m] <- dvh.matrix[n, m]/tot.volume # elements of the matrix as relative volume
}
}
DVHList[,1]<-dvh.matrix[,1] # sets the dose column
return(DVHList)
}
# function for converting a matrix of differential DVH
# into matrix of cumulative DVH (relative volume)
cum.dvh <- function(dvh.matrix, relative=TRUE)
{
dvh.size <- dim(dvh.matrix)
DVHList <- matrix(nrow=dvh.size[1] + 1, ncol=dvh.size[2]) # create the matrix of cumulative DVHs
for (m in 2:dvh.size[2]) { # loop for columns (volumes)
tot.volume <- sum(dvh.matrix[,m]) # calculate the total volume
if (relative==TRUE) {
for (n in 1:dvh.size[1]) { # loop for rows
DVHList[n+1, m] <- (tot.volume - sum(dvh.matrix[c(1:n),m]))/tot.volume # elements of the matrix as relative volume
}
DVHList[1,m] <- 1 # first element is 1 by default if relative==TRUE
} else {
for (n in 1:dvh.size[1]) { # loop for rows
DVHList[n+1, m] <- tot.volume - sum(dvh.matrix[c(1:n),m]) # elements of the matrix as relative volume
}
DVHList[1,m] <- tot.volume # first element is total volume by default
}
}
DVHList[1,1]<-0
# DVHList[2:nrow(DVHList),1]<-dvh.matrix[,1]
for (n in 2:nrow(DVHList)) DVHList[n,1]<-2*dvh.matrix[n-1,1]-DVHList[n-1,1]
return(DVHList)
}
# return matrix without dvhmatrix class structure
if ((dvh.type=="differential") && (vol.distr=="absolute")) final.matrix<-diff
if ((dvh.type=="differential") && (vol.distr=="relative")) final.matrix<-rel.diff.dvh(diff)
if ((dvh.type=="cumulative") && (vol.distr=="absolute")) final.matrix<-cum.dvh(dvh.matrix=diff, relative=FALSE)
if ((dvh.type=="cumulative") && (vol.distr=="relative")) final.matrix<-cum.dvh(dvh.matrix=diff, relative=TRUE)
if (createObj==FALSE) return(final.matrix)
# return matrix within a dvhmatrix object
if (createObj==TRUE) return(new("dvhmatrix", dvh=final.matrix, dvh.type=dvh.type, vol.distr=vol.distr, volume=TotalVol))
}
#' mean of DVHs
#' @param dvh A \code{dvhmatrix} object
#' @description Function that gives the value of the mean dose of a \code{dvhmatrix} class object
#' @return a vector with the means of doses in the \code{dvhmatrix} object
#' @export
#' @examples ## generate a dataset of DVHs
#' a<-DVH.generate(dvh.number = 100)
#' m<-DVH.mean(a)
DVH.mean<-function(dvh) {
if (class(dvh)!="dvhmatrix") stop("dvh MUST be a dvhmatrix class object")
if (dvh@dvh.type=="cumulative") dvh<-DVH.cum.to.diff(dvh = dvh)
return(DVH.eud(dvh = dvh, a = 1)) # use a = 1 that corresponds to mean dose
}
#' Converts absolute \code{dvhmatrix} class objects to relative
#' @param dvh A \code{dvhmatrix} class object
#' @description Function that converts an object of \code{dvhmatrix} class where DVH are stored in absoulte mode into
#' another \code{dvhmatrix} class object where DVH are stored in relative mode.
#' @return A \code{dvhmatrix} class object.
#' @export
#' @examples ## generate a dataset of absolute DVHs
#' a<-DVH.generate(dvh.number = 10, vol.distr = "absolute")
#' b<-DVH.relative(dvh = a)
DVH.relative<-function(dvh) {
if (dvh@vol.distr=="absolute")
for (n in 1:(ncol(dvh@dvh) - 1)) dvh@dvh[,n+1]<-dvh@dvh[,n+1]/dvh@volume[n]
dvh@vol.distr<-"relative"
return(dvh)
}
#' Converts relative \code{dvhmatrix} class objects to absolute
#' @param dvh A \code{dvhmatrix} class object
#' @description Function that converts an object of \code{dvhmatrix} class where DVH are stored in relative mode into
#' another \code{dvhmatrix} class object where DVH are stored in absolute mode.
#' @return A \code{dvhmatrix} class object.
#' @export
#' @examples ## generate a dataset of relative DVHs
#' a<-DVH.generate(dvh.number = 10, vol.distr = "relative")
#' b<-DVH.absolute(dvh = a)
DVH.absolute<-function(dvh) {
if (dvh@vol.distr=="relative")
for (n in 1:(ncol(dvh@dvh) - 1)) dvh@dvh[,n+1]<-dvh@dvh[,n+1]*dvh@volume[n]
dvh@vol.distr<-"absolute"
return(dvh)
}
#' Calculates Equivalent Uniform Dose for a \code{dvhmatrix} object
#' @param dvh A \code{dvhmatrix} class object
#' @param a Value for parallel-serial correlation in radiobiological response
#' @description Function that calculates the value of Equivalent Uniform Dose (EUD) for a \code{dvhmatrix} object. The
#' formula to compute the equivalent uniform dose
#' is \deqn{EUD=\left (\sum_{i=1}^{ }D_{i}^{a}\cdot \frac{V_{i}}{V_{tot}} \right )^{1/a}}{EUD=(Sum(Di^a * Vi/Vtot))^(1/a)}
#' where \emph{i} is the counter of volume bins in the structure and \emph{a} is the parameter that relates
#' the organs response to the \emph{serial} or \emph{parallel} physiological function of the organ. Usually \emph{parallel}
#' structures are organs that can suffer depletion of functional subunits up to a given threshold, as liver, lungs or kidneys.
#' They are more sentitive to mean dose rather than hot spots and maximum dose levels delivered to the structure. On the other hand
#' \emph{serial} structures are organs more sensitive to hot-spots delivered in the volume, such as spinal cord.
#' @return A vector containing the values of EUD(s) for the given DVH(s)
#' @export
#' @references Niemierko A. \emph{Reporting and analyzing dose distributions: a concept of equivalent uniform dose.} Med Phys. 1997 Jan;24(1):103-10. PubMed PMID: 9029544.
#' @useDynLib moddicom
DVH.eud<-function(dvh, a = 1) {
dvh<-DVH.cum.to.diff(dvh = dvh)
dvh<-DVH.relative(dvh = dvh)
Ncol<-ncol(dvh@dvh) - 1
Nrow<-nrow(dvh@dvh)
ceud<-rep.int(x = 0, times = Ncol)
dosebin<-dvh@dvh[,1]
volumebin<-dvh@dvh[,2:(Ncol + 1)]
result<-.C("cEUD", as.double(dosebin), as.double(volumebin), as.double(a), as.integer(Nrow),
as.integer(Ncol), as.double(ceud))
return(result[[6]])
}
#' Merge two different \code{dvhmatrix} class objects into one
#' @param receiver The \code{dvhmatrix} object that will receive the \code{addendum} object.
#' @param addendum The \code{dvhmatrix} object to be merged with \code{addendum}.
#' @description This function can e used to create a \code{dvhmatrix} class object from two different
#' ones. The slots of the resulting object will be the same in the \code{receiver} one,
#' so some convertions can be automatically realized to create the homogeneous final result.
#' @return A \code{dvhmatrix} class object.
#' @export
#' @examples ## creates two different dvhmatrx objects
#' a<-DVH.generate(dvh.number = 100, dvh.type="differential", vol.distr = "relative")
#' b<-DVH.generate(dvh.number = 100, dvh.type="cumulative", vol.distr = "absolute")
#' ab<-DVH.merge(receiver = a, addendum = b)
DVH.merge<-function(receiver=NULL, addendum=NULL) {
# check dvhmatrix class
if ((class(receiver)!="dvhmatrix") || (class(addendum)!="dvhmatrix"))
stop("BOTH receiver AND addendum MUST be dvhmatrix class objects")
# creates list of dose bins
dbin<-list(receiver=receiver@dvh[3,1]-receiver@dvh[2,1], addendum=addendum@dvh[3,1]-addendum@dvh[2,1])
if (receiver@vol.distr=="relative") { # relative state of receiver
rel<-TRUE
vol.distr.fun.name <- "DVH.relative"
}
else {
rel<-FALSE
vol.distr.fun.name <- "DVH.absolute"
}
# STEP (0): check the relative/absolute distribution
if (receiver@vol.distr!=addendum@vol.distr) {
vol.distr.fun <- match.fun(FUN = vol.distr.fun.name)
addendum <- vol.distr.fun(addendum)
}
# STEP (1): identify conversion function if DVH types are different and convert addendum according receiver type
if (receiver@dvh.type!=addendum@dvh.type) {
if (receiver@dvh.type=="cumulative") conv.funct.name<-"DVH.diff.to.cum"
if (receiver@dvh.type=="differential") conv.funct.name<-"DVH.cum.to.diff"
conv.funct<-match.fun(FUN=conv.funct.name)
# correct the addendum if receiver has absolute distribution
addendum@dvh<-conv.funct(dvh=addendum@dvh)
}
# STEP (2): check the dose-bin for interpolating addendum or receiver according the lowest dbin
if (dbin$receiver!=dbin$addendum) {
warning("Different dose bins between receiver and addendum: linear interpolation performed.")
dose.seq<-seq(from = min(receiver@dvh[,1], addendum@dvh[,1]),
to = max(receiver@dvh[,1], addendum@dvh[,1]),
by = min(diff(x = receiver@dvh[,1]), diff(x = addendum@dvh[,1])))
if (ncol(receiver@dvh) == 2) { # single DVH receiver
apf.receiver <- approxfun(x = receiver@dvh[,1], y = receiver@dvh[,2])
temp.receiver <- cbind(dose.seq, apf.receiver(dose.seq))
} else { # multiple DVH receiver
apf.receiver <- apply(X = receiver@dvh[,2:ncol(receiver@dvh)], MARGIN = 2, FUN = approxfun, x = receiver@dvh[,1]) # generate list of approxfun
temp.receiver <- cbind(dose.seq, sapply(X = apf.receiver, FUN = function(x) return(x(dose.seq)), USE.NAMES = FALSE) ) # generate new dvh matrix
}
if (ncol(addendum@dvh) == 2) { # single DVH addendum
apf.addendum <- approxfun(x = addendum@dvh[,1], y = addendum@dvh[,2])
temp.addendum <- cbind(dose.seq, apf.addendum(dose.seq))
} else { # multiple DVH addendum
apf.addendum <- apply(X = addendum@dvh[,2:ncol(addendum@dvh)], MARGIN = 2, FUN = approxfun, x = addendum@dvh[,1]) # generate list of approxfun
temp.addendum <- cbind(dose.seq, sapply(X = apf.addendum, FUN = function(x) return(x(dose.seq)), USE.NAMES = FALSE) ) # generate new dvh matrix
}
temp.receiver[which(is.na(temp.receiver))]<-0
temp.addendum[which(is.na(temp.addendum))]<-0
} else {
temp.receiver<-receiver@dvh
temp.addendum<-addendum@dvh
}
# STEP (3): check if the two dvhs have equal length
if (nrow(temp.receiver)>nrow(temp.addendum)) {
for (n in 1:(nrow(temp.receiver)-nrow(temp.addendum))) temp.addendum<-rbind(temp.addendum, 0)
temp.addendum[,1]<-temp.receiver[,1]
}
if (nrow(temp.receiver)<nrow(temp.addendum)) {
for (n in 1:(nrow(temp.addendum)-nrow(temp.receiver))) temp.receiver<-rbind(temp.receiver, 0)
temp.receiver[,1]<-temp.addendum[,1]
}
# STEP (4): joins the two dvhs
result.DVH<-unname(cbind(temp.receiver, temp.addendum[,2:ncol(temp.addendum)]))
return(new("dvhmatrix", dvh=result.DVH, dvh.type=receiver@dvh.type, vol.distr=receiver@vol.distr,
volume=c(receiver@volume, addendum@volume)))
}
#' method for plotting the Log-Likelihood values in \code{DoseVolumeModel} class object
#' @description This method plots the results of fitting process for a \emph{Dose/Volume} model. The \emph{x} axis represents either the Dose,
#' or the Volume, the \emph{y} axis represents the Log-Likelihood result of the modeing process.
#' @param x The \code{DoseVolumeModel} class object
#' @param ShowFittingValue A logical value for showing the fitting point in the Log-Likelihood plot
#' @param ... Other parameters to be passed to \code{plot} function
#' @details The \code{DoseVolumeModel} class object contains the results of a modeling process across different values of Dose or
#' Volumes (depending on which kind of fitting was chosen during fitting process). The plot doesn't show the impact of other covariates
#' that could be part of the model itself.
#' @exportMethod plot
setMethod('plot', signature(x = 'DoseVolumeModel'),
function(x, ShowFittingValue = TRUE, ...) {
M <- slot(object = x, name = 'output.matrix')
tp <- slot(object = x, name = 'fitted.parameter')
value <- slot(object = x, name = 'fitted.value')
if (tp == 'Vdose') xlbl <- 'V-Dose [Gy]'
if (tp == 'Dvolume') xlbl <- 'D-Volume [cc]'
plot(x = M[,1], y = M[,2], xlab = xlbl, ylab = 'Log-Likelihood', type = 'l', ...)
if (ShowFittingValue == TRUE) {
abline(h = min(M[,2]), lty = 2)
abline(v = value, lty = 2)
points(x = value, y = min(M[,2]), ...)
axis(side = 1, at = round(x = value, digits = 2))
}
}
)
# method for class dvhmatrix
# enhance: is a vector of number of DVHs to be enhanced in the plot
# mean: plots the mean DVH overlapped to the background
# elements: if equal to a string "all" (default value) plots all the DVHs, if it is a vector plots only the
# chosen DVHs
#' method for plotting \code{dvhmatrix} class objects
#' @description This method overrides the default method \code{plot} and handles the plots of \code{dvhmatrix}
#' class objects.
#' @param x The \code{dvhmatrix} object to be plotted.
#' @param elements A vector containing the indeces of the DVHs to be plotted.
#' @param enhance A vector containing the indeces of the DVHs curves to be enhanced in the plot by changing de default color.
#' @param mean.dvh A \code{logical} value, if \code{TRUE} the mean dvh is plotted.
#' @param median.dvh A \code{logical} value, if \code{TRUE} the median dvh is plotted.
#' @param mean.median.alone A \code{logical} value, if \code{TRUE} only the mean and/or median DVH(s) is/are plotted.
#' Mean or median DVHs are plotted according the values in \code{mean.dvh} and \code{median.dvh}.
#' @param el.color The color for plotting DVHs. The default value is \code{"black"}.
#' @param en.color The color for plotting the enhanced DVHs. The default value is \code{"red"}.
#' @param mean.color The color for the mean dvh plot. The default value is \code{"blue"}.
#' @param median.color The color for the median dvh plot. The default value is \code{"red"}.
#' @param lwd The size of the lines in the plot. Default value is 1.
#' @param C.I.dvh Plot the confidence interval of mean and/or median DVH(s) according the values in \code{mean.dvh} and
#' \code{median.dvh}.
#' @param C.I.dvh.width The width of confidence interval of mean, median DVH(s) and overall dvh series.
#' @param C.I.dvh.range Plot the confidence interval of the whole DVHs series according the value in \code{C.I.dvh.width}.
#' @param n.boot The number of bootstrap repetitions for calculating the confidence interval of mean and median DVHs.
#' @param ... Other parameters to be passed to \code{plot} function.
#' @details The dvh are stored in the slot \code{dvh} of a \code{dvhmatrix} class object. Of course it is possible for7
#' users to use these matrices for plotting the dvhs considering that the structure of the matrix is referenced
#' under \code{\link{dvhmatrix-class}}. In order to simplify the plotting operations this method has been coded
#' allowing users to automatically obtain some useful features in the graphs, as the mean and median dvh curves,
#' and the confidence intervals for mean, median and overall dvh matrix. The confidence intervals of mean and
#' median are calculated by bootstrapping the whole \code{dvhmatrix} and calculating a given number of mean and
#' median dvh values that finally are reduced by using a quantile function to create the final plot.
#' @exportMethod plot
#' @examples ## create a dvhmatrix object
#' b <- DVH.generate(dvh.number = 100, dvh.type = "cumulative", vol.distr = "absolute")
#' plot(x = b, mean.dvh = TRUE, median.dvh = TRUE, mean.median.alone = TRUE,
#' C.I.dvh = TRUE, C.I.dvh.range = TRUE, C.I.dvh.width = .67, lwd = 2)
setMethod("plot", signature(x="dvhmatrix"),
function(x, elements=NULL, enhance=NULL, mean.dvh=FALSE, median.dvh=FALSE,
mean.median.alone=FALSE, el.color="black", en.color="red", mean.color="blue", median.color="red",
lwd=1, C.I.dvh=FALSE, C.I.dvh.width=.95, C.I.dvh.range=FALSE, n.boot=2000, ...){
dvh <- slot(x,"dvh")
dvh.type <- slot(x, "dvh.type")
vol.distr <- slot(x, "vol.distr")
# define the lables for the y axis
if ((dvh.type=="cumulative") && (vol.distr=="relative")) ylab<-"Volume [*100%]"
if ((dvh.type=="cumulative") && (vol.distr=="absolute")) ylab<-"Volume [cc]"
if ((dvh.type=="differential") && (vol.distr=="relative")) ylab<-"dVolume/dDose [*100%/Gy]"
if ((dvh.type=="differential") && (vol.distr=="absolute")) ylab<-"dVolume/dDose [cc/Gy]"
# max value of y axes
ymax<-max(dvh[,2:ncol(dvh)], na.rm = TRUE)
# direct plot for a single DVH
if (ncol(dvh)==2) plot(x=dvh[,1], y=dvh[,2], type="l", ylab=ylab, xlab="Dose [Gy]", lwd=lwd, ...)
# plot for more than one DVH
if ((ncol(dvh)>2) && (mean.median.alone==FALSE)) {
# creates the vector of elements to be plotted and enhanced
if (is.null(elements)) elements<-c(2:ncol(dvh)) else elements<-elements+1
if (!is.null(enhance)) for (n in 1:length(enhance)) elements<-elements[elements!=(enhance[n]+1)]
# open the plot panel and raws the 1st plot
if (length(elements)>0) plot(x=dvh[,1], y=dvh[,elements[1]], ylim=c(0, ymax), col=el.color,
lwd=lwd, type="l", xlab="Dose [Gy]", ylab=ylab, ...)
else
if (length(elements)>0) lines(x=dvh[,1], y=dvh[,elements[1]], ylim=c(0, ymax), col=el.color,
lwd=lwd)
# plots the remaining DVHs into "elements"
if (length(elements)==0) {
start.enhance<-2
plot(x=dvh[,1], y=dvh[,enhance[1]+1], col=en.color, lwd=lwd, type="l", ...)
} else start.enhance<-1
if (length(elements)>1)
for (n in 2:length(elements)) lines(x=dvh[,1], y=dvh[,elements[n]], col=el.color, lwd=lwd)
# plots DVHs to be enhanced
if (length(enhance)>0)
for (n in start.enhance:length(enhance)) lines(x=dvh[,1], y=dvh[,enhance[n]+1], col=en.color, lwd=lwd)
}
# plot the mean DVH
if (ncol(dvh)>2) {
mean.DVH<-apply(X=dvh[,2:ncol(dvh)], MARGIN=1, FUN=mean)
median.DVH<-apply(X=dvh[,2:ncol(dvh)], MARGIN=1, FUN=median)
}
if (mean.median.alone==FALSE) {
if ((mean.dvh==TRUE) && (ncol(dvh)>2))
lines(x=dvh[,1], y=mean.DVH, col=mean.color, lwd=lwd)
if ((median.dvh==TRUE) && (ncol(dvh)>2))
lines(x=dvh[,1], y=median.DVH , col=median.color, lwd=lwd)
}
if (mean.median.alone==TRUE) {
# plot both mean and median
if ((mean.dvh==TRUE) && (median.dvh==TRUE)) {
plot(x=dvh[,1], y=mean.DVH, ylim=c(0, ymax), col=mean.color,
lwd=lwd, type="l", xlab="Dose [Gy]", ylab=ylab, ...)
lines(x=dvh[,1], y=median.DVH , col=median.color, lwd=lwd)
}
# plot only mean
if ((mean.dvh==TRUE) && (median.dvh==FALSE)) {
plot(x=dvh[,1], y=mean.DVH, ylim=c(0, ymax), col=mean.color,
lwd=lwd, type="l", xlab="Dose [Gy]", ylab=ylab, ...)
}
# plot only median
if ((mean.dvh==FALSE) && (median.dvh==TRUE)) {
plot(x=dvh[,1], y=median.DVH, ylim=c(0, ymax), col=median.color,
lwd=lwd, type="l", xlab="Dose [Gy]", ylab=ylab, ...)
}
}
# plot of C.I. area, both 95% C.I. of the mean dose and 95% with quantiles are calculated
if ((C.I.dvh==TRUE) && (ncol(dvh)>2)) {
message("sampling DVHs for C.I. calculation...")
DVH.CI<-DVH.baseStat(dvh = x, C.I.width = C.I.dvh.width, n.boot = n.boot)
if ((mean.dvh==TRUE)&&(median.dvh==FALSE)){
lowerCI.mean <- DVH.CI@dvh[,3]
higherCI.mean<- DVH.CI@dvh[,4]
xpoly.mean<-c(dvh[,1], rev(dvh[,1])) # vector of x coordinates of polygon of CI
ypoly.mean<-c(lowerCI.mean, rev(higherCI.mean)) # vector of y coordinates of polygon of CI of mean
# draws the polygon of mean CI
polygon(x=xpoly.mean, y=ypoly.mean, col=adjustcolor(col=mean.color, alpha.f=.25))
}
if ((mean.dvh==FALSE)&&(median.dvh==TRUE)){
lowerCI.median <- DVH.CI@dvh[,6]
higherCI.median<- DVH.CI@dvh[,7]
xpoly.median<-c(dvh[,1], rev(dvh[,1])) # vector of x coordinates of polygon of CI
ypoly.median<-c(lowerCI.median, rev(higherCI.median)) # vector of y coordinates of polygon of CI of median
# draws the polygon of median CI
polygon(x=xpoly.median, y=ypoly.median, col=adjustcolor(col=median.color, alpha.f=.25))
}
if ((median.dvh==TRUE)&&(mean.dvh==TRUE)){
lowerCI.mean <- DVH.CI@dvh[,3]
higherCI.mean<- DVH.CI@dvh[,4]
xpoly.mean<-c(dvh[,1], rev(dvh[,1])) # vector of x coordinates of polygon of CI
ypoly.mean<-c(lowerCI.mean, rev(higherCI.mean)) # vector of y coordinates of polygon of CI of mean
lowerCI.median <- DVH.CI@dvh[,6]
higherCI.median<- DVH.CI@dvh[,7]
xpoly.median<-c(dvh[,1], rev(dvh[,1])) # vector of x coordinates of polygon of CI
ypoly.median<-c(lowerCI.median, rev(higherCI.median)) # vector of y coordinates of polygon of CI of median
# draws the polygon of mean CI
polygon(x=xpoly.mean, y=ypoly.mean, col=adjustcolor(col=mean.color, alpha.f=.25))
# draws the polygon of median CI
polygon(x=xpoly.median, y=ypoly.median, col=adjustcolor(col=median.color, alpha.f=.25))
}
}
# plot the C.I. of data range of DVHs
if (C.I.dvh.range==TRUE) {
# quantiles for DVHs data range
lowerCI <- c()
higherCI <- c()
for (n in 1:nrow(dvh)) {
lowerCI <- c(lowerCI, quantile(x=dvh[n,2:ncol(dvh)], probs=((1-C.I.dvh.width)/2)))
higherCI<- c(higherCI, quantile(x=dvh[n,2:ncol(dvh)], probs=((1+C.I.dvh.width)/2)))
}
xpoly<-c(dvh[,1], rev(dvh[,1])) # vector of x coordinates of polygon of CI
ypoly<-c(lowerCI, rev(higherCI)) # vector of y coordinates of polygon of CI
polygon(x=xpoly, y=ypoly, col=adjustcolor(col="grey", alpha.f=.25)) # draws the polygon
}
}
)
#' Function for extracting the V-Dose from cumulative DVH(s)
#' @description Function for calculating the value of V-Dose of a \code{dvhmatrix} object
#' @param dvh A \code{dvhmatrix} class object
#' @param Dose The dose for calculating the V-Dose value(s)
#' @details V-Dose is, in a given Dose Volume Histogram, the value of the volume of the structure receiving
#' \strong{at least} the chosen level of dose.
#' @return A vector containing the value(s) of V-Dose(s).
#' @export
#' @examples ## create a dvhmatrix class object
#' a<-DVH.generate(dvh.number = 100)
#' DVH.Vdose(dvh = a, Dose = 50)
DVH.Vdose <- function(dvh, Dose) {
dvh<-DVH.diff.to.cum(dvh = dvh)
dv<-dvh@dvh[,1] # vector of doses
if (ncol(dvh@dvh) > 2) {
apf <- apply(X = dvh@dvh[,2:ncol(dvh@dvh)], MARGIN = 2, FUN = approxfun, x = dv) # generate list of approxfun
return( sapply(X = apf, FUN = function(x) return(x(Dose))) )
} else { # option with 1 single DVH
apf <- approxfun(x = dv, y = dvh@dvh[,2])
return(apf(Dose))
}
}
#' Function for extracting the D-Volume from cumulative DVH(s)
#' @description Function for calculating the value of D-Volume of a \code{dvhmatrix} object
#' @param dvh A \code{dvhmatrix} class object
#' @param Volume Volume for calculating the D-Volume value(s). Default is 0.001, corresponding to about maximum dose
#' delivered to the structure.
#' @details D-Volume is, in a given Dose Volume Histogram, the value of the dose received
#' by a given volume of the structure(s) of interest.
#' @return A vector containing the value(s) of D-Volume(s).
#' @export
#' @examples ## create a dvhmatrix class object
#' a<-DVH.generate(dvh.number = 100)
#' DVH.Dvolume(dvh = a, Volume = 0.5)
DVH.Dvolume <- function(dvh, Volume=0.001) {
dvh<-DVH.diff.to.cum(dvh = dvh)
dv<-dvh@dvh[,1] # vector of doses
if (ncol(dvh@dvh) > 2) {
apf <- apply(X = dvh@dvh[,2:ncol(dvh@dvh)], MARGIN = 2, FUN = approxfun, y = dv) # generate list of approxfun
DVol <- sapply(X = apf, FUN = function(x) return(x(Volume)))
DVol[which(is.na(DVol))] <- 0
return(DVol)
} else { # option with 1 single DVH
apf <- approxfun(x = dvh@dvh[,2], y = dv)
DVol <- apf(Volume)
if(is.na(DVol)) DVol <- 0
return(DVol)
}
}
#' Function for calculating the mean and median dvh with related confidence intervals
#' @description This function calculates the mean dvh from a \code{dvhmatrix} class object. The mean dvh is
#' calculated with its confidence interval that is given by a bootstrapped dvh series from
#' the dvh given in the \code{dvh} object.
#' @param dvh A \code{dvhmatrix} class object
#' @param C.I.width The width of confidence interval
#' @param n.boot The number of bootstrapped dvhs for computing the quantile in C.I. calculation
#' @return A \code{dvhmatrix} object where 7 columns: Dose, mean DVH, low C.I. of mean DVH, high C.I. of mean DVH,
#' median DVH, low C.I. of median DVH, high C.I. of median DVH.
#' @export
#' @useDynLib moddicom
#' @examples ## creates a dvhmatrix class object with 100 DVHs
#' a<-DVH.generate(dvh.number = 100, dvh.type = "cumulative", vol.distr = "relative")
#' b<-DVH.baseStat(dvh = a)
#' ## plot the dvhmatrix object "b", showing both mean,
#' ## median DVH and corresponding confidence intarvals
#' plot(b)
DVH.baseStat<-function(dvh, C.I.width = .95, n.boot = 2000) {
# calculate mean dvh
mean.dvh<-apply(X = dvh@dvh[,2:ncol(dvh@dvh)], MARGIN = 1, FUN = mean)
median.dvh<-apply(X = dvh@dvh[,2:ncol(dvh@dvh)], MARGIN = 1, FUN = median)
# vectorized DVH
Vdvh<-as.vector(dvh@dvh[,2:ncol(dvh@dvh)])
# vector for sampled dvh
sampledvh<-rep.int(x = 0, times = nrow(x = dvh@dvh) * (ncol(x = dvh@dvh) - 1))
# create the mean dvh vector
meanV<-rep.int(x = 0, times = nrow(dvh@dvh) * n.boot)
medianV<-meanV
stepMedian<-rep.int(x = 0, times = ncol(dvh@dvh) - 1)
# number of histograms
Ndvh <- ncol(dvh@dvh) - 1
result<-(.C("meanmediandvh", as.double(Vdvh), as.integer(nrow(dvh@dvh)), as.integer(n.boot),
as.double(meanV), as.double(sampledvh), as.integer(Ndvh), as.double(medianV), as.double(stepMedian)))
# create dvh matrix
meanDvh<-matrix(data = result[[4]], nrow = nrow(dvh@dvh))
# generate the mean CI DVH
meanCI<-apply(X = meanDvh, MARGIN = 1, FUN = quantile, probs = c((1-C.I.width)/2, (1+C.I.width)/2))
# generate the median CI DVH
medianDvh<-matrix(data = result[[7]], nrow = nrow(dvh@dvh))
medianCI<-apply(X = medianDvh, MARGIN = 1, FUN = quantile, probs = c((1-C.I.width)/2, (1+C.I.width)/2))
# return(new("dvhmatrix", dvh = cbind(dvh@dvh[,1], medianDvh), dvh.type = dvh@dvh.type,
# vol.distr = dvh@vol.distr, volume = rep.int(x = mean(dvh@volume), times = n.boot)))
result.DVH<-new("dvhmatrix", dvh = cbind(dvh@dvh[,1], mean.dvh, meanCI[1,], meanCI[2,], median.dvh, medianCI[1,], medianCI[2,]), dvh.type = dvh@dvh.type,
vol.distr = dvh@vol.distr, volume = c(rep.int(x = mean(dvh@volume), times = 3),
rep.int(x = median(dvh@volume), times = 3)))
colnames(x = result.DVH@dvh)<-c("Dose [Gy]", "mean", paste((1-C.I.width)/2,"% C.I. of mean"), paste((1+C.I.width)/2,"% C.I. of mean"),
"median", paste((1-C.I.width)/2,"% C.I. of median"), paste((1+C.I.width)/2,"% C.I. of median"))
return(result.DVH)
}
#' Function that returns the linear quadratic correction of dose bin in a \code{dvhmatrix} object
#' @description This function appies the linear-quadratic correction for the dose in a \code{dvmatrix} object by
#' converting each dose bins as function of a given \eqn{\alpha\beta} ratio (default value is 3, valid for late outcome).
#' @param dvh A \code{dvhmatrix} object
#' @param ref.frac Reference fractionation, default value 2 for direct computing EQD2
#' @param nf Fractions number for the input dvh
#' @param alphabeta \eqn{\alpha\beta} ratio for computing LQ cvonversion, default value 3
#' @return A \code{dvhmatrix} object with converted dose bins
#' @export
DVH.lq.correct<-function(dvh, ref.frac = 2, nf, alphabeta = 3, reverse = TRUE) {
dvh@dvh[,1] <- dvh@dvh[,1] * (alphabeta + dvh@dvh[,1]/nf) / (alphabeta + ref.frac)
return(dvh)
}
#' Function for fitting Vdose and Dvolume to a given \code{dvhmatrix} object and outcome vector
#' @description This function can be used to find the best \emph{Vdose} or \emph{Dvolume} values that describe
#' a given clinically observed outcome
#' @details The \emph{Vdose} is the value of the volume of a structure in a dose volume histogram that receives
#' a given level of dose. It is a widely used indicator and predictor of possible outcome when DVHs
#' have to be clinically evaluated. There are many papers that exploit the correct use of \emph{Vdose} for
#' clinical evaluation of treatment plans. The \code{DR.fit.DoseVolume} function used inside the \pkg{moddicom} package
#' allows to fit clinical data with \code{dvhmatrix} objects in order to detect the dose-volume relationship achievable
#' by dosimetric data. The dose-volume-outcome fitting is performed by iterative search, so, as far as the author know,
#' this is the first example of dose-volume-outcome fitting function that allows to find out \emph{the best fitting value}
#' of \emph{Vdose} or \emph{Dvolume}, rather than using an empyrical search from pre-fetermined values as usually reported in literature.
#' When setting the \code{model} parameter the users have to chose between a \code{\link[stats]{glm}} model and a \code{\link[survival]{coxph}} model
#' by writing the formula of the model \emph{without} the \code{dvhmatrix} object
#' @param dvh A \code{dvhmatrix} object
#' @param outcome A vector of binary values (if not it will coerced to binary values) showing the observed outcome
#' @param model The model that will be used to fit the dose-volume-response data (see examples). Models with available
#' the \code{\link[stats]{logLik}} method can be used belonging to classes \code{lm}, \code{glm} and \code{coxph}
#' @param type A character value representing the two type of dose-volume fitting as \code{Vdose} or \code{Dvolume}
#' @param CI Logcal Value for calculating the confidence interval of the dose-volume fitting parameter and for the \eqn{\alpha\beta} if modLQ is optioned \code{TRUE}
#' @param CI.width Width of confidence interval
#' @param epsilon Error limit for calculating the Log Likelihood function in determining the confidence interval of dose-volume parameter
#' @examples ## NOT RUN DR.fit.DoseVolume(dvh = D, outcome = Dout, model = glm(formula = Dout ~ 1 , family = binomial(link = 'logit')), type = 'Vdose')
#' @references \emph{Special Considerations Regarding Absorbed-Dose and Dose-Volume Prescribing and Reporting in IMRT}. J ICRU. 2010 Apr;10(1):27-40. doi: 10.1093/jicru/ndq008. PubMed PMID: 24173325.
#' @references Graham MV, Purdy JA, Emami B, Harms W, Bosch W, Lockett MA, Perez CA. \emph{Clinical dose-volume histogram analysis for pneumonitis after 3D treatment for non-small cell lung cancer (NSCLC)}. Int J Radiat Oncol Biol Phys. 1999 Sep 1;45(2):323-9. PubMed PMID: 10487552.
#' @export
DR.fit.DoseVolume<-function(dvh, outcome, model, type = c("Vdose", "Dvolume"), CI = FALSE, CI.width = 0.95, epsilon = 1e-6) {
type <- match.arg(arg = type)
if (dvh@dvh.type == "differential") dvh<-DVH.diff.to.cum(dvh = dvh) # transform in cumulative if differential
if (type == "Vdose") { # create the approxfun and the fitting Vdose function (FUN)
value<-seq(from = 1, to = max(floor(dvh@dvh[,1])), by = .5) # create the dose vector
apf <- apply(X = dvh@dvh[, 2:ncol(dvh@dvh)], MARGIN = 2, FUN = approxfun, x = dvh@dvh[, 1]) # approxfun list
update.model<-function(model, Vdose) { # function for interactively update the model
Vdose <<- Vdose # .GlobalEnv needed for scoping of variable
new.model<-update(object = model, formula. = ". ~ . + Vdose")
return(new.model)
}
f <- function(val) return(sapply(X = apf, FUN = function(x) return(x(val)))) # function that calculates the Vdoses
output.matrix<-c()
for (n in 1:length(value)) {
Vd<-f(value[n])
temp.model<-update.model(model = model, Vdose = Vd)
coeff<-summary(temp.model)$coefficients[,4]
if ( length(coeff) == length(temp.model$coefficients) ) # check length to prevent joining models with only one covariate
output.matrix<-rbind(output.matrix, c(value[n], -logLik(temp.model), coeff)) else
warning(paste("Vdose =", value[n], "Gy doesn't fit a model"))
}
obj.FUN<-function(x) { # objective function for optimizing the best model
Vdose<<-f(x[1])
new.model<-update(object = model, formula. = ". ~ . + Vdose")
return(-logLik(new.model))
}
opt.model<-nlminb(start = output.matrix[which(output.matrix[,2] == min(output.matrix[,2])), 1], objective = obj.FUN)
optimized.model<-update.model(model = model, Vdose = f(opt.model$par[1]))
coeff<-summary(optimized.model)$coefficients[,4]
output.matrix<-rbind(output.matrix, c(opt.model$par[1], -logLik(optimized.model), coeff))
output.matrix<-output.matrix[order(output.matrix[,1]),]
rm(Vdose, envir = .GlobalEnv)
fit.par<-opt.model$par[1]
}
if (type == "Dvolume") { # create the approxfun and the fitting Dvolume function (FUN)
maxVol<-min(apply(X = dvh@dvh[,2:ncol(dvh@dvh)], FUN = max, MARGIN = 2)) # the minimum of max Volume
minVol<-max(apply(X = dvh@dvh[,2:ncol(dvh@dvh)], FUN = min, MARGIN = 2)) # the maximum of min Volume
# browser()
# dVol<-(maxVol - minVol)/200 # delta Volume
value <- dvh@dvh[2:(nrow(dvh@dvh) - 1), 1] # create volume vector removing extreme values
apf <- apply(X = dvh@dvh[2:(nrow(dvh@dvh) - 1), 2:ncol(dvh@dvh)], MARGIN = 2, FUN = approxfun, x = value)
update.model<-function(model, Dvolume) {
Dvolume <<- Dvolume
new.model<-update(object = model, formula. = ". ~ . + Dvolume")
return(new.model)
}
browser()
f <- function(val) return(sapply(X = apf, FUN = function(x) return(x(val))))
output.matrix<-c()
for (n in 1:length(value)) {
Dv<-f(value[n])
temp.model<-update.model(model = model, Dvolume = Dv)
coeff<-summary(temp.model)$coefficients[,4]
if ( length(coeff) == length(temp.model$coefficients) ) # check length to prevent joining models with only one covariate
output.matrix<-rbind(output.matrix, c(value[n], -logLik(temp.model), coeff)) else
warning(paste("Dvolume =", value[n], "doesn't fit a model"))
}
obj.FUN<-function(x) { # objective function for optimizing the best model
Dvolume<<-f(x[1])
new.model<-update(object = model, formula. = ". ~ . + Dvolume")
return(-logLik(new.model))
}
opt.model<-nlminb(start = output.matrix[which(output.matrix[,2] == min(output.matrix[,2])), 1], objective = obj.FUN)
optimized.model<-update.model(model = model, Dvolume = f(opt.model$par[1]))
output.matrix<-rbind(output.matrix, c(opt.model$par[1], -logLik(optimized.model), summary(optimized.model)$coefficients[,4]))
output.matrix<-output.matrix[order(output.matrix[,1]),]
rm(Dvolume, envir = .GlobalEnv)
fit.par<-opt.model$par[1]
}
names(fit.par)<-type
colnames(output.matrix)<-c("Dose", "-LL", names(summary(optimized.model)$coefficients[,4]))
# calculating confidence interval for Dvolume and Vdose
CI.val <- NULL
if (CI == TRUE) {
bound <- qchisq(CI.width,1)/2
### confidence interval for Vdose
if (type == "Vdose") {
#### start caluclation for low bound of CI ####
start.value<-fit.par # start with optimal parameter value
delta.value<-(fit.par - value[1])/2 # delta for steps going down in the CI boundary
while ( delta.value > epsilon ) {
start.value <- start.value - delta.value # create step for the value to be fitted
while (( start.value >= 0 ) && ((logLik(object = update.model(model = model, Vdose = f(start.value))) -logLik(optimized.model) + bound) > 0 ))
start.value <- start.value - delta.value
start.value <- start.value + delta.value
delta.value <- delta.value / 2
}
low.bound.CI <- start.value - delta.value / 2
#### start calculation for highbound of CI ####
start.value<-fit.par # start with optimal parameter value
delta.value<-(value[length(value)] - fit.par)/2 # delta for steps going up in the CI boundary
while ( delta.value > epsilon ) {
start.value <- start.value + delta.value # create step for the value to be fitted
while (( start.value <= value[length(value)]) && ((logLik(object = update.model(model = model, Vdose = f(start.value))) -logLik(optimized.model) + bound) > 0 ))
start.value <- start.value + delta.value
start.value <- start.value - delta.value
delta.value <- delta.value / 2
}
high.bound.CI <- start.value + delta.value / 2
}
### confidence interval for Dvolume
if (type == "Dvolume") {
#### start caluclation for low bound of CI ####
start.value<-fit.par # start with optimal parameter value
delta.value<-(fit.par - value[1])/2 # delta for steps going down in the CI boundary
while ( delta.value > epsilon ) {
start.value <- start.value - delta.value # create step for the value to be fitted
while (( start.value >= 0 ) && ((logLik(object = update.model(model = model, Dvolume = f(start.value))) -logLik(optimized.model) + bound) > 0))
start.value <- start.value - delta.value
start.value <- start.value + delta.value
delta.value <- delta.value / 2
}
low.bound.CI <- start.value - delta.value / 2
#### start calculation for highbound of CI ####
start.value<-fit.par # start with optimal parameter value
delta.value<-(value[length(value)] - fit.par)/2 # delta for steps going up in the CI boundary
while ( delta.value > epsilon ) {
start.value <- start.value + delta.value # create step for the value to be fitted
while (( start.value <= value[length(value)]) && ((logLik(object = update.model(model = model, Dvolume = f(start.value))) -logLik(optimized.model) + bound) > 0 ))
start.value <- start.value + delta.value
start.value <- start.value - delta.value
delta.value <- delta.value / 2
}
high.bound.CI <- start.value + delta.value / 2
}
CI.val <- c(low.bound.CI, high.bound.CI)
names(CI.val) <- c(paste((1-CI.width)/2, "%"), paste((1+CI.width)/2, "%"))
}
result <- new('DoseVolumeModel')
result@output.matrix <- output.matrix
result@fitted.model <- optimized.model
result@fitted.value <- fit.par
result@CI <- CI.val
result@fitted.parameter <- type
# return(list(output.matrix = output.matrix, optimized.model = optimized.model, fit.par = fit.par, CI = CI.val))
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.