diff --git a/R/braincousens.R b/R/braincousens.R index cb8e45f..2d241fd 100644 --- a/R/braincousens.R +++ b/R/braincousens.R @@ -126,7 +126,34 @@ fctName, fctText) derDose <- tempVal*tempVal1*parmVec[1]/EDdose-parmVec[5]/tempVal2 EDder <- derParm/derDose - + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (EDhelper), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper(pv, respl, reference, type) + tv0 <- (100 - p0) / 100 + helpEqn0 <- function(dose) { + ev <- exp(pv[1] * (log(dose) - log(pv[4]))) + pv[5] * (1 + ev * (1 - pv[1])) - (pv[3] - pv[2]) * ev * pv[1] / dose + } + maxAt0 <- uniroot(helpEqn0, interval)$root + eqn0 <- function(dose) { + tv0 * (1 + exp(pv[1] * (log(dose) - log(pv[4])))) - + (1 + pv[5] * dose / (pv[3] - pv[2])) + } + uniroot(eqn0, lower = maxAt0, upper = upper)$root + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + return(list(EDp, EDder[notFixed])) } diff --git a/R/fplogistic.R b/R/fplogistic.R index 3e66771..a40c325 100644 --- a/R/fplogistic.R +++ b/R/fplogistic.R @@ -152,6 +152,24 @@ fctName, fctText) denVal <- parmVec[1] * p1 * (logEDp)^(p1-1) + parmVec[4] * p2 * (logEDp)^(p2-1) derVec <- (EDp+1) * c(logEDp^p1, logEDp^p2) / denVal EDder <- c(derVec[1], 0, 0, derVec[2]) + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (EDhelper2), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper2(pv, respl, reference, type, pv[1] > 0) + invfp(log((100 - p0) / p0), pv[1], pv[4]) + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + if (loged) { EDder <- EDder / EDp diff --git a/R/llogistic.R b/R/llogistic.R index fa46d80..618e43b 100644 --- a/R/llogistic.R +++ b/R/llogistic.R @@ -170,6 +170,26 @@ fctName, fctText) EDp*c(-log(expTerm-1)/(parmVec[1]^2), 0, 0, 1/parmVec[4], expTerm*tempVal/(parmVec[5]^2)*(1/parmVec[1])*((expTerm-1)^(-1))) + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (EDhelper), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper(pv, respl, reference, type) + tv0 <- log((100 - p0) / 100) + et0 <- exp(-tv0 / pv[5]) + if (is.na(et0) || et0 <= 1) return(Inf) + pv[4] * (et0 - 1)^(1 / pv[1]) + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } } return(list(EDp, EDder[notFixed])) diff --git a/R/llogistic2.R b/R/llogistic2.R index 85b9395..d041e89 100644 --- a/R/llogistic2.R +++ b/R/llogistic2.R @@ -241,6 +241,25 @@ fctName, fctText) 0, 0, 1, tempVal1^(1/parmVec[5]-1)/(parmVec[1]*parmVec[5]*(tempVal1^(1/parmVec[5]-1)))) + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (EDhelper), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper(pv, respl, reference, type) + tv1 <- 100 / (100 - p0) + tv2 <- log(tv1^(1 / pv[5]) - 1) + pv[4] + tv2 / pv[1] + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + lEDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + return(list(lEDp, lEDder[notFixed])) } diff --git a/R/lnormal.R b/R/lnormal.R index 9d89b55..aeeaffb 100644 --- a/R/lnormal.R +++ b/R/lnormal.R @@ -205,6 +205,30 @@ fctName, fctText, loge = FALSE) } EDp <- EDfct(parmVec[1], parmVec[2], parmVec[3], parmVec[4]) EDder <- attr(EDfct(parmVec[1], parmVec[2], parmVec[3], parmVec[4]), "gradient") + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (absToRel), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- absToRel(pv, respl, type) + p0 <- 100 - p0 # reversal for absolute type + pProp0 <- 1 - (100 - p0) / 100 + if (!loge) { + pv[4] * exp(qnorm(pProp0) / pv[1]) + } else { + pv[4] + qnorm(pProp0) / pv[1] + } + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + return(list(EDp, EDder[notFixed])) } diff --git a/R/logistic.R b/R/logistic.R index d23ff3e..6a9fe5a 100644 --- a/R/logistic.R +++ b/R/logistic.R @@ -100,9 +100,20 @@ fctName, fctText) } ## Defining the ED function - edfct <- function(parm, p, ...) + edfct <- function(parm, respl, reference = "control", type = "relative", ...) { parmVec[notFixed] <- parm + + ## Convert absolute response level to relative. + ## Note: unlike log-logistic models where b < 0 means decreasing, + ## the logistic model has b < 0 = increasing. EDhelper's p-swap + ## (for b < 0, relative type) would be wrong here, so we perform + ## only the absolute-to-relative conversion inline. + if (identical(type, "absolute")) { + p <- 100 * ((parmVec[3] - respl) / (parmVec[3] - parmVec[2])) + } else { + p <- respl + } ## deriv(~e + log((100/(100-p))^(1/f) - 1) / b, c("b", "c", "d", "e", "f"), function(b,c,d,e,f){}) ## evaluated at the R prompt @@ -127,6 +138,26 @@ fctName, fctText) EDp <- as.numeric(EDcalc) EDder <- attr(EDcalc, "gradient") + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion, where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- 100 * ((pv[3] - respl) / (pv[3] - pv[2])) + .expr2 <- 100 / p0 + .expr4 <- .expr2^(1 / pv[5]) + .expr5 <- .expr4 - 1 + pv[4] + log(.expr5) / pv[1] + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + return(list(EDp, EDder[notFixed])) } diff --git a/R/weibull1.R b/R/weibull1.R index 01e336e..dbd203a 100644 --- a/R/weibull1.R +++ b/R/weibull1.R @@ -143,7 +143,25 @@ fctName, fctText) EDp <- exp(tempVal/parmVec[1] + log(parmVec[4])) EDder <- EDp*c(-tempVal/(parmVec[1]^2), 0, 0, 1/parmVec[4]) - + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The analytical derivatives above miss the chain-rule contribution from + ## the absolute-to-relative conversion (EDhelper), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper(pv, respl, reference, type) + tv0 <- log(-log((100 - p0) / 100)) + exp(tv0 / pv[1] + log(pv[4])) + } + .eps <- .Machine$double.eps + for (.i in c(2, 3)) { + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + EDder[.i] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + return(list(EDp, EDder[notFixed])) } diff --git a/R/weibull2.R b/R/weibull2.R index 00f37d1..2c49e32 100644 --- a/R/weibull2.R +++ b/R/weibull2.R @@ -130,6 +130,7 @@ fctName, fctText) edfct <- function(parm, p, reference, type, ...) { parmVec[notFixed] <- parm + respl <- p # save original response level p <- absToRel(parmVec, p, type) @@ -139,7 +140,37 @@ fctName, fctText) p <- 100 - p } - weibull1(fixed, names)$edfct(parm, p, reference, "relative", ...) + result <- weibull1(fixed, names)$edfct(parm, p, reference, "relative", ...) + + ## Fix: correct c and d derivatives for absolute type using central differences. + ## The delegation to weibull1 with type="relative" produces zero derivatives + ## for c and d, missing the chain-rule contribution from the + ## absolute-to-relative conversion (absToRel) where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- absToRel(pv, respl, type) + # Replicate weibull2's reversal (for b > 0 and absolute type) + if (pv[1] > 0 && identical(reference, "control")) p0 <- 100 - p0 + # Replicate weibull1's EDhelper swap (for b < 0 and relative type) + if (pv[1] < 0 && identical(reference, "control")) p0 <- 100 - p0 + tv0 <- log(-log((100 - p0) / 100)) + exp(tv0 / pv[1] + log(pv[4])) + } + .eps <- .Machine$double.eps + .nfIdx <- which(notFixed) + for (.i in c(2, 3)) { + if (!notFixed[.i]) next + .h <- if (abs(parmVec[.i]) > sqrt(.eps)) abs(parmVec[.i]) * .eps^(1/3) else .eps^(1/3) + .pvUp <- replace(parmVec, .i, parmVec[.i] + .h) + .pvDn <- replace(parmVec, .i, parmVec[.i] - .h) + .pos <- which(.nfIdx == .i) + if (length(.pos) == 1L) { + result[[2]][.pos] <- (.edval(.pvUp) - .edval(.pvDn)) / (2 * .h) + } + } + } + + result }