Merge pull request #31 from JensBaalkilde/master

Version 2.7.3:
Updates to drmHetVar, bmdHetVar bmdHetVarMA and general clean-up.

Author: JensBaalkilde

.DS_Store

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

.Rbuildignore

@@ -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

.github/workflows/r.yml

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

DESCRIPTION

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

NAMESPACE

@@ -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/AIC.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/BIC.drcOrdinal.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/PAV.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/bmd.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)