"llogistic2" <- function(
fixed = c(NA, NA, NA, NA, NA), names = c("b", "c", "d", "e", "f"), ss = c("1", "2", "3"), ssfct = NULL,
fctName, fctText)
{
## Matching 'adjust' argument
ss <- match.arg(ss)
## Checking arguments
numParm <- 5
if (!is.character(names) | !(length(names) == numParm)) {stop("Not correct 'names' argument")}
if ( !(length(fixed) == numParm) ) {stop("Not correct 'fixed' argument")}
## Handling 'fixed' argument
notFixed <- is.na(fixed)
parmVec <- rep(0, numParm)
parmVec[!notFixed] <- fixed[!notFixed]
## Defining the basic non-linear function
bfct <- function(x, parm)
{
parm[2] + (parm[3]-parm[2])/((1+(x/exp(parm[4]))^parm[1]))^parm[5]
}
## Defining the non-linear function
fct <- function(dose, parm)
{
parmMat <- matrix(parmVec, nrow(parm), numParm, byrow=TRUE)
parmMat[, notFixed] <- parm
cParm <- parmMat[, 2]
cParm + (parmMat[, 3] - cParm)/((1+exp(parmMat[, 1]*(log(dose)-parmMat[, 4])))^parmMat[, 5])
}
## Defining the self starter function
## Defining self starter function based on supplied function
if (!is.null(ssfct))
{
ssfct <- ssfct
} else {
## Version 1 (default)
if (ss == "1")
{
ssfct <- function(dframe)
{
x <- dframe[, 1]
y <- dframe[, 2]
startVal <- rep(0, numParm)
startVal[3] <- max(y) + 0.001 # the d parameter
startVal[2] <- min(y) - 0.001 # the c parameter
startVal[5] <- 1 # better choice may be possible!
if (length(unique(x))==1) {return((c(NA, NA, startVal[3], NA, NA))[notFixed])}
# only estimate of upper limit if a single unique dose value
indexT2 <- (x > 0 & is.finite(x))
# if (!any(indexT2)) {return((rep(NA, numParm))[notFixed])} # for negative dose value
x2 <- x[indexT2]
y2 <- y[indexT2]
startVal[c(1,4)] <- find.be2(x2, y2, startVal[2] - 0.001, startVal[3])
startVal[4] <- log(startVal[4])
# 0.001 to avoid 0 in the denominator
return(startVal[notFixed])
}
}
## Version 2
if (ss == "2")
{
ssfct <- function(dframe)
{
x <- dframe[, 1]
y <- dframe[, 2]
dVal <- ifelse(notFixed[3], 1.01*max(y), fixed[3])
cVal <- ifelse(notFixed[2], 0.99*min(y), fixed[2])
fVal <- 1 # need not be updated with value in 'fixed[5]'
# better choice than 1 may be possible!
if ( length(unique(x)) == 1 ) {return((c(NA, NA, dVal, NA, NA))[notFixed])}
# only estimate of upper limit if a single unique dose value
# Cutting away response values close to d
indexT1a <- x > 0
x2 <- x[indexT1a]
y2 <- y[indexT1a]
beVec <- find.be2(x2, y2, cVal, dVal)
bVal <- beVec[1]
eVal <- log(beVec[2])
return(as.vector(c(bVal, cVal, dVal, eVal, fVal)[notFixed]))
}
}
## Version 3
if (ss == "3")
{
ssfct <- function(dframe)
{
x <- dframe[, 1]
y <- dframe[, 2]
cVal <- ifelse(notFixed[2], 0.99 * min(y), fixed[2])
dVal <- ifelse(notFixed[3], 1.01 * max(y), fixed[3])
fVal <- 1 # need not be updated with value in 'fixed[5]'
if ( length(unique(x)) == 1 ) {return((c(NA, NA, dVal, NA, NA))[notFixed])}
# only estimate of upper limit if a single unique dose value
beVec <- find.be1(x, y, cVal, dVal)
bVal <- beVec[1]
eVal <- log(beVec[2])
return(as.vector(c(bVal, cVal, dVal, eVal, fVal)[notFixed]))
}
}
}
## Finding b and e based on stepwise increments
find.be1 <- function(x, y, c, d)
{
unix <- unique(x)
uniy <- tapply(y, x, mean)
lenx <- length(unix)
j <- 2
for (i in 2:lenx)
{
crit1 <- (uniy[i] > (d + c)/2) && (uniy[i-1] < (d + c)/2)
crit2 <- (uniy[i] < (d + c)/2) && (uniy[i-1] > (d + c)/2)
if (crit1 || crit2) break
j <- j + 1
}
eVal <- (unix[j] + unix[j-1])/2
bVal <- -sign(uniy[j] - uniy[j-1]) # -(uniy[j] - uniy[j-1]) / (unix[j] - unix[j-1])
return(as.vector(c(bVal, eVal)))
}
## Finding b and e based on linear regression
find.be2 <- function(x, y, c, d)
{
logitTrans <- log((d - y)/(y - c))
lmFit <- lm(logitTrans ~ log(x))
coefVec <- coef(lmFit)
bVal <- coefVec[2]
eVal <- exp(-coefVec[1]/bVal)
return(as.vector(c(bVal, eVal)))
}
## Defining names
names <- names[notFixed]
##Defining the first derivatives (in the parameters)
deriv1 <- function(dose, parm)
{
parmMat <- matrix(parmVec, nrow(parm), numParm, byrow = TRUE)
parmMat[, notFixed] <- parm
t1 <- parmMat[, 3] - parmMat[, 2]
t2 <- exp(parmMat[, 1]*(log(dose) - parmMat[, 4]))
t5 <- (1 + t2)^parmMat[, 5]
cbind( -t1 * xlogx(dose/exp(parmMat[, 4]), parmMat[, 1], parmMat[, 5] + 1) * parmMat[, 5],
1 - 1/t5,
1/t5,
t1 * parmMat[, 5] * divAtInf(t2, (1 + t2)^(parmMat[, 5] + 1)) * parmMat[, 1],
-t1 * divAtInf(log(1+t2), t5) )[, notFixed]
}
deriv2 <- NULL
##Defining the first derivative (in the dose)
derivx <- function(x, parm)
{
parmMat <- matrix(parmVec, nrow(parm), numParm, byrow = TRUE)
parmMat[, notFixed] <- parm
temp1 <- exp(parmMat[, 1]*(log(x) - parmMat[, 4])) # x/parmMat[, 4]
temp2 <- 1 + temp1
temp3 <- parmMat[, 5]*(temp2^(parmMat[, 5] - 1))*temp1*parmMat[, 1]/x
temp4 <- temp2^(2*parmMat[, 5])
(-(parmMat[, 3] - parmMat[, 2])*temp3)/temp4
}
# ## Setting the limits
# if (length(lowerc) == numParm) {lowerLimits <- lowerc[notFixed]} else {lowerLimits <- lowerc}
# if (length(upperc) == numParm) {upperLimits <- upperc[notFixed]} else {upperLimits <- upperc}
## Defining the ED function
## (returning ED values and corresponding standard errors on log scale)
edfct <- function(parm, respl, reference, type, ...)
{
parmVec[notFixed] <- parm
# if (type == "absolute")
# {
# p <- 100*((parmVec[3] - respl)/(parmVec[3] - parmVec[2]))
# } else {
# p <- respl
# }
# if ( (parmVec[1] < 0) && (reference == "control") )
# {
# p <- 100 - p
# }
p <- EDhelper(parmVec, respl, reference, type)
tempVal1 <- 100/(100-p)
tempVal2 <- log(tempVal1^(1/parmVec[5]) - 1)
lEDp <- parmVec[4] + tempVal2 / parmVec[1]
### updated llostic2.R function in package
if (type=="relative"){
lEDder <-
c(-tempVal2/(parmVec[1]^2),
0, 0, 1,
tempVal1^(1/parmVec[5]-1)/(parmVec[1]*parmVec[5]*(tempVal1^(1/parmVec[5]-1))))
}else{
## ask Yan for her derivative derivation--change 0,0 part
## note below derivative w.r.t c and d are only for f==1 case!!!
lEDder <-
c(-tempVal2/(parmVec[1]^2),
tempVal1/(parmVec[1]*(parmVec[3]-parmVec[2])),
100/(parmVec[1]*(parmVec[3]-parmVec[2])*p),
1,
tempVal1^(1/parmVec[5]-1)/(parmVec[1]*parmVec[5]*(tempVal1^(1/parmVec[5]-1))))
}
return(list(lEDp, lEDder[notFixed]))
}
## Defining the inverse function
invfct <- function(y, parm)
{
parmVec[notFixed] <- parm
exp(log(((parmVec[3] - parmVec[2])/(y - parmVec[2]))^(1/parmVec[5]) - 1)/parmVec[1] + parmVec[4])
}
## Defining functions returning lower and upper limit and monotonicity
lowerAs <- pickParm(parmVec, notFixed, 2)
upperAs <- pickParm(parmVec, notFixed, 3)
monoton <- monoParm(parmVec, notFixed, 1, -1)
## Setting function details
if (missing(fctName))
{
fctName <- as.character(match.call()[[1]])
}
if (missing(fctText))
{
fctText <- "Log-logistic (log(ED50) as parameter)"
}
## Returning the function with self starter and names
returnList <-
list(fct = fct, ssfct = ssfct, names = names, deriv1 = deriv1, deriv2 = deriv2, derivx = derivx,
edfct = edfct, bfct = bfct, inversion = invfct,
# lowerc=lowerLimits, upperc=upperLimits,
name = fctName,
text = fctText,
noParm = sum(is.na(fixed)),
lowerAs = lowerAs,
upperAs = upperAs,
monoton = monoton)
class(returnList) <- "llogistic"
invisible(returnList)
}
lowupFixed <- function(modelStr, upper)
{
paste(modelStr, "with lower limit at 0 and upper limit at", upper)
}
lowFixed <- function(modelStr)
{
paste(modelStr, "with lower limit at 0")
}
upFixed <- function(modelStr, upper)
{
paste(modelStr, "with upper limit at", upper)
}
"LL2.2" <-
function(upper = 1, fixed = c(NA, NA), names = c("b", "e"), ...)
{
## Checking arguments
numParm <- 2
if (!is.character(names) | !(length(names)==numParm)) {stop("Not correct 'names' argument")}
if (!(length(fixed)==numParm)) {stop("Not correct length of 'fixed' argument")}
return( llogistic2(fixed = c(fixed[1], 0, upper, fixed[2], 1),
names = c(names[1], "c", "d", names[2], "f"),
fctName = as.character(match.call()[[1]]),
fctText = lowupFixed("Log-logistic (log(ED50) as parameter)", upper),
...) )
}
"LL2.3" <-
function(fixed = c(NA, NA, NA), names = c("b", "d", "e"), ...)
{
## Checking arguments
numParm <- 3
if (!is.character(names) | !(length(names)==numParm)) {stop("Not correct 'names' argument")}
if (!(length(fixed)==numParm)) {stop("Not correct length of 'fixed' argument")}
return( llogistic2(fixed = c(fixed[1], 0, fixed[2:3], 1),
names = c(names[1], "c", names[2:3], "f"),
fctName = as.character(match.call()[[1]]),
fctText = lowFixed("Log-logistic (log(ED50) as parameter)"),
...) )
}
"LL2.3u" <-
function(upper = 1, fixed = c(NA, NA, NA), names = c("b", "c", "e"), ...)
{
## Checking arguments
numParm <- 3
if (!is.character(names) | !(length(names)==numParm)) {stop("Not correct 'names' argument")}
if (!(length(fixed)==numParm)) {stop("Not correct length of 'fixed' argument")}
return( llogistic2(fixed = c(fixed[1:2], upper, fixed[3], 1),
names = c(names[1:2], "d", names[3], "f"),
fctName = as.character(match.call()[[1]]),
fctText = upFixed("Log-logistic (log(ED50) as parameter)", upper),
...) )
}
"LL2.4" <-
function(fixed = c(NA, NA, NA, NA), names = c("b", "c", "d", "e"), ...)
{
## Checking arguments
numParm <- 4
if (!is.character(names) | !(length(names)==numParm)) {stop("Not correct names argument")}
if (!(length(fixed)==numParm)) {stop("Not correct length of 'fixed' argument")}
return( llogistic2(fixed = c(fixed, 1), names = c(names, "f"),
fctName = as.character(match.call()[[1]]), ...) )
}
"LL2.5" <-
function(fixed = c(NA, NA, NA, NA, NA), names = c("b", "c", "d", "e", "f"), ...)
{
return( llogistic2(fixed = fixed, names = names,
fctName = as.character(match.call()[[1]]),
fctText = "Generalised log-logistic (log(ED50) as parameter)", ...) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.