Merge pull request #35 from DoseResponse/devel

Version 2.7.4

Author: JensBaalkilde

.DS_Store

@@ -4,4 +4,6 @@ cache/
 README.Rmd
 README.md
 README_files/
-^\\.github$
+^\.github$
+docs/
+example_extracts/

.Rbuildignore

@@ -2,4 +2,8 @@
 .Rhistory
 .RData
 .Ruserdata
-cache/
+cache/
+bmd.Rproj
+inst/doc
+inst/dev
+example_extracts/

.gitignore

@@ -1,7 +1,7 @@
 Package: bmd
 Type: Package
 Title: Benchmark dose estimation for dose-response data
-Version: 2.7.3
+Version: 2.7.4
 Date: 2025-03-24
 Author: Signe M.Jensen, Christian Ritz and Jens Riis Baalkilde
 Maintainer: Signe M. Jensen <smj@plen.ku.dk>
@@ -25,14 +25,18 @@ Suggests:
     tidyr,
     numDeriv,
     drcData,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    knitr,
+    rmarkdown
 Remotes:
     doseresponse/drc,
     doseresponse/drcData
 License: file LICENSE
 Encoding: UTF-8
 LazyData: true
 Config/testthat/edition: 3
+Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 Depends: 
     R (>= 3.5)
+VignetteBuilder: knitr

DESCRIPTION

@@ -1,29 +1,78 @@
-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,
-       qplotDrc, qplotBmd, MACurve, monotonicityTest, trendTest, bmdHetVar, drmHetVar, drmMMRE, bmdHetVarMA)
+# Generated by roxygen2: do not edit by hand
 
-## 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(vcov, drcMMRE)
-S3method(print, drcMMRE)
-S3method(coef, drcHetVar)
+S3method(AIC,drcHetVar)
+S3method(AIC,drcOrdinal)
+S3method(BIC,drcHetVar)
+S3method(BIC,drcOrdinal)
+S3method(coef,drcHetVar)
+S3method(logLik,drcHetVar)
+S3method(logLik,drcOrdinal)
+S3method(plot,bmd)
+S3method(plot,drcHetVar)
+S3method(plot,drcOrdinal)
+S3method(print,bmd)
+S3method(print,bmdHetVar)
+S3method(print,bmdOrdinal)
+S3method(print,drcHetVar)
+S3method(print,drcMMRE)
+S3method(print,drcOrdinal)
+S3method(vcov,drcMMRE)
+export(BCa)
+export(MACurve)
+export(PAV)
+export(bmd)
+export(bmdBoot)
+export(bmdHetVar)
+export(bmdHetVarMA)
+export(bmdIso)
+export(bmdIsoBoot)
+export(bmdMA)
+export(bmdMACurve)
+export(bmdOrdinal)
+export(bmdOrdinalMA)
+export(drmHetVar)
+export(drmMMRE)
+export(drmOrdinal)
+export(getStackingWeights)
+export(monotonicityTest)
+export(qplotBmd)
+export(qplotDrc)
+export(trendTest)
+import(dplyr)
+import(drc)
+import(ggplot2)
+importFrom(graphics,lines)
+importFrom(stats,AIC)
+importFrom(stats,BIC)
+importFrom(stats,aggregate)
+importFrom(stats,approx)
+importFrom(stats,as.formula)
+importFrom(stats,coef)
+importFrom(stats,complete.cases)
+importFrom(stats,confint)
+importFrom(stats,constrOptim)
+importFrom(stats,df.residual)
+importFrom(stats,dnorm)
+importFrom(stats,fitted)
+importFrom(stats,lm)
+importFrom(stats,logLik)
+importFrom(stats,model.frame)
+importFrom(stats,model.matrix)
+importFrom(stats,optim)
+importFrom(stats,pnorm)
+importFrom(stats,predict)
+importFrom(stats,qchisq)
+importFrom(stats,qnorm)
+importFrom(stats,qt)
+importFrom(stats,quantile)
+importFrom(stats,rbinom)
+importFrom(stats,resid)
+importFrom(stats,residuals)
+importFrom(stats,rnorm)
+importFrom(stats,sd)
+importFrom(stats,uniroot)
+importFrom(stats,update)
+importFrom(stats,var)
+importFrom(stats,vcov)
+importFrom(utils,setTxtProgressBar)
+importFrom(utils,txtProgressBar)

NAMESPACE

@@ -1,3 +1,14 @@
+
+#' AIC Method for drcOrdinal Objects
+#'
+#' Calculates the Akaike Information Criterion for drcOrdinal model objects.
+#'
+#' @param object A drcOrdinal model object
+#' @param ... Additional arguments (not used)
+#' @param k The penalty per parameter to be used; default is 2
+#'
+#' @return The AIC value
+#' @export
 AIC.drcOrdinal <- function(object, ..., k = 2) {
   dots <- list(...)
   if (!is.null(dots$epsilon)){
@@ -8,4 +19,4 @@ AIC.drcOrdinal <- function(object, ..., k = 2) {
 
   n.parameters <- sum(sapply(object$drmList, function(mod) length(mod$coefficients)))
   -2 * logLik(object, epsilon = epsilon) + k * n.parameters
-}
+}

R/AIC.drcOrdinal.R

@@ -1,25 +1,82 @@
+#' Calculate BCa (Bias-Corrected and Accelerated) Bootstrap Confidence Interval
+#'
+#' @description
+#' This function calculates the BCa confidence interval for a given statistic using bootstrap and jackknife estimates.
+#' 
+#' @details
+#' It computes the BCa confidence interval for a statistic using bootstrap and jackknife estimates.
+#' The BCa method adjusts for both bias and skewness in the bootstrap distribution and generally
+#' provides better coverage than standard bootstrap confidence intervals.
+#'
+#' @param obs numeric; The observed value of the statistic computed from the original sample
+#' @param data data.frame or matrix; The original dataset from which bootstrap samples are drawn
+#' @param bootSample numeric vector; Bootstrap replicates of the statistic. Must be of length R 
+#'        where R is the number of bootstrap replicates
+#' @param bootjack numeric vector; Jackknife replicates of the statistic. Must be of length n 
+#'        where n is the number of observations in data
+#' @param level numeric; Confidence level between 0 and 1 (e.g., 0.95 for 95% confidence interval)
+#'
+#' @return A numeric vector of length 2 containing:
+#' \itemize{
+#'   \item lower: The lower confidence limit
+#'   \item upper: The upper confidence limit
+#' }
+#'
+#'
+#' @examples \dontrun{
+#' # Example with mean as the statistic
+#' data <- rnorm(100)
+#' obs_mean <- mean(data)
+#' 
+#' # Generate bootstrap samples
+#' R <- 1000
+#' boot_means <- replicate(R, mean(sample(data, replace = TRUE)))
+#' 
+#' # Generate jackknife samples
+#' jack_means <- sapply(1:length(data), 
+#'                      function(i) mean(data[-i]))
+#' 
+#' # Calculate 95% CI
+#' BCa(obs_mean, data, boot_means, jack_means, 0.95)
+#' }
+#'
+#' @references
+#' Efron, B. (1987). Better Bootstrap Confidence Intervals. 
+#' _Journal of the American Statistical Association_, 82(397), 171-185.
+#'
+#' DiCiccio, T. J., & Efron, B. (1996). Bootstrap confidence intervals.
+#' _Statistical Science_, 11(3), 189-212.
+#'
+#' @export
 BCa <- function(obs, data, bootSample, bootjack, level){
-R <- length(bootSample)
-b <- qnorm((sum(bootSample > obs)+sum(bootSample==obs)/2)/R)
+  R <- length(bootSample)
+  b <- qnorm((sum(bootSample > obs)+sum(bootSample==obs)/2)/R)
 
-n <- nrow(data) 
-n1 <- n-1 
-obsn <- obs*n
-pv <- i <- 0 
-while(i < n){
-  i = i+1 
-pv[i] = obsn-n1*bootjack[i]
+  n <- nrow(data) 
+  n1 <- n-1 
+  obsn <- obs*n
+  pv <- i <- 0 
+  while(i < n){
+    i = i+1 
+    pv[i] = obsn-n1*bootjack[i]
+  }
+  pv<-pv[!is.na(pv)]
+  
+  je <- mean(pv)-pv
+  a <- sum(je^3)/(6*sum(je^2))^(3/2)
+  
+  alpha <- (1-level)*2 
+  z <- qnorm(c(alpha/2,1-alpha/2)) # Std. norm. limits
+  p <- pnorm((z-b)/(1-a*(z-b))-b) # correct & convert to proportions
+  
+  quantile(bootSample,p=p) # ABC percentile lims.      
 }
-pv<-pv[!is.na(pv)]
 
-je <- mean(pv)-pv
-a <- sum(je^3)/(6*sum(je^2))^(3/2)
 
-alpha <- (1-level)*2 
-z <- qnorm(c(alpha/2,1-alpha/2)) # Std. norm. limits
-p <- pnorm((z-b)/(1-a*(z-b))-b) # correct & convert to proportions
 
-quantile(bootSample,p=p) # ABC percentile lims.      
-}
+
+
+
+
 
 

R/BCa.R

@@ -1,3 +1,88 @@
+#' BIC Method for drcOrdinal Objects
+#'
+#' @description
+#' S3 method to calculate the Bayesian Information Criterion (BIC) for fitted 
+#' ordinal dose-response models. This method extends the generic `BIC()` function 
+#' to handle `drcOrdinal` objects which contain multiple fitted dose-response models 
+#' for different response categories.
+#'
+#' @param object An object of class "drcOrdinal" containing fitted ordinal 
+#'        dose-response models
+#' @param ... Additional arguments passed to the method. Currently supports:
+#'        \itemize{
+#'          \item \code{epsilon}: Small positive value to avoid log(0) in likelihood 
+#'                calculations (default: 1e-16)
+#'        }
+#'
+#' @details
+#' The BIC is calculated using the standard formula:
+#' \deqn{BIC = k \log(n) - 2 \log(L)}{BIC = k * log(n) - 2 * log(L)}
+#' where:
+#' \itemize{
+#'   \item \eqn{k}{k} is the total number of parameters across all models in the ordinal fit
+#'   \item \eqn{n}{n} is the total number of observations (sum of weights)
+#'   \item \eqn{L}{L} is the likelihood of the fitted model
+#' }
+#' 
+#' For ordinal dose-response models, the total number of parameters is the sum 
+#' of parameters from all individual models in the `drmList` component, and the 
+#' sample size is calculated as the sum of weights in the data.
+#' 
+#' The `epsilon` parameter is used in the `logLik()` calculation to prevent 
+#' numerical issues when probabilities approach zero.
+#'
+#' @return A numeric value representing the BIC of the fitted ordinal dose-response model.
+#'         Lower values indicate better model fit when comparing models.
+#'
+#' @section Model Comparison:
+#' BIC penalizes model complexity more heavily than AIC, making it useful for:
+#' \itemize{
+#'   \item Comparing ordinal dose-response models with different numbers of parameters
+#'   \item Model selection in situations where simpler models are preferred
+#'   \item Avoiding overfitting in dose-response modeling
+#' }
+#'
+#' @note
+#' This method assumes that the `drcOrdinal` object contains:
+#' \itemize{
+#'   \item A `drmList` component with fitted dose-response models
+#'   \item A `data` component with the original data including weights
+#'   \item A `weights` component specifying the column name for observation weights
+#' }
+#'
+#' @seealso 
+#' \code{\link{BIC}} for the generic BIC function,
+#' \code{\link{AIC.drcOrdinal}} for the corresponding AIC method,
+#' \code{\link{logLik.drcOrdinal}} for the log-likelihood method
+#'
+#' @examples
+#' \dontrun{
+#' # Assuming you have a fitted drcOrdinal object
+#' ordinal_model <- drm(response ~ dose, 
+#'                      weights = count,
+#'                      data = ordinal_data, 
+#'                      fct = cumulative.logit())
+#' 
+#' # Calculate BIC
+#' bic_value <- BIC(ordinal_model)
+#' 
+#' # Compare models
+#' model1_bic <- BIC(ordinal_model1)
+#' model2_bic <- BIC(ordinal_model2)
+#' 
+#' # Lower BIC indicates better model
+#' if(model1_bic < model2_bic) {
+#'   cat("Model 1 is preferred\n")
+#' } else {
+#'   cat("Model 2 is preferred\n")
+#' }
+#' 
+#' # Use custom epsilon for numerical stability
+#' bic_custom <- BIC(ordinal_model, epsilon = 1e-12)
+#' }
+#'
+#' @method BIC drcOrdinal
+#' @export
 BIC.drcOrdinal <- function(object, ...){
   dots <- list(...)
   if (!is.null(dots$epsilon)){