Version 2.7.3

.Rbuildignore

@@ -3,4 +3,5 @@
 cache/
 README.Rmd
 README.md
-README_files/
+README_files/
+^\\.github$

.github/workflows/r.yml

@@ -1,3 +1,11 @@
+# This workflow uses actions that are not certified by GitHub.
+# They are provided by a third-party and are governed by
+# separate terms of service, privacy policy, and support
+# documentation.
+#
+# See https://github.com/r-lib/actions/tree/master/examples#readme for
+# additional example workflows available for the R community.
+
 on:
   push:
     branches: [main, master]
@@ -43,7 +51,6 @@ jobs:
           key: ${{ runner.os }}-r-${{ matrix.config.r }}-${{ hashFiles('DESCRIPTION') }}
           restore-keys: |
             ${{ runner.os }}-r-${{ matrix.config.r }}-
-
       - name: Install dependencies
         run: |
           install.packages(c("sandwich", "CVXR", "multcomp", "gridExtra", "isotone", 
@@ -68,4 +75,4 @@ jobs:
         uses: actions/upload-artifact@v4
         with:
           name: ${{ runner.os }}-r${{ matrix.config.r }}-results
-          path: check
+          path: check

DESCRIPTION

@@ -1,33 +1,38 @@
 Package: bmd
 Type: Package
 Title: Benchmark dose estimation for dose-response data
-Version: 2.7.1
+Version: 2.7.3
 Date: 2025-03-24
 Author: Signe M.Jensen, Christian Ritz and Jens Riis Baalkilde
 Maintainer: Signe M. Jensen <smj@plen.ku.dk>
-Description: Benchmark dose analysis for continuous, quantal, count and ordinal dose-response data
-Imports: 
-    drc,
+Description: Provides tools for benchmark dose analysis of dose-response data. Supports continuous, quantal, count, and ordinal responses.
+Imports:
+    drc, 
     ggplot2,
-    dplyr
+    dplyr,
+    stats
 Suggests: 
     sandwich,
     CVXR,
     multcomp,
     gridExtra,
     isotone,
     reshape2,
-    dplyr,
     car,
     Matrix,
-    RLRsim,
     scales,
+    metafor,
+    tidyr,
+    numDeriv,
     drcData,
     testthat (>= 3.0.0)
 Remotes:
-    doseresponse/drcData,
-    doseresponse/drc
-License: GPL
+    doseresponse/drc,
+    doseresponse/drcData
+License: file LICENSE
 Encoding: UTF-8
 LazyData: true
 Config/testthat/edition: 3
+RoxygenNote: 7.3.2
+Depends: 
+    R (>= 3.5)

NAMESPACE

@@ -1,15 +1,29 @@
 import(drc, ggplot2, dplyr)
+importFrom("graphics", "lines")
+importFrom("stats", "aggregate", "approx", "as.formula", "coef",
+           "complete.cases", "confint", "constrOptim", "df.residual",
+           "dnorm", "fitted", "lm", "model.frame", "model.matrix",
+           "optim", "pnorm", "predict", "qchisq", "qnorm", "qt",
+           "quantile", "rbinom", "resid", "residuals", "rnorm", "sd",
+           "uniroot", "update", "var", "vcov", "AIC", "BIC", "logLik")
+importFrom("utils", "setTxtProgressBar", "txtProgressBar")
 export(bmd, bmdBoot, bmdIso, bmdIsoBoot, PAV, bmdMA, bootDataGen, bmdMACurve, BCa, invBmd, expandBinomial, 
        getStackingWeights, drmOrdinal, bmdOrdinal, bmdOrdinalMA,
-       expandOrdinal, bootDataGenOrdinal,
-       qplotDrc, qplotBmd, MACurve, monotonicityTest, trendTest, bmdHetVar, drmHetVar)
+       qplotDrc, qplotBmd, MACurve, monotonicityTest, trendTest, bmdHetVar, drmHetVar, drmMMRE, bmdHetVarMA)
 
 ## S3 methods
 S3method(logLik, drcOrdinal)
+S3method(logLik, drcHetVar)
 S3method(AIC, drcOrdinal)
+S3method(AIC, drcHetVar)
 S3method(BIC, drcOrdinal)
+S3method(BIC, drcHetVar)
 S3method(plot, drcOrdinal)
 S3method(print, drcOrdinal)
 S3method(print, bmdOrdinal)
+S3method(print, drcHetVar)
 S3method(plot, bmd)
-S3method(plot, drcHetVar)
+S3method(plot, drcHetVar)
+S3method(vcov, drcMMRE)
+S3method(print, drcMMRE)
+S3method(coef, drcHetVar)

R/AIC.drcOrdinal.R

@@ -1,4 +1,11 @@
-AIC.drcOrdinal <- function(object, epsilon = 10^(-16)){
+AIC.drcOrdinal <- function(object, ..., k = 2) {
+  dots <- list(...)
+  if (!is.null(dots$epsilon)){
+    epsilon <- dots$epsilon
+  } else {
+    epsilon <- 1e-16
+  }
+
   n.parameters <- sum(sapply(object$drmList, function(mod) length(mod$coefficients)))
-  - 2 * logLik(object, epsilon) + 2 * n.parameters
-}
+  -2 * logLik(object, epsilon = epsilon) + k * n.parameters
+}

R/BIC.drcOrdinal.R

@@ -1,4 +1,11 @@
-BIC.drcOrdinal <- function(object, epsilon = 10^(-16)){
+BIC.drcOrdinal <- function(object, ...){
+  dots <- list(...)
+  if (!is.null(dots$epsilon)){
+    epsilon <- dots$epsilon
+  } else {
+    epsilon <- 1e-16
+  }
+
   n.parameters <- sum(sapply(object$drmList, function(mod) length(mod$coefficients)))
   n.obs <- sum(object$data[[object$weights]])
   n.parameters * log(n.obs) - 2 * logLik(object, epsilon) 

R/PAV.R

@@ -1,4 +1,5 @@
-PAV<-function(object,data,type){
+PAV<-function(formula,data,type){
+  object <- formula
   if( identical(type,"binomial")){
     N <- length( data[, paste(object[[3]]) ])
     Events <- data[, strsplit(as.character(object[[2]]),"/")[[2]] ]

R/bmd.R

@@ -17,7 +17,11 @@ bmd<-function(object, bmr, backgType = c("modelBased", "absolute", "hybridSD", "
     backgType <- "modelBased"
     } 
   if (missing(backgType)) {
-    stop(paste("backgType is missing", sep=""))
+    if(!(def %in% c("hybridExc", "hybridAdd"))){
+      backgType <- "modelBased"
+    } else {
+      backgType <- "hybridSD"
+    }
   }
   if (!(def %in% c("excess", "additional", "relative", "extra", "added", "hybridExc", "hybridAdd", "point"))) {
     stop(paste("Could not recognize def", sep=""))
@@ -29,16 +33,19 @@ bmd<-function(object, bmr, backgType = c("modelBased", "absolute", "hybridSD", "
   level <- 1-2*(1-level)
 
   interval <- match.arg(interval)
+  if(inherits(object, "drcMMRE") & interval != "delta"){
+    stop("only delta type confidence interval supported for object of type \"drcMMRE\"")
+  }
   if(interval == "sandwich"){
     sandwich.vcov <- TRUE
     interval <- "delta"
   }
-  if(sandwich.vcov & !require("sandwich")){
+  if(sandwich.vcov & !requireNamespace("sandwich")){
     stop('package "sandwich" must be installed to compute sandwich confidence intervals')
   }
   respTrans <- match.arg(respTrans)
 
-  if(class(object$fct) == "braincousens" & is.null(object$fct$fixed)){
+  if(inherits(object$fct, "braincousens") & is.null(object$fct$fixed)){
     if(object$fct$name == "BC.4"){
       object$fct$fixed <- c(NA, 0, NA, NA, NA)
     } else if(object$fct$name == "BC.5"){
@@ -82,8 +89,12 @@ bmd<-function(object, bmr, backgType = c("modelBased", "absolute", "hybridSD", "
       }
       dBmdVal <- EDeval[[2]] + bmrScaledList$dBmrScaled[,1] / fctDerivx(bmdVal, t(parmMat))[1]
       bmdSEVal <- sqrt(dBmdVal %*% varCov %*% dBmdVal)
-      intMat <- drc:::confint.basic(matrix(c(bmdVal, bmdSEVal), ncol = 2), 
-                                    level = level, object$"type", df.residual(object), FALSE)
+      if(inherits(object, "drcMMRE")){
+        intMat <- matrix(qnorm(c((1-level)/2, 1-(1-level)/2), mean = bmdVal, sd = bmdSEVal), ncol = 2)
+      } else {
+        intMat <- drc:::confint.basic(matrix(c(bmdVal, bmdSEVal), ncol = 2), 
+                                      level = level, object$"type", df.residual(object), FALSE)
+      }
     } else if(interval == "inv"){
       if(!identical(respTrans, "none")){stop("inverse regression interval not available for transformed response.")}
       slope <- drop(ifelse(object$curve[[1]](0)-object$curve[[1]](Inf)>0,"decreasing","increasing"))

R/bmdBoot.R

@@ -38,7 +38,7 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
 
   respTrans <- match.arg(respTrans)
 
-  if(class(object$fct) == "braincousens" & is.null(object$fct$fixed)){
+  if(inherits(object$fct, "braincousens") & is.null(object$fct$fixed)){
     if(object$fct$name == "BC.4"){
       object$fct$fixed <- c(NA, 0, NA, NA, NA)
     } else if(object$fct$name == "BC.5"){
@@ -57,16 +57,17 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
   get.drm.list <- function(tmp.data){
     if(ncol(object$parmMat) == 1){
       drm.list.tmp <- lapply(tmp.data, function(x){
-        try(eval(substitute(drm(formula0, data = x, type = object$type, fct = object[["fct"]]),
+        try(eval(substitute(drm(formula0, data = x, type = object$type, fct = object[["fct"]],
+                                control = drmc(noMessage = TRUE)),
                             list(formula0 = object$call$formula)
                             )), TRUE)
       }
       )
     } else if(is.null(object$call$pmodels)){
       drm.list.tmp <- lapply(tmp.data, function(x){
-        if(object$type != "binomial"){
-          x[[as.character(object$call$curveid)]] <- x[[paste0("orig.", as.character(object$call$curveid))]]
-        }
+        # if(object$type != "binomial"){
+        #   x[[as.character(object$call$curveid)]] <- x[[paste0("orig.", as.character(object$call$curveid))]]
+        # }
         try(
           eval(substitute(drm(object$call$formula, weights = weights0, curveid = curveid0,
                               data = x, type = object$type, fct = object$fct, control = drmc(noMessage = TRUE)),
@@ -77,11 +78,11 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
       })
     } else {
       drm.list.tmp <- lapply(tmp.data, function(x){
-        if(object$type != "binomial"){
-          x[[as.character(object$call$curveid)]] <- x[[paste0("orig.", as.character(object$call$curveid))]]
-        }
+        # if(object$type != "binomial"){
+        #   x[[as.character(object$call$curveid)]] <- x[[paste0("orig.", as.character(object$call$curveid))]]
+        # }
         try(
-          eval(substitute(drm(formula0, weights = weights0, curveid = curveid0,pmodels = pmodels0,
+          eval(substitute(drm(formula0, weights = weights0, curveid = curveid0, pmodels = pmodels0,
                               data = x, type = object$type, fct = object$fct, control = drmc(noMessage = TRUE)),
                           list(formula0 = object$call$formula,
                                weights0 = object$call$weights,
@@ -138,7 +139,7 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
     if(identical(bootInterval,"BCa")){
       jackData <- list()
       for(i in 1:(dim(object$data)[1])){
-        jackData[[i]] <- object$data[-i,]
+        jackData[[i]] <- object$origData[-i,]
       }
       # bootJack.drm.tmp <- lapply(jackData, function(x){
       #   try(drm(object$call$formula, data = x, fct = object[["fct"]]),TRUE)

R/bmdHetVar.R

@@ -1,12 +1,9 @@
-bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybridPercentile"), backg = NA, def = c("hybridExc", "hybridAdd"), interval = c("boot", "none"), R = 1000, level = 0.95, progressInfo = TRUE, display = TRUE){
+bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybridPercentile"), backg = NA, def = c("hybridExc", "hybridAdd"), interval = c("boot", "none"), R = 1000, level = 0.95, bootType = "nonparametric", progressInfo = TRUE, display = TRUE){
   ### Assertions ###
   # object
-  if(!identical(class(object),"drcHetVar")){
+  if(!inherits(object,"drcHetVar")){
     stop('object must be a dose-response model with a heterogeneous variance structure of class "drcHetVar" ')
   }
-  if(length(unique(object$model$dataList$curveid)) != 1){
-    stop("dose-response models with multiple curves not supported for heteroscedasticity analysis")
-  }
 
   # bmr
   if(missing(bmr)){
@@ -35,12 +32,21 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     stop('Could not recognize def. Options are "hybridExc" or "hybridAdd"')
   }
 
+  if(length(bootType) != 1){
+    if(!identical(bootType, c("nonparametric", "semiparametric", "parametric"))){
+    }
+    bootType <- "nonparametric" # default
+  }
+  if(!bootType %in% c("nonparametric", "semiparametric", "parametric")){
+    stop('"bootType" not recognised. Options are: "nonparametric", "semiparametric" and "parametric"')
+  }
+
   level <- 1-2*(1-level)
 
   # SLOPE
-  slope <- drop(ifelse(object$model$curve[[1]](0)-object$model$curve[[1]](Inf)>0,"decreasing","increasing"))
-  if(is.na(object$model$curve[[1]](0)-object$model$curve[[1]](Inf))){
-    slope <- drop(ifelse(object$model$curve[[1]](0.00000001)-object$model$curve[[1]](100000000)>0,"decreasing","increasing"))
+  slope <- drop(ifelse(object$curve(0)-object$curve(Inf)>0,"decreasing","increasing"))
+  if(is.na(object$curve(0)-object$curve(Inf))){
+    slope <- drop(ifelse(object$curve(0.00000001)-object$curve(100000000)>0,"decreasing","increasing"))
   }
 
   # sigmaFun
@@ -53,7 +59,7 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
       if(is.na(backg)){
         stop('backgType = absolute, but backg not supplied')
       }
-      p0 <- 1 - pnorm((backg - object$model$curve[[1]](0)) / sigmaFun0(0))
+      p0 <- 1 - pnorm((backg - object$curve(0)) / sigmaFun0(0))
     }
     if(identical(backgType, "hybridPercentile")) {
       p0 <- ifelse(is.na(backg),1-0.9,1-backg)
@@ -66,17 +72,17 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     bmrScaled <- switch(
       def,
       hybridExc = function(x){ sigmaFun0(x) * 
-          (qnorm(1 - p0) - qnorm(1 - p0 - (1 - p0)*bmr)) + object$model$curve[[1]](0)},
+          (qnorm(1 - p0) - qnorm(1 - p0 - (1 - p0)*bmr)) + object$curve(0)},
       hybridAdd = function(x){ sigmaFun0(x) * 
-          (qnorm(1 - p0) - qnorm(1 - (p0 + bmr))) + object$model$curve[[1]](0)}
+          (qnorm(1 - p0) - qnorm(1 - (p0 + bmr))) + object$curve(0)}
     ) 
   } else {
     # BACKGROUND
     if (identical(backgType,"absolute")) {
       if(is.na(backg)){
         stop('backgType = absolute, but backg not supplied')
       }
-      p0 <- pnorm((backg - object$model$curve[[1]](0)) / sigmaFun0(0))
+      p0 <- pnorm((backg - object$curve(0)) / sigmaFun0(0))
     }
     if(identical(backgType, "hybridPercentile")) {
       p0 <- ifelse(is.na(backg),0.1,backg)
@@ -89,15 +95,15 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     bmrScaled <- switch(
       def,
       hybridExc = function(x){ sigmaFun0(x) * 
-          (qnorm(p0) - qnorm(bmr + (1-bmr) * p0)) + object$model$curve[[1]](0)},
+          (qnorm(p0) - qnorm(bmr + (1-bmr) * p0)) + object$curve(0)},
       hybridAdd = function(x){ sigmaFun0(x) * 
-          (qnorm(p0) - qnorm(bmr + p0)) + object$model$curve[[1]](0)}
+          (qnorm(p0) - qnorm(bmr + p0)) + object$curve(0)}
     ) 
   }
 
   # BMD ESTIMATION
-  f0 <- function(x) object$model$curve[[1]](x) - bmrScaled(x)
-  interval0 <- range(object$model$dataList$dose, na.rm = TRUE)
+  f0 <- function(x) object$curve(x) - bmrScaled(x)
+  interval0 <- range(object$dataList$dose, na.rm = TRUE)
   uniroot0 <- try(uniroot(f = f0, interval = interval0), silent = TRUE)
 
   if(inherits(uniroot0, "try-error")){
@@ -113,13 +119,17 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     BMDL <- NA
     BMDU <- NA
   } else {
-    bootDataList <- bootDataGen(object$model, R=R, bootType="nonparametric",aggregated=TRUE)
+    # drc_obj <- eval(substitute(drm(formula0, data = object$data, fct = fct0, type = "continuous", control = drmc(maxIt = 1, noMessage = TRUE)),
+    #                            list(formula0 = object$formula,
+    #                                 fct0 = object$fct
+    #                            )))
+    # bootDataList <- bootDataGen(drc_obj, R=R, bootType="nonparametric",aggregated=FALSE)
+    bootDataList <- bootDataGenHetVar(object, R = R, bootType = bootType)
 
     bmdHetVarBoot <- function(bootData){
-      bootMod <- update(object$model, data = bootData)
-      bootModHetVar <- drmHetVar(bootMod, object$var.formula)
-      bootBmdEst <- bmdHetVar(object = bootModHetVar, bmr = bmr, 
-                              backgType = backgType, backg = backg, def = def, interval = "none", display = FALSE)$Results[,1]
+      bootModHetVar <- drmHetVar(formula = object$formula, var.formula = object$var.formula, data = bootData, fct = object$fct)
+      bootBmdEst <- bmdHetVar(object = bootModHetVar, bmr = bmr, backgType = backgType, backg = backg,
+                              def = def, interval = "none", display = FALSE)$Results[,1]
       bootBmdEst
     }
 
@@ -147,7 +157,7 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
   }
 
   resMat <- matrix(c(bmdEst, BMDL), nrow = 1, ncol = 2, dimnames = list(NULL, c("BMD", "BMDL")))
-  bmrScaled <- matrix(object$model$curve[[1]](bmdEst), nrow = 1, ncol = 1, dimnames = list("", "bmrScaled"))
+  bmrScaled <- matrix(object$curve(bmdEst), nrow = 1, ncol = 1, dimnames = list("", "bmrScaled"))
   bmdInterval <- matrix(c(BMDL, BMDU), nrow = 1, ncol = 2, dimnames = list("", c("Lower", "Upper")))
 
   # GATHER RESULTS
@@ -159,6 +169,6 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
                bmrScaled = bmrScaled,
                interval = bmdInterval,
                model = object)
-  class(resBMD) <- "bmdHetVar"
+  class(resBMD) <- c("bmdHetVar", "bmd")
   invisible(resBMD) 
-}
+}