Merge pull request #36 from DoseResponse/devel

Version 2.7.5
- Added bootstrap functionality for bmd estimation with drmMMRE models.
- Updated documentation for bmdMA (details for count data were missing for def argument).
- Added delta method confidence interval for bmdHetVar based on numerical differentiation.
- Fixed issue with bmd estimation based on BC model, where maximum of dose-response curve was reached for very small dose values.

Author: JensBaalkilde

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bmd
 Type: Package
 Title: Benchmark dose estimation for dose-response data
-Version: 2.7.4
+Version: 2.7.5
 Date: 2025-03-24
 Author: Signe M.Jensen, Christian Ritz and Jens Riis Baalkilde
 Maintainer: Signe M. Jensen <smj@plen.ku.dk>

NAMESPACE

@@ -16,6 +16,7 @@ S3method(print,bmdOrdinal)
 S3method(print,drcHetVar)
 S3method(print,drcMMRE)
 S3method(print,drcOrdinal)
+S3method(vcov,drcHetVar)
 S3method(vcov,drcMMRE)
 export(BCa)
 export(MACurve)

R/bmd.R

@@ -11,7 +11,7 @@
 #' Details on the implemented definitions and methods can be found in Crump
 #' (2002)
 #' 
-#' @param object object of class \code{drc}
+#' @param object object of class \code{drc} or \code{drcMMRE}
 #' @param bmr numeric value of benchmark response level for which to calculate
 #' the benchmark dose
 #' @param backgType character string specifying how the background level is
@@ -169,6 +169,16 @@
 #' ## BMD using the hybrid method, background risk is 2 SD, hybrid definition using excess risk
 #' bmd(ryegrass.m1, 0.05, backg = 2, backgType = "hybridSD", def = "hybridAdd", display = TRUE)
 #' 
+#' ## BMD on meta-analytic random effects model
+#' set.seed(1)
+#' data0 <- data.frame(x = rep(drcData::ryegrass$conc, 2),
+#'                     y = rep(drcData::ryegrass$rootl, 2) +
+#'                       c(rnorm(n = nrow(drcData::ryegrass), mean = 2, sd = 0.5),
+#'                         rnorm(n = nrow(drcData::ryegrass), mean = 2.7, sd = 0.7)),
+#'                     EXP_ID = rep(as.character(1:2), each = nrow(drcData::ryegrass)))
+#' 
+#' modMMRE <- drmMMRE(y~x, exp_id = EXP_ID, data = data0, fct = LL.4())
+#' bmd(modMMRE, bmr = 0.1, backgType = "modelBased", def = "relative")
 #' 
 #' @export
 bmd<-function(object, bmr, backgType = c("modelBased", "absolute", "hybridSD", "hybridPercentile"),

R/bmd.edfct.R

@@ -213,7 +213,17 @@ bmd.edfct <- function(object){
             expVal <- exp(parmVec[1]*(log(dose)-log(parmVec[4])))
             parmVec[5]*(1+expVal*(1-parmVec[1]))-(parmVec[3]-parmVec[2])*expVal*parmVec[1]/dose
           }
-          maxAt <- uniroot(helpEqn, interval)$root
+          maxAt <- try(uniroot(helpEqn, interval)$root, TRUE)
+          
+          # solution if maxAt fails
+          if(inherits(maxAt, "try-error")){
+            find_maxAt_tries <- 1
+            while(inherits(maxAt, "try-error") & (find_maxAt_tries<10)){
+              interval[1] <- interval[1]/1e2
+              maxAt <- try(uniroot(helpEqn, interval)$root, TRUE)
+              find_maxAt_tries <- find_maxAt_tries+1
+            }
+          }
 
           eqn <- function(dose) {tempVal*(1+exp(parmVec[1]*(log(dose)-log(parmVec[4]))))-(1+parmVec[5]*dose/(parmVec[3]-parmVec[2]))}
           EDp <- uniroot(eqn, lower = maxAt, upper = upper)$root

R/bmdBoot.R

@@ -5,25 +5,27 @@
 #' 
 #' BMDL is defined as the 5th percentile in the bootstrap distribution.
 #' 
-#' Bootstrapping with the argument boot = "nonparametric" is done by sampling
+#' Bootstrapping with the argument bootType = "nonparametric" is done by sampling
 #' with replacement from the original data set. Bootstrapping with the argument
-#' boot = "parametric" is done by sampling from norm(mean(Y_i),sd(Y_0)),
+#' bootType = "parametric" is done by sampling from norm(mean(Y_i),sd(Y_0)),
 #' assuming equal variance between groups, in case of continuous data. For
 #' binomial data, each bootstrap data set is sampled from binom(N_i,Y_i/N_i).
 #' In case of Y_i = 0 or Y_i = N_i shrinkage is used to avoid that the
 #' resampling always produces 0 or 1, respectively. In this case data is
 #' sampled from binom(N_i,(Y_i+1/4)/(N_i+1/2)). Bootstrapping with argument
-#' boot = "semiparametric" is done by sampling with replacement from the
-#' residuals.
+#' bootType = "semiparametric" is done by sampling with replacement from the
+#' residuals. Bootstrapping with argument bootType = "wild" is done by resampling 
+#' with replacement from the residuals multiplied by a random sign (-1 or +1).
 #' 
-#' All sampling is made within dose groups.
+#' All sampling is made within dose groups. When a meta-analytic random effects 
+#' model is supplied, sampling is made within dose groups within each experiment.
 #' 
-#' @param object object of class \code{drc}
+#' @param object object of class \code{drc} or \code{drcMMRE}
 #' @param bmr numeric value of benchmark response level for which to calculate
 #' the benchmark dose
 #' @param R number of bootstrap samples. Default is 1000
 #' @param bootType character string specifying type of bootstrap used.
-#' "nonparametric" (default), "semiparametric" or "parametric". "Semiparametric
+#' "nonparametric" (default), "semiparametric", "parametric" or "wild". "semiparametric" and "wild"
 #' is only available for continuous data and "nonparametric" is at present the
 #' only option for count data. See details below
 #' @param bmdType Type of estimate for BMD. Default is "orig" the bmd estimate
@@ -102,10 +104,10 @@
 #' and adjusted)
 #' @param display logical. If TRUE the results are displayed; otherwise they
 #' are not
-#' @param level numeric value specifying the levle of the confidence interval
+#' @param level numeric value specifying the level of the confidence interval
 #' underlying BMDL. Default is 0.95
-#' @return A list of three elements: Results contain the estimated BMD and
-#' BMDL, bootEst is a vector af all of the estimated BMD values from each
+#' @return A list of four elements: Results contain the estimated BMD and
+#' BMDL, Boot.samples.used is the number of samples used, bootEst is a vector af all of the estimated BMD values from each
 #' bootstrap sample, Interval gives BMDL and BMDU, which is identical to the
 #' confidence interval for the percentile interval approach.
 #' @author Signe M. Jensen
@@ -125,6 +127,17 @@
 #' ## BMD from the same definitions but using parametric bootstrap
 #' bmdBoot(ryegrass.m1, 0.05, backgType = "hybridSD", def = "hybridAdd", bootType="parametric",R = 50)
 #' 
+#' ## BMD on meta-analytic random effects model
+#' set.seed(1)
+#' data0 <- data.frame(x = rep(drcData::ryegrass$conc, 2),
+#'                     y = rep(drcData::ryegrass$rootl, 2) +
+#'                       c(rnorm(n = nrow(drcData::ryegrass), mean = 2, sd = 0.5),
+#'                         rnorm(n = nrow(drcData::ryegrass), mean = 2.7, sd = 0.7)),
+#'                     EXP_ID = rep(as.character(1:2), each = nrow(drcData::ryegrass)))
+#' 
+#' modMMRE <- drmMMRE(y~x, exp_id = EXP_ID, data = data0, fct = LL.4())
+#' bmdBoot(modMMRE, bmr = 0.1, backgType = "modelBased", def = "relative", R = 50, bootType = "wild")
+#' 
 #' @export
 bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "orig",
                     backgType = c("modelBased", "absolute", "hybridSD", "hybridPercentile"),
@@ -183,7 +196,16 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
               interval = "delta", display = FALSE)
 
   get.drm.list <- function(tmp.data){
-    if(ncol(object$parmMat) == 1){
+    if(inherits(object, "drcMMRE")){
+      drm.list.tmp <- lapply(tmp.data, function(x){
+        try(eval(substitute(drmMMRE(formula0, exp_id = exp_id0, data = x, type = object$type, fct = object[["fct"]]),
+                            list(formula0 = object$call$formula,
+                                 exp_id0 = object$call$exp_id)
+        )), TRUE)
+      }
+      )
+    }
+    else 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"]],
                                 control = drmc(noMessage = TRUE)),
@@ -193,9 +215,6 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
       )
     } 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))]]
-        # }
         try(
           eval(substitute(drm(object$call$formula, weights = weights0, curveid = curveid0,
                               data = x, type = object$type, fct = object$fct, control = drmc(noMessage = TRUE)),
@@ -206,9 +225,6 @@ 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))]]
-        # }
         try(
           eval(substitute(drm(formula0, weights = weights0, curveid = curveid0, pmodels = pmodels0,
                               data = x, type = object$type, fct = object$fct, control = drmc(noMessage = TRUE)),
@@ -233,11 +249,14 @@ bmdBoot <- function(object, bmr, R=1000, bootType="nonparametric", bmdType = "or
   }
 
   if (object$type %in% c("binomial","continuous")) {
-    
-    tmp.data <- bootDataGen(object,R,bootType,aggregated=FALSE)
+    if(inherits(object, "drcMMRE")){
+      tmp.data <- bootDataGenMMRE(object,R,bootType)
+    } else {
+      tmp.data <- bootDataGen(object,R,bootType,aggregated=FALSE)
+    }
 
     drm.list.tmp <- get.drm.list(tmp.data)
-    list.condition <- sapply(drm.list.tmp, function(x) class(x)=="drc")
+    list.condition <- sapply(drm.list.tmp, inherits, "drc")
     drm.list  <- drm.list.tmp[list.condition]
 
     bmd.list.try <- lapply(drm.list,get.bmd)

R/bmdHetVar.R

@@ -49,10 +49,12 @@
 #' is the normal distribution function and sigma(BMD) is the SD at the
 #' benchmark dose.
 #' @param interval character string specifying the type of confidence interval
-#' to use: "boot" (default) or "none"
+#' to use: "boot" (default), "delta" or "none"
 #' 
 #' "boot" - BMDL is based on percentile bootstrapping.
 #' 
+#' "delta" - BMDL is based on the delta method. 
+#' 
 #' "none" - no confidence interval is computed.
 #' @param R number of bootstrap samples. Ignored if \code{interval = "none"}
 #' @param level numeric value specifying the levle of the confidence interval
@@ -113,7 +115,7 @@
 #'           def = "hybridExc", R = 50, level = 0.95, progressInfo = TRUE, display = TRUE)
 #' 
 #' @export
-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){
+bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybridPercentile"), backg = NA, def = c("hybridExc", "hybridAdd"), interval = c("boot", "delta", "none"), R = 1000, level = 0.95, bootType = "nonparametric", progressInfo = TRUE, display = TRUE){
   ### Assertions ###
   # object
   if(!inherits(object,"drcHetVar")){
@@ -158,14 +160,19 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
 
   level <- 1-2*(1-level)
 
+  # Model parameters
+  curveParInd <- 1:length(object$curvePar)
+  sigmaParInd <- (length(object$curvePar)+1):(length(object$curvePar)+length(object$sigmaPar))
+  
   # SLOPE
   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
-  sigmaFun0 <- object$sigmaFun # sigmaFun(object, var.formula)
+  # functions
+  curveFun0 <- object$funList$curveFun
+  sigmaFun0 <- object$funList$sigmaFun # sigmaFun(object, var.formula)
 
   # bmrScaled
   if(slope == "increasing"){
@@ -174,50 +181,51 @@ 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$curve(0)) / sigmaFun0(0))
+      p0 <- function(par) 1 - pnorm((backg - curveFun0(0, par[curveParInd])) / sigmaFun0(0, sigmaPar = par[sigmaParInd], curvePar = par[curveParInd]))
     }
     if(identical(backgType, "hybridPercentile")) {
-      p0 <- ifelse(is.na(backg),1-0.9,1-backg)
+      p0 <- function(par) ifelse(is.na(backg),1-0.9,1-backg)
     }
     if (identical(backgType,"hybridSD")) {
-      p0 <- ifelse(is.na(backg), 1-pnorm(2), 1-pnorm(backg))
+      p0 <- function(par) ifelse(is.na(backg), 1-pnorm(2), 1-pnorm(backg))
     }
 
     # BMRSCALED
     bmrScaled <- switch(
       def,
-      hybridExc = function(x){ sigmaFun0(x) * 
-          (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$curve(0)}
+      hybridExc = function(x, par){ sigmaFun0(x) * 
+          (qnorm(1 - p0(par)) - qnorm(1 - p0(par) - (1 - p0)*bmr)) + curveFun0(0, par[curveParInd])},
+      hybridAdd = function(x, par){ sigmaFun0(x) * 
+          (qnorm(1 - p0(par)) - qnorm(1 - (p0(par) + bmr))) + curveFun0(0, par[curveParInd])}
     ) 
   } else {
     # BACKGROUND
     if (identical(backgType,"absolute")) {
       if(is.na(backg)){
         stop('backgType = absolute, but backg not supplied')
       }
-      p0 <- pnorm((backg - object$curve(0)) / sigmaFun0(0))
+      p0 <- function(par) pnorm((backg - curveFun0(0, par[curveParInd])) / sigmaFun0(0, sigmaPar = par[sigmaParInd], curvePar = par[curveParInd]))
     }
     if(identical(backgType, "hybridPercentile")) {
-      p0 <- ifelse(is.na(backg),0.1,backg)
+      p0 <- function(par) ifelse(is.na(backg),0.1,backg)
     }
     if (identical(backgType,"hybridSD")) {
-      p0 <- ifelse(is.na(backg), pnorm(-2), pnorm(-backg))
+      p0 <- function(par) ifelse(is.na(backg), pnorm(-2), pnorm(-backg))
     }
 
     # BMRSCALED
     bmrScaled <- switch(
       def,
-      hybridExc = function(x){ sigmaFun0(x) * 
-          (qnorm(p0) - qnorm(bmr + (1-bmr) * p0)) + object$curve(0)},
-      hybridAdd = function(x){ sigmaFun0(x) * 
-          (qnorm(p0) - qnorm(bmr + p0)) + object$curve(0)}
+      hybridExc = function(x, par){ sigmaFun0(x, sigmaPar = par[sigmaParInd], curvePar = par[curveParInd]) * 
+          (qnorm(p0(par)) - qnorm(bmr + (1-bmr) * p0(par))) + curveFun0(0, par[curveParInd])},
+      hybridAdd = function(x, par){ sigmaFun0(x, sigmaPar = par[sigmaParInd], curvePar = par[curveParInd]) * 
+          (qnorm(p0(par)) - qnorm(bmr + p0(par))) + curveFun0(0, par[curveParInd])}
     ) 
   }
 
   # BMD ESTIMATION
-  f0 <- function(x) object$curve(x) - bmrScaled(x)
+  par0 <- c(object$curvePar, object$sigmaPar)
+  f0 <- function(x) curveFun0(x, par0[curveParInd]) - bmrScaled(x, par0)
   interval0 <- range(object$dataList$dose, na.rm = TRUE)
   uniroot0 <- try(uniroot(f = f0, interval = interval0), silent = TRUE)
 
@@ -228,17 +236,34 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     bmdEst <- uniroot0$root
   }
 
+  
   # INTERVAL
   interval <- match.arg(interval)
   if(identical(interval, "none")){
     BMDL <- NA
     BMDU <- NA
-  } else {
-    # 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)
+    SDbmd <- NA
+  } else if(identical(interval, "delta")){
+    if(!requireNamespace("numDeriv")){
+      stop('package "numDeriv" must be installed to compute delta confidence intervals with bmdHetVar')
+    }
+    
+    getBmdEst <- function(par){
+      f0 <- function(x) curveFun0(x, par[curveParInd]) - bmrScaled(x, par)
+      interval0 <- range(object$dataList$dose, na.rm = TRUE)
+      uniroot0 <- try(uniroot(f = f0, interval = interval0), silent = TRUE)
+      bmdEst <- as.numeric(uniroot0$root)
+    }
+    
+    dBmd <- try(numDeriv::grad(getBmdEst, par0))
+    
+    Vbmd <- dBmd %*% vcov(object) %*% dBmd
+    SDbmd <- sqrt(Vbmd)
+    
+    BMDL <- qnorm(c(1-level)/2, mean = bmdEst, sd = SDbmd) 
+    BMDU <- qnorm(1-(1-level)/2, mean = bmdEst, sd = SDbmd) 
+    
+  } else if(identical(interval, "boot")) {
     bootDataList <- bootDataGenHetVar(object, R = R, bootType = bootType)
 
     bmdHetVarBoot <- function(bootData){
@@ -265,9 +290,11 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
     if(length(boot0) == 0){ 
       BMDL <- NA 
       BMDU <- NA
+      SDbmd <- NA
     } else {
       BMDL <- quantile(boot0,p=c(1-level), na.rm = TRUE) # ABC percentile lims.  
       BMDU <- quantile(boot0,p=c(level), na.rm = TRUE)
+      SDbmd <- sd(boot0)
     }
   }
 
@@ -282,6 +309,7 @@ bmdHetVar <- function(object, bmr, backgType = c("absolute", "hybridSD", "hybrid
 
   resBMD<-list(Results = resMat,
                bmrScaled = bmrScaled,
+               SE = SDbmd,
                interval = bmdInterval,
                model = object)
   class(resBMD) <- c("bmdHetVar", "bmd")

R/bmdMA.R

@@ -57,10 +57,10 @@
 #' dose-response models), 2 SD for "hybridSD" background and 0.9 for
 #' "hybridpercentile"
 #' @param def character string specifying the definition of the benchmark dose
-#' to use in the calculations. "excess" , "additional" and "point" are for
-#' binomial response whereas "relative", "extra", "added", "hybridExc" (excess
+#' to use in the calculations. "excess", "additional" and "point" are for
+#' binomial response. "relative", "extra", "added", "hybridExc" (excess
 #' hybrid), "hybridAdd" (additional hybrid), and "point" are for continuous
-#' response
+#' response. "relative", "extra", and "point" are for count response data.
 #' 
 #' "excess" - BMR is defined as: BMR = (f(BMD) - p0)/(1 - p0).  Works for
 #' binomial response. BMR should be between 0 and 1.

R/bootDataGen.R

@@ -239,6 +239,23 @@ bootDataGen <- function(object, R=1000, bootType="nonparametric",aggregated=TRUE
       }
     }
   }
+  else if(bootType=="wild"){
+    if(object$type=="binomial"){
+      stop(paste("wild is not possible for binomial data", sep=""))
+    }
+    if(object$type=="continuous"){
+      data.st<-object$data
+      
+      tmp.data <- list()
+      for(i in 1:R){
+        random.sign <- sample(c(-1,1), size = length(resid(object)), replace = TRUE)
+        sampled <- data.frame(y = fitted(object)+sample(resid(object),replace=TRUE)*random.sign, 
+                              dose = object$data[,as.character(object$call$formula[[3]])])
+        colnames(sampled) <- c(as.character(object$call$formula[[2]]), as.character(object$call$formula[[3]]))
+        tmp.data[[i]] <- sampled
+      }
+    }
+  }
 tmp.data      
 }