From 5dcd7a254237a2a1c280b306f8cd811d47d9c7ae Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 13 Apr 2026 08:25:58 +0000 Subject: [PATCH 1/4] Initial plan From e441eaf19da319b1fdb1cbcf15356d5866130f41 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 13 Apr 2026 08:42:14 +0000 Subject: [PATCH 2/4] Fix incorrect standard errors for absolute ED calculations in model-level edfct functions When type='absolute', the edfct functions in model files set derivatives for c (lower limit) and d (upper limit) to 0. This is incorrect because the absolute-to-relative conversion makes p depend on c and d, requiring chain-rule corrections. Use central differences to compute correct c and d derivatives for absolute type in all affected model files. Fixes: braincousens.R, fplogistic.R, llogistic.R, llogistic2.R, lnormal.R, weibull1.R, weibull2.R, logistic.R Agent-Logs-Url: https://github.com/hreinwald/drc/sessions/7ea3ad01-6a69-4d8a-8c99-5b419f07b01f Co-authored-by: hreinwald <115988583+hreinwald@users.noreply.github.com> --- R/braincousens.R | 29 ++++++++++++++++++++++++++++- R/fplogistic.R | 18 ++++++++++++++++++ R/llogistic.R | 20 ++++++++++++++++++++ R/llogistic2.R | 19 +++++++++++++++++++ R/lnormal.R | 24 ++++++++++++++++++++++++ R/logistic.R | 23 ++++++++++++++++++++++- R/weibull1.R | 20 +++++++++++++++++++- R/weibull2.R | 30 +++++++++++++++++++++++++++++- 8 files changed, 179 insertions(+), 4 deletions(-) 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..4f18dc4 100644 --- a/R/logistic.R +++ b/R/logistic.R @@ -100,9 +100,10 @@ fctName, fctText) } ## Defining the ED function - edfct <- function(parm, p, ...) + edfct <- function(parm, respl, reference = "control", type = "relative", ...) { parmVec[notFixed] <- parm + p <- EDhelper(parmVec, respl, reference, type) ## 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 +128,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 (EDhelper), where p depends on c and d. + if (identical(type, "absolute")) { + .edval <- function(pv) { + p0 <- EDhelper(pv, respl, reference, type) + .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..e08d5e9 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,34 @@ 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) + 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 } From 94ec2d468c3b165878747b41f681f4b955d665f7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 13 Apr 2026 08:53:32 +0000 Subject: [PATCH 3/4] Fix weibull2.R edfct0: add EDhelper swap for b<0 to match weibull1 delegation path Agent-Logs-Url: https://github.com/hreinwald/drc/sessions/7ea3ad01-6a69-4d8a-8c99-5b419f07b01f Co-authored-by: hreinwald <115988583+hreinwald@users.noreply.github.com> --- R/weibull2.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/R/weibull2.R b/R/weibull2.R index e08d5e9..2c49e32 100644 --- a/R/weibull2.R +++ b/R/weibull2.R @@ -149,7 +149,10 @@ fctName, fctText) 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])) } From aef7e9f971f94618b3874a4a4f26a46ec48eebdc Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 13 Apr 2026 10:48:39 +0000 Subject: [PATCH 4/4] Fix logistic.R edfct: replace EDhelper with inline conversion to avoid incorrect p-swap The logistic model has opposite b-sign convention from log-logistic: b < 0 means increasing (not decreasing). EDhelper swaps p for b < 0, which is correct for log-logistic but wrong for logistic. Replace with inline absolute-to-relative conversion that skips the p-swap. Agent-Logs-Url: https://github.com/hreinwald/drc/sessions/801bd853-cb45-4ae2-8e75-c28dae361ad5 Co-authored-by: hreinwald <115988583+hreinwald@users.noreply.github.com> --- R/logistic.R | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/R/logistic.R b/R/logistic.R index 4f18dc4..6a9fe5a 100644 --- a/R/logistic.R +++ b/R/logistic.R @@ -103,7 +103,17 @@ fctName, fctText) edfct <- function(parm, respl, reference = "control", type = "relative", ...) { parmVec[notFixed] <- parm - p <- EDhelper(parmVec, respl, reference, type) + + ## 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 @@ -130,10 +140,10 @@ fctName, fctText) ## 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. + ## the absolute-to-relative conversion, where p depends on c and d. if (identical(type, "absolute")) { .edval <- function(pv) { - p0 <- EDhelper(pv, respl, reference, type) + p0 <- 100 * ((pv[3] - respl) / (pv[3] - pv[2])) .expr2 <- 100 / p0 .expr4 <- .expr2^(1 / pv[5]) .expr5 <- .expr4 - 1