Fixes to MMRE functionality. bmdBoot now works for drmMMRE models.

Author: JensBaalkilde

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/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/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      
 }
 

R/bootDataGenMMRE.R

@@ -0,0 +1,144 @@
+#' Generate Bootstrap Data for meta-analytic BMD Analysis
+#'
+#' @description
+#' Helper function for `bmdBoot` that generates bootstrap datasets from fitted 
+#' dose-response models. Supports multiple bootstrap methods (nonparametric, 
+#' parametric, semiparametric and wild) and handles different response types 
+#' (binomial, continuous).
+#'
+#' @param object A fitted dose-response model object of class "drcMMRE" (typically from `drmMMRE()`) 
+#'        containing the original data and model specifications
+#' @param R integer; Number of bootstrap replicates to generate (default: 1000)
+#' @param bootType character; Type of bootstrap resampling method. Options are:
+#'        \itemize{
+#'          \item \code{"nonparametric"}: Resamples observations within dose groups
+#'          \item \code{"parametric"}: Generates new data from fitted distributions
+#'          \item \code{"semiparametric"}: Uses fitted values plus resampled residuals
+#'          \item \code{"wild"}: Uses fitted values plus resampled residuals
+#'        }
+#'
+#' @details
+#' The function implements different bootstrap strategies based on data type and method:
+#' 
+#' **Nonparametric Bootstrap:**
+#' \itemize{
+#'   \item \strong{Binomial}: Expands binomial data to individual observations, 
+#'         resamples within dose groups, then optionally re-aggregates
+#'   \item \strong{Continuous/Count}: Resamples observations within dose groups 
+#'         from the original dataset
+#' }
+#' 
+#' **Parametric Bootstrap:**
+#' \itemize{
+#'   \item \strong{Binomial}: Generates new binomial observations using estimated 
+#'         success probabilities, with continuity correction (adds 0.25/0.5) for 
+#'         boundary cases
+#'   \item \strong{Continuous}: Generates normal random variables using 
+#'         dose-specific means and standard deviations from original data
+#' }
+#' 
+#' **Semiparametric Bootstrap:**
+#' \itemize{
+#'   \item \strong{Continuous only}: Uses fitted values plus resampled residuals
+#'   \item \strong{Binomial}: Not supported (throws error)
+#' }
+#'
+#' @return A list of length R containing bootstrap datasets. Each element is a 
+#'         data.frame with the same structure as the original data, containing:
+#'         \itemize{
+#'           \item Dose variable (same name as original)
+#'           \item Response variable(s) (same name(s) as original)
+#'           \item For binomial data: number of successes and total observations
+#'           \item For multi-curve data: curve identifier (if present)
+#'         }
+#'
+#' @section Data Type Handling:
+#' 
+#' **Binomial Data:**
+#' \itemize{
+#'   \item Handles both aggregated (n successes out of N trials) and expanded formats
+#'   \item Preserves dose group structure during resampling
+#'   \item Applies continuity correction in parametric bootstrap
+#' }
+#' 
+#' **Continuous Data:**
+#' \itemize{
+#'   \item Maintains dose group structure
+#'   \item Preserves within-group variability patterns
+#'   \item Uses original data (`origData`) when available
+#' }
+#' 
+#'
+#' @note
+#' This is an internal helper function for `bmdBoot`. It assumes the input object 
+#' has the standard structure from `drmMMRE()` fitting, including components like 
+#' `data`, `origData`, `call`, `type`, etc.
+#' 
+#' **Important considerations:**
+#' \itemize{
+#'   \item Semiparametric bootstrap requires model residuals and fitted values
+#'   \item Parametric bootstrap assumes distributional assumptions are met
+#'   \item Large R values may require substantial memory for complex datasets
+#' }
+#'
+#' @seealso 
+#' \code{\link{bmdBoot}} for the main bootstrap BMD function,
+#'
+#' @examples
+#' \dontrun{
+#' # Typically called internally by bmdBoot, but can be used directly:
+#' 
+#' 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())
+#' boot_data <- bootDataGenMMRE(modMMRE, R = 1000, bootType = "parametric")
+#' 
+#' # Access first bootstrap sample
+#' first_sample <- boot_data[[1]]
+#' } 
+bootDataGenMMRE <- function(object, R=1000, bootType=c("nonparametric", "semiparametric", "parametric", "wild")){
+  if(!inherits(object, "drcMMRE")){
+    stop("bootDataGenMMRE only works for object of type drcMMRE")
+  }
+  
+  R <- as.integer(R)
+  if(is.na(R)){
+    stop("R must be a positive integer")
+  }
+  if(R<=0){
+    stop("R must be a positive integer")
+  }
+  
+  bootType <- match.arg(bootType)
+  
+  # resample for each experiment independently, then collect individual resampled datasets
+  boot.data.per.exp_id <- lapply(names(object$objList), 
+                                 function(x){
+                                   object0 <- object$objList[[x]]
+                                   object0$call$formula <- object$call$formula
+                                   boot.data.per.exp_id0 <- bootDataGen(object0,R=R, bootType=bootType, aggregated = FALSE)
+                                   # add exp_id to data
+                                   boot.data.per.exp_id0 <- lapply(boot.data.per.exp_id0,
+                                                                   function(z) {
+                                                                     z[[as.character(object$call$exp_id)]] <- rep(x, nrow(z))
+                                                                     z
+                                                                     }
+                                   )
+                                   boot.data.per.exp_id0
+                                   }
+                                 )
+  # collect individual resampled datasets
+  tmp.data <- lapply(
+    1:R,
+    function(r) do.call(rbind, lapply(boot.data.per.exp_id, function(list) list[[r]]))
+    )
+  
+  tmp.data
+}
+
+

R/drmMMRE.R

@@ -17,7 +17,7 @@
 #' 
 #' @param formula a symbolic description of the model to be fit of the form
 #' 'response ~ dose'
-#' @param exp_id the name of the column in the data set that specifies the
+#' @param exp_id the column in the data set that specifies the
 #' hierarchical structure of the data
 #' @param data a data frame containing the variables in the model.
 #' @param fct a list with three or more elements specifying the non-linear
@@ -34,7 +34,7 @@
 #' structure of class \code{drcMMRE}.
 #' 
 #' The primary objective is to use this model for benchmark dose estimation
-#' based on dose-response data with a heterogeneous variance structure.
+#' based on dose-response data with a hierarchical variance structure.
 #' @author Signe M. Jensen and Jens Riis Baalkilde
 #' @keywords models nonlinear
 #' @examples
@@ -90,11 +90,6 @@ drmMMRE <- function(formula, exp_id, data, fct, type = c("continuous", "binomial
   drc.objList <- lapply(exp_id_unique,
                         function(exp_id){
                           drm(formula, data = subset(data, data[[exp_id_char]] == exp_id), fct = fct, type = type)
-                          # eval(substitute(drm(formula = formula0, data = subset(data, data[[exp_id_char]] == x), fct = fct0, type = type0),
-                          #            list(formula0 = formula,
-                          #                 data0 = subset(data, data[[exp_id_char]] == x),
-                          #                 fct0 = fct,
-                          #                 type0 = type)))
                         })
   names(drc.objList) <- exp_id_unique
 
@@ -108,7 +103,7 @@ drmMMRE <- function(formula, exp_id, data, fct, type = c("continuous", "binomial
                                           }
                                         ))
 
-  # Fit MA_MA_model
+  # Fit MV_MA_model
   MV_MA_model<-tryCatch( expr={metafor::rma.mv(Estimate
                                       ,(block_diag_matrix+t(block_diag_matrix))/2
                                       , mods = ~0+Coef
@@ -144,13 +139,9 @@ drmMMRE <- function(formula, exp_id, data, fct, type = c("continuous", "binomial
   object$summary <- NULL
   object$start <- NULL
   for(exp_i in exp_id_unique){
-    object$predres[data[[exp_id_char]] == exp_i,] <- drc.objList[[exp_i]]$predres
+    object$predres[data[[exp_id_char]] == exp_i,] <- drc.objList[[as.character(exp_i)]]$predres
   }
   object$call <- call_expr
-  # 
-  # object$df.residual <- NA
-  # object$sumList$df.residual <- NA
-  # object$deriv1 <- NULL
 
   # add MV_MA_model and drc.objList
   object$objList <- drc.objList
@@ -165,7 +156,7 @@ drmMMRE <- function(formula, exp_id, data, fct, type = c("continuous", "binomial
 #'
 #' @description
 #' S3 method to extract the variance-covariance matrix from fitted drcMMRE objects
-#' by delegating to the underlying model-averaged model.
+#' by delegating to the underlying multivariate meta-analytic model.
 #'
 #' @param object An object of class "drcMMRE"
 #' @param ... Additional arguments (currently unused)