Fix model-level edfct derivatives for absolute ED type

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]))
     }
 

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

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]))

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]))
     }
 

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]))
     }
 

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]))
     }
 

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]))
     }
 

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
     }