Nothing
#' Estimation of the equivalent dose (ED) value for the SAR protocol
#'
#' Internal function called by \link{analyse_TL.SAR}. \cr
#' This function estimates the equivalent dose (ED) based on a doses vector and a Lx/Tx vector provided. \cr
#' See details for more information.
#'
#' @param doses
#' \link{numeric} (\bold{required}): doses vector
#' @param LnTn
#' \link{numeric} (\bold{required}): Ln/Tn.
#' @param LnTn.error
#' \link{numeric} (\bold{required}): Error for the Ln/Tn.
#' @param LxTx
#' \link{numeric} (\bold{required}): Lx/Tx vector
#' @param LxTx.error
#' \link{numeric} (\bold{required}): Error for the Lx/Tx vector
#' @param fitting.parameters
#' \link{list} (with default): fitting parameters. See details.
#'
#' @details
#' This function estimates the equivalent dose based on the doses vector, Ln/Tn and the Lx/Tx matrix provided. \cr
#' Different fitting methods are available (\code{LIN}, \code{EXP}, \code{EXP+LIN} or \code{EXP+EXP}).
#' Moreover, the fitting can be weigthed or not. \cr
#'
#' \bold{Fitting parameters} \cr
#' The fitting parameters are: \cr
#' \describe{
#' \item{\code{method}}{
#' \link{character}: Fitting method (\code{LIN}, \code{EXP}, \code{EXP+LIN} or \code{EXP+EXP}).}
#' \item{\code{fit.weighted}}{
#' \link{logical}: If the fitting is weighted or not.}
#' \item{\code{fit.use.slope}}{
#' \link{logical}: If the slope of the Q growth curve is reused for the sublinearity correction.}
#' \item{\code{fit.rDoses.min}}{
#' \link{numeric}: Lowest regenerative dose used for the fitting.}
#' \item{\code{fit.rDoses.max}}{
#' \link{numeric}: Highest regenerative dose used for the fitting.}
#' }
#'
#' @return
#' The function provides an \linkS4class{TLum.Results} object containing: \cr
#' \describe{
#' \item{\code{GC}}{
#' \link{list}: fitting curve.}
#' \item{\code{Q}}{
#' \link{numeric}: equivalent dose estimation}
#' \item{\code{Q.error}}{
#' \link{numeric}: Error for the equivalent dose estimation}
#' \item{\code{summary}}{
#' \link{list}: parameters of the fitting result.}
#' }
#'
#' @details
#' \bold{Warning}: This function is an internal function and should not be used except for development purposes.
#' Internal functions can be heavily modified and even renamed or removed in new version of the package.
#'
#' @seealso
#' \link{analyse_TL.SAR}.
#'
#' @author David Strebler, University of Cologne (Germany).
#'
#' @export calc_TL.SAR.fit
calc_TL.SAR.fit <- function(
doses,
LnTn,
LnTn.error,
LxTx,
LxTx.error,
fitting.parameters=list(fit.method="LIN",
fit.weighted=FALSE)
){
fit.method <- fitting.parameters$fit.method
fit.weighted <- fitting.parameters$fit.weighted
# Integrity Check ---------------------------------------------------------
if (missing(LnTn)){
stop("[calc_TL.SARfit] Error: LnTn is missing.")
}
if (missing(LnTn.error)){
stop("[calc_TL.SARfit] Error: LnTn.error is missing.")
}
if (missing(LxTx)){
stop("[calc_TL.SARfit] Error: LxTx is missing.")
}
if (missing(LxTx.error)){
stop("[calc_TL.SARfit] Error: LxTx.error is missing.")
}
if (missing(doses)){
stop("[calc_TL.SARfit] Error: doses is missing.")
}
# ------------------------------------------------------------------------------
w <- 1/(LxTx.error^2)
w[!is.finite(w)] <- 0.0001
LxTx[!is.finite(LxTx)] <- 0.0001
data <- data.frame(x=doses,
y=LxTx)
#LIN
function.LIN <- function(a,b,x){a+b*x}
#EXP
function.EXP <- function(a,b,c,x){a*(1-exp(-(x+c)/b))}
#EXP+LIN
function.EXPLIN <- function(a,b,c,g,x){a*(1-exp(-(x+c)/b)+(g*x))}
#EXP+EXP
function.EXPEXP <- function(a1,a2,b1,b2,x){(a1*(1-exp(-(x)/b1)))+(a2*(1-exp(-(x)/b2)))}
if(fit.method == "LIN"){
# Weighted
if(fit.weighted){
fit <- lm(formula = y ~ x,
data = data,
weights=w)
# Not weighted
}else{
fit <- lm(formula = y ~ x,
data = data)
}
a <- summary(fit,correlation = TRUE)$coefficients[1,"Estimate"]
a.error <- summary(fit,correlation = TRUE)$coefficients[1,"Std. Error"]
b <- summary(fit,correlation = TRUE)$coefficients[2,"Estimate"]
b.error <- summary(fit,correlation = TRUE)$coefficients[2,"Std. Error"]
r <- abs(summary(fit,correlation = TRUE)$correlation[1,2])
s <- list(a = a,
a.error = a.error,
b = b,
b.error = b.error,
r = r)
Q <- (LnTn - a)/b
#Q.error <- sqrt((a.error/a)^2 + (b.error/b)^2 + 2*(a.error/a)*(b.error/b)*r ) * Q
#Q.error <- sqrt((a.error/a)^2 + (b.error/b)^2 ) * Q
a2 <- abs(LnTn - a)
a2.error <- sqrt(sum(LnTn.error^2,a.error^2,na.rm = TRUE))
Q.error <- sqrt(sum((a2.error/a2)^2,(b.error/b)^2,abs(2*(a2.error/a2)*(b.error/b)*r), na.rm = TRUE))*Q
}else if(fit.method == "EXP"){
# initial parameters
# a
a.i<-max(LxTx)
# b
temp.lm <- try(lm(log(data$y)~data$x),
silent = TRUE)
if(class(temp.lm) =="try-error"){
b.i<-1
}else{
b.i <- as.numeric(1/temp.lm$coefficients[2])
}
#c
temp.lm <- try(lm(data$y~data$x),
silent = TRUE)
if(class(temp.lm) =="try-error"){
c.i<-0
}else{
c.i <- as.numeric(abs(temp.lm$coefficients[1]/temp.lm$coefficients[2]))
}
#Fitting
# Weighted
if(fit.weighted){
fit<-try(nls(y~function.EXP(a,b,c,x),
data=data,
start=c(a=a.i,
b=b.i,
c=c.i),
weights=w,
trace=FALSE,
algorithm="port",
nls.control(maxiter=500)),
silent=TRUE)
# Not weighted
}else{
fit<-try(nls(y~function.EXP(a,b,c,x),
data=data,
start=c(a=a.i,
b=b.i,
c=c.i),
trace=FALSE,
algorithm="port",
nls.control(maxiter=500)),
silent=TRUE)
}
if (class(fit)=="try-error"){
print("pas ok")
fit <- NA
s <- list(a = NA,
a.error = NA,
b = NA,
b.error = NA,
c = NA,
c.error = NA)
Q <- 0
Q.error <- NA
}else{
print("ok")
res <- summary(fit)$coefficients[,"Estimate"]
res.error.r <- summary(fit)$coefficients[,"Std. Error"]
res.error <- abs(res*res.error)
s <- list(a = res["a"],
a.error = res.error["a"],
b = res["b"],
b.error = res.error["b"],
c = res["c"],
c.error = res.error["c"])
Q <- -res["c"]-res["b"]*log(1-LnTn/res["b"])
Q.error <- 1 # temp
}
}else{
stop("[calc_TL.SARfit] Error: fit.method not supported.")
}
result <- list(GC=fit,
Q=Q,
Q.error=Q.error,
summary=s)
new.originator <- as.character(match.call()[[1]])
new.plotData <- list()
new.TLum.Results.calc_TL.SAR.fit <- set_TLum.Results(originator= new.originator,
data = result,
plotData = new.plotData)
return (new.TLum.Results.calc_TL.SAR.fit)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.