diff --git a/.github/workflows/r.yml b/.github/workflows/r.yml
index ad77320..347ecba 100644
--- a/.github/workflows/r.yml
+++ b/.github/workflows/r.yml
@@ -44,13 +44,13 @@ jobs:
           http-user-agent: ${{ matrix.config.http-user-agent }}
           use-public-rspm: true
 
-      - name: Cache R packages
-        uses: actions/cache@v3
-        with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-r-${{ matrix.config.r }}-${{ hashFiles('DESCRIPTION') }}
-          restore-keys: |
-            ${{ runner.os }}-r-${{ matrix.config.r }}-
+#      - name: Cache R packages
+#        uses: actions/cache@v3
+#        with:
+#          path: ${{ env.R_LIBS_USER }}
+#          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", 
diff --git a/DESCRIPTION b/DESCRIPTION
index 1124bb5..5970216 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: bmd
 Type: Package
 Title: Benchmark dose estimation for dose-response data
-Version: 2.7.5
+Version: 2.7.6
 Date: 2025-03-24
 Author: Signe M.Jensen, Christian Ritz and Jens Riis Baalkilde
 Maintainer: Signe M. Jensen <smj@plen.ku.dk>
@@ -36,7 +36,7 @@ Encoding: UTF-8
 LazyData: true
 Config/testthat/edition: 3
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
 Depends: 
     R (>= 3.5)
 VignetteBuilder: knitr
+Config/roxygen2/version: 8.1.0
diff --git a/NAMESPACE b/NAMESPACE
index b0daac8..b6cf581 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -43,37 +43,42 @@ 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)
+importFrom(stats,
+  AIC,
+  BIC,
+  aggregate,
+  approx,
+  as.formula,
+  coef,
+  complete.cases,
+  confint,
+  constrOptim,
+  df.residual,
+  dnorm,
+  fitted,
+  lm,
+  logLik,
+  model.frame,
+  model.matrix,
+  optim,
+  pnorm,
+  predict,
+  qchisq,
+  qnorm,
+  qt,
+  quantile,
+  rbinom,
+  resid,
+  residuals,
+  rnorm,
+  sd,
+  uniroot,
+  update,
+  var,
+  vcov
+)
+importFrom(utils,
+  packageVersion,
+  setTxtProgressBar,
+  txtProgressBar
+)
diff --git a/R/bmd-package.R b/R/bmd-package.R
index dc939fc..c70cc7e 100644
--- a/R/bmd-package.R
+++ b/R/bmd-package.R
@@ -11,7 +11,7 @@
 #' @importFrom stats lm model.frame model.matrix optim pnorm predict
 #' @importFrom stats qchisq qnorm qt quantile rbinom resid residuals
 #' @importFrom stats rnorm sd uniroot update var vcov AIC BIC logLik
-#' @importFrom utils setTxtProgressBar txtProgressBar
+#' @importFrom utils setTxtProgressBar txtProgressBar packageVersion
 #'
 #' @name bmd-package
 #' @aliases bmd-package
diff --git a/R/bmdMA.R b/R/bmdMA.R
index a067435..b45b469 100644
--- a/R/bmdMA.R
+++ b/R/bmdMA.R
@@ -285,8 +285,8 @@ bmdMA <- function(modelList, modelWeights, bmr,
       modelWeights0 <- modelWeights
     }
     
-    if(identical(modelList[[1]]$type,"continuous")){
-      my.fun<-function(x,y){drm(y$call$formula, data = x, fct = y[["fct"]])}
+    if(identical(modelList[[1]]$type,"continuous") | identical(modelList[[1]]$type,"Poisson")){
+      my.fun<-function(x,y){drm(y$call$formula, data = x, fct = y[["fct"]], type = y[["type"]])}
       
       if(identical(type,"Kang")){
         maBMD <- sum(modelWeights0 * sapply(bmdList, function(x){x$Results[,1]}))
@@ -321,7 +321,7 @@ bmdMA <- function(modelList, modelWeights, bmr,
         
         bmdMAboot <- function(data){
           bootModelList <- lapply(modelList, function(model) try(
-            eval(substitute(drm(formula = formula0, data = data, fct = model$fct, weights = weights0, start = start0,
+            eval(substitute(drm(formula = formula0, data = data, fct = model$fct, type = model$type, weights = weights0, start = start0,
                                 control = drmc(noMessage = TRUE)),
                             list(formula0 = model$call$formula, 
                                  weights0 = model$call$weights,
@@ -723,6 +723,7 @@ bmdMA <- function(modelList, modelWeights, bmr,
         }
       }
     }
+    
   } 
   
   if (nCurves > 1){
@@ -755,7 +756,7 @@ bmdMA <- function(modelList, modelWeights, bmr,
         modelWeights0 <- modelWeights
       }
       
-      if(identical(modelList[[1]]$type,"continuous")){
+      if(identical(modelList[[1]]$type,"continuous") | identical(modelList[[1]]$type,"Poisson")){
         if(identical(type,"Kang")){
           maBMD <- colSums(modelWeights0 * t(sapply(bmdList, function(x) x$Results[,1])))
           maBMDL <- colSums(modelWeights0 * t(sapply(bmdList, function(x) x$interval[,1])))
@@ -1274,6 +1275,7 @@ bmdMA <- function(modelList, modelWeights, bmr,
           }
         }
       }
+      
     } else {
       # CURVES FITTED INDEPENDENTLY
       modelListList <- lapply(1:length(modelList[[1]]$objList), function(i) lapply(modelList, function(object) object$objList[[i]]))
diff --git a/R/getStackingWeights.R b/R/getStackingWeights.R
index 41eb0f3..775729d 100644
--- a/R/getStackingWeights.R
+++ b/R/getStackingWeights.R
@@ -57,8 +57,8 @@ computeWeightsFromSplit <- function(trainData, validateData, modelList){
     objective <- CVXR::Minimize(sum((predMatrix %*% alphaHat - validateData[[as.character(modelList[[1]]$call$formula[[2]][[2]])]])^2))
   }
   problem <- CVXR::Problem(objective, constraints = list(alphaHat <= 1, alphaHat >= 0,sum(alphaHat) == 1))
-  result <- CVXR::solve(problem)
-  res <- result$getValue(alphaHat)
+  opt_value <- CVXR::psolve(problem) # solves the constrained optimisation problem
+  res <- CVXR::value(alphaHat) # returns the value of alpha under the constrained optimisation
   
   # Initialise weights to zero
   tmpWeights <- numeric(length(modelList))
@@ -245,6 +245,10 @@ getStackingWeights <- function(modelList, nSplits = 2){
     stop('package "CVXR" must be installed to estimate stacking weights')
   }
   
+  if (packageVersion("CVXR") < "1.8.0") {
+    stop("Please update CVXR. Version 1.8.0 or newer is required.")
+  }
+  
   if(nSplits %in% c("LOO")){
     nSplits <- ifelse(modelList[[1]]$type == "binomial", sum(modelList[[1]]$data$weights), modelList[[1]]$sumList$lenData)
   } else if(!is.numeric(nSplits)){
diff --git a/R/qplotDrc.R b/R/qplotDrc.R
index 9428557..3025cea 100644
--- a/R/qplotDrc.R
+++ b/R/qplotDrc.R
@@ -395,7 +395,7 @@ qplotDrc <- function(x, add = FALSE, level = NULL, type = c("average", "all", "b
       obsLayer +
       scale_x_continuous(trans = xtrans, limits = xLimits) +
       scale_y_continuous(trans = ytrans, limits = yLimits) +
-      labs(x = xlab, y = ylab, col = "", fill = "", shape = "", linetype = "")
+      labs(x = xlab, y = ylab, col = NULL, fill = NULL, shape = NULL, linetype = NULL)
   } else {
     list(
       confBandLayer = confBandLayer,
diff --git a/man/bmdMA.Rd b/man/bmdMA.Rd
index f9d9f8d..7cff06e 100644
--- a/man/bmdMA.Rd
+++ b/man/bmdMA.Rd
@@ -65,10 +65,10 @@ dose-response models), 2 SD for "hybridSD" background and 0.9 for
 "hybridpercentile"}
 
 \item{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.
diff --git a/tests/testthat/test-bmd.R b/tests/testthat/test-bmd.R
index dd9ade0..b0219a6 100644
--- a/tests/testthat/test-bmd.R
+++ b/tests/testthat/test-bmd.R
@@ -1523,8 +1523,8 @@ test_that("bmd function works on drcMMRE object", {
   bmdMMRE <- bmd(modMMRE, bmr = 0.1, backgType = "modelBased", def = "relative", display = FALSE)
   
   expect_true(all(!is.na(bmdMMRE$Results[, "BMD"])))
-  expect_equal(bmdMMRE$Results[, "BMD"], 1.66913593445629)
-  expect_equal(bmdMMRE$bmrScaled[,1], 9.15712352078559)
+  expect_equal(bmdMMRE$Results[, "BMD"], 1.66913593445629, tolerance = 1e-6)
+  expect_equal(bmdMMRE$bmrScaled[,1], 9.15712352078559, tolerance = 1e-6)
   expect_equal(unname(bmdMMRE$bmrScaled[,1]), drop(modMMRE$curve[[1]](bmdMMRE$Results[, "BMD"])))
   expect_equal(bmdMMRE$interval[1,], c(Lower = 1.3166277025622, Upper = 2.02164416635037), tolerance = 1e-4)
   expect_equal(bmdMMRE$SE[,"SE"], 0.214309787885148, tolerance = 1e-4) 
diff --git a/tests/testthat/test-bmdMA.R b/tests/testthat/test-bmdMA.R
index 691fb7b..c426dd4 100644
--- a/tests/testthat/test-bmdMA.R
+++ b/tests/testthat/test-bmdMA.R
@@ -9,6 +9,8 @@
 #   - correct bmd estimate (all definitions)
 # - TCDD model (binomial)
 #   - correct bmd estimate (excess + additional)
+# - Lemna model (count)
+#   - correct bmd estimate (all definitions)
 # - S.alba model (continuous with multiple curves)
 #   - correct bmd estimate (point, extra, hybridExc)
 # - Decreasing binomial model with multiple curves
@@ -205,6 +207,7 @@ test_that("bmdMA handles seed correctly when using Stacking weights",{
   
   expect_equal(twoNormalVariables.seed1, c(oneNormalVariable.seed1, secondNormalVariable.seed1))
   expect_equal(bmdMAStackingWeights.seed123.inside$modelWeights, bmdMAStackingWeights.seed123.outside$modelWeights, tolerance = 1e-4)
+  expect_equal(bmdMAStackingWeights.seed123.inside$modelWeights, c(9.23081721161721e-07, 3.65643449053248e-07, 0.33333311279605, 0.66666559847878), tolerance = 1e-6)
   
   # second set of seed (156, 999)
   set.seed(156, kind = "Mersenne-Twister", normal.kind = "Inversion")
@@ -617,6 +620,196 @@ test_that("bmdMA function computes BMD (additional) correctly for TCDD models",
 })
 
 
+# lemna results -----------------------------------------------------------
+
+test_that("bmdMA function computes BMD (point) correctly for lemna models", {
+  # data and fitted models
+  data0 <- drcData::lemna
+  object.LL <- drm(frond.num ~ conc, data = data0, fct = LL.3(), type = "Poisson") 
+  object.LN <- drm(frond.num ~ conc, data = data0, fct = LN.3(), type = "Poisson")  
+  object.W1 <- drm(frond.num ~ conc, data = data0, fct = W1.3(), type = "Poisson") 
+  object.W2 <- drm(frond.num ~ conc, data = data0, fct = W2.3(), type = "Poisson") 
+  modelList0 <- list(object.LL, object.LN, object.W1, object.W2)
+  
+  # results
+  resultKang <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "Kang", display = FALSE)
+  resultBuckland <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "Buckland", display = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBoot <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "bootstrap", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurve <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "curve", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBootBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "bootstrap", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurveBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 52, def = "point", backgType = "modelBased", type = "curve", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  
+  # Expected results based on manual calculation (checked in v2.7.6)
+  # Kang
+  expect_true(!is.na(resultKang$Results[1, "BMD_MA"]))
+  expect_equal(resultKang$Results[1, "BMD_MA"], 4.20666707444774)
+  expect_equal(resultKang$SE[1,1], NA)
+  expect_equal(unname(resultKang$interval[1,]), c(0.735395381212274,7.67793876768321), tolerance = 1e-6)
+  
+  # Buckland
+  expect_true(!is.na(resultBuckland$Results[1, "BMD_MA"]))
+  expect_equal(resultBuckland$Results[1, "BMD_MA"], 4.20666707444774)
+  expect_equal(resultBuckland$SE[1,1], 2.11445972559811, tolerance = 1e-5)
+  expect_equal(unname(resultBuckland$interval[1,]), c(0.728690325754875,7.68464382314061), tolerance = 1e-5)
+  
+  # Boot
+  expect_true(!is.na(resultBoot$Results[1, "BMD_MA"]))
+  expect_equal(resultBoot$Results[1, "BMD_MA"], 4.20666707444774)
+  expect_equal(resultBoot$Boot.samples.used, 50)
+  expect_equal(unname(resultBoot$interval[1,]), c(3.51653291754643,5.06353543927267), tolerance = 1e-4)
+  
+  # Curve
+  expect_true(!is.na(resultCurve$Results[1, "BMD_MA"]))
+  expect_equal(resultCurve$Results[1, "BMD_MA"], 4.19604169207832)
+  expect_equal(resultCurve$Boot.samples.used, 50)
+  expect_equal(unname(resultCurve$interval[1,]), c(3.47683102683814,4.99243797156008), tolerance = 1e-4)
+  
+  # BootBCa
+  expect_true(all(!is.na(resultBootBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultBootBCa$Results[, "BMD_MA"]), c(4.20666707444774))
+  expect_equal(resultBootBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultBootBCa$Results[,"BMDL_MA"]), c(3.56675301982365), tolerance = 1e-4)
+  expect_equal(unname(resultBootBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+  
+  # CurveBCa
+  expect_true(all(!is.na(resultCurveBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultCurveBCa$Results[, "BMD_MA"]), c(4.19604169207832))
+  expect_equal(resultCurveBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultCurveBCa$Results[,"BMDL_MA"]), c(3.52753413584386), tolerance = 1e-4)
+  expect_equal(unname(resultCurveBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+  
+})
+
+test_that("bmdMA function computes BMD (extra) correctly for lemna models", {
+  # data and fitted models
+  data0 <- drcData::lemna
+  object.LL <- drm(frond.num ~ conc, data = data0, fct = LL.3(), type = "Poisson") 
+  object.LN <- drm(frond.num ~ conc, data = data0, fct = LN.3(), type = "Poisson")  
+  object.W1 <- drm(frond.num ~ conc, data = data0, fct = W1.3(), type = "Poisson") 
+  object.W2 <- drm(frond.num ~ conc, data = data0, fct = W2.3(), type = "Poisson") 
+  modelList0 <- list(object.LL, object.LN, object.W1, object.W2)
+  
+  # results
+  resultKang <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "Kang", display = FALSE)
+  resultBuckland <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "Buckland", display = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBoot <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "bootstrap", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurve <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "curve", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBootBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "bootstrap", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurveBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "extra", backgType = "modelBased", type = "curve", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  
+  # Expected results based on manual calculation (checked in v2.6.7)
+  # Kang
+  expect_true(!is.na(resultKang$Results[1, "BMD_MA"]))
+  expect_equal(resultKang$Results[1, "BMD_MA"], 0.853320135656905)
+  expect_equal(resultKang$SE[1,1], NA)
+  expect_equal(unname(resultKang$interval[1,]), c(-1.09581446544656,2.80245473676037), tolerance = 1e-5)
+  
+  # Buckland
+  expect_true(!is.na(resultBuckland$Results[1, "BMD_MA"]))
+  expect_equal(resultBuckland$Results[1, "BMD_MA"], 0.853320135656905)
+  expect_equal(resultBuckland$SE[1,1], 1.2172299131057, tolerance = 1e-5)
+  expect_equal(unname(resultBuckland$interval[1,]), c(-1.14884490174883,2.85548517306264), tolerance = 1e-5)
+  
+  # Boot
+  expect_true(all(!is.na(resultBoot$Results[, "BMD_MA"])))
+  expect_equal(unname(resultBoot$Results[, "BMD_MA"]), c(0.853320135656905))
+  expect_equal(resultBoot$Boot.samples.used, 50)
+  expect_equal(unname(resultBoot$interval[,"BMDL_MA"]), c(0.511946652480417), tolerance = 1e-4)
+  expect_equal(unname(resultBoot$interval[,"BMDU_MA"]), c(1.15648988912451), tolerance = 1e-4)
+  
+  # Curve
+  expect_true(all(!is.na(resultCurve$Results[, "BMD_MA"])))
+  expect_equal(unname(resultCurve$Results[, "BMD_MA"]), c(0.85650746114556))
+  expect_equal(resultCurve$Boot.samples.used, 50)
+  expect_equal(unname(resultCurve$interval[,"BMDL_MA"]), c(0.564987166410153), tolerance = 1e-4)
+  expect_equal(unname(resultCurve$interval[,"BMDU_MA"]), c(1.29046349839692), tolerance = 1e-4)
+  
+  # BootBCa
+  expect_true(all(!is.na(resultBootBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultBootBCa$Results[, "BMD_MA"]), c(0.853320135656905))
+  expect_equal(resultBootBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultBootBCa$Results[,"BMDL_MA"]), c(0.534491473150658), tolerance = 1e-4)
+  expect_equal(unname(resultBootBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+  
+  # CurveBCa
+  expect_true(all(!is.na(resultCurveBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultCurveBCa$Results[, "BMD_MA"]), c(0.85650746114556))
+  expect_equal(resultCurveBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultCurveBCa$Results[,"BMDL_MA"]), c(0.548606564117917), tolerance = 1e-4)
+  expect_equal(unname(resultCurveBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+})
+
+test_that("bmdMA function computes BMD (relative) correctly for lemna models", {
+  # data and fitted models
+  data0 <- drcData::lemna
+  object.LL <- drm(frond.num ~ conc, data = data0, fct = LL.3(), type = "Poisson") 
+  object.LN <- drm(frond.num ~ conc, data = data0, fct = LN.3(), type = "Poisson")  
+  object.W1 <- drm(frond.num ~ conc, data = data0, fct = W1.3(), type = "Poisson") 
+  object.W2 <- drm(frond.num ~ conc, data = data0, fct = W2.3(), type = "Poisson") 
+  modelList0 <- list(object.LL, object.LN, object.W1, object.W2)
+  
+  # results
+  resultKang <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "Kang", display = FALSE)
+  resultBuckland <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "Buckland", display = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBoot <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "bootstrap", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurve <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "curve", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultBootBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "bootstrap", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
+  resultCurveBCa <- bmdMA(modelList0, modelWeights = "AIC", bmr = 0.1, def = "relative", backgType = "modelBased", type = "curve", bootInterval = "BCa", R = 50, display = FALSE, progressInfo = FALSE)
+  
+  # Expected results based on manual calculation (checked in v2.6.7)
+  # Kang
+  expect_true(!is.na(resultKang$Results[1, "BMD_MA"]))
+  expect_equal(resultKang$Results[1, "BMD_MA"], 0.853320135656905)
+  expect_equal(resultKang$SE[1,1], NA)
+  expect_equal(unname(resultKang$interval[1,]), c(-0.923557891551217,2.63019816286503), tolerance = 1e-5)
+  
+  # Buckland
+  expect_true(!is.na(resultBuckland$Results[1, "BMD_MA"]))
+  expect_equal(resultBuckland$Results[1, "BMD_MA"], 0.853320135656905)
+  expect_equal(resultBuckland$SE[1,1], 1.11688376601134, tolerance = 1e-5)
+  expect_equal(unname(resultBuckland$interval[1,]), c(-0.983790177750073,2.69043044906388), tolerance = 1e-5)
+  
+  # Boot
+  expect_true(all(!is.na(resultBoot$Results[, "BMD_MA"])))
+  expect_equal(unname(resultBoot$Results[, "BMD_MA"]), c(0.853320135656905))
+  expect_equal(resultBoot$Boot.samples.used, 50)
+  expect_equal(unname(resultBoot$interval[,"BMDL_MA"]), c(0.511946652480417), tolerance = 1e-4)
+  expect_equal(unname(resultBoot$interval[,"BMDU_MA"]), c(1.15648988912451), tolerance = 1e-4)
+  
+  # Curve
+  expect_true(all(!is.na(resultCurve$Results[, "BMD_MA"])))
+  expect_equal(unname(resultCurve$Results[, "BMD_MA"]), c(0.85650746114556))
+  expect_equal(resultCurve$Boot.samples.used, 50)
+  expect_equal(unname(resultCurve$interval[,"BMDL_MA"]), c(0.564987166410153), tolerance = 1e-4)
+  expect_equal(unname(resultCurve$interval[,"BMDU_MA"]), c(1.29046349839692), tolerance = 1e-4)
+  
+  # BootBCa
+  expect_true(all(!is.na(resultBootBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultBootBCa$Results[, "BMD_MA"]), c(0.853320135656905))
+  expect_equal(resultBootBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultBootBCa$Results[,"BMDL_MA"]), c(0.534491473150658), tolerance = 1e-4)
+  expect_equal(unname(resultBootBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+  
+  # CurveBCa
+  expect_true(all(!is.na(resultCurveBCa$Results[, "BMD_MA"])))
+  expect_equal(unname(resultCurveBCa$Results[, "BMD_MA"]), c(0.85650746114556))
+  expect_equal(resultCurveBCa$Boot.samples.used, 50)
+  expect_equal(unname(resultCurveBCa$Results[,"BMDL_MA"]), c(0.548606564117917), tolerance = 1e-4)
+  expect_equal(unname(resultCurveBCa$interval[,"BMDU_MA"]), c("Not available for BCa bootstrap"))
+})
+
 
 
 # S.alba models -----------------------------------------------------------
@@ -782,6 +975,7 @@ test_that("bmdMA function handles modelWeights argument on S.alba data with mult
   expect_equal(bmdMABICWeights$interval[1,2], sum(bmduVals[1,] * BICWeights0))
   expect_equal(bmdMABICWeights$interval[2,2], sum(bmduVals[2,] * BICWeights0))
   
+  expect_equal(stackingWeights0, c(0.302283264429429, 5.25311121531472e-06, 0.447718064962311, 0.249993417497045), tolerance = 1e-6)
   expect_equal(bmdMAStackingWeights$modelWeights, stackingWeights0)
   expect_equal(bmdMAStackingWeights$Results[1,1], sum(bmdVals[1,] * stackingWeights0), tolerance = 1e-4)
   expect_equal(bmdMAStackingWeights$Results[2,1], sum(bmdVals[2,] * stackingWeights0), tolerance = 1e-4)
@@ -1131,10 +1325,13 @@ test_that("bmdMA function computes BMD (point with stacking weights) correctly f
   expect_equal(resultKang$Boot.samples.used, NA)
   expect_equal(unname(resultKang$interval[,"BMDL_MA"]), c(6.06814221070238,21.8988805204411), tolerance = 1)
   expect_equal(unname(resultKang$interval[,"BMDU_MA"]), c(27.231549080278,38.7342484788404), tolerance = 1)
+  expect_equal(resultKang$modelWeights, c(6.99464226206443e-07, 0.44271627558062, 2.55154979958813e-06, 0.557280473405355), tolerance = 0.05)
   # resultBoot
   expect_true(all(!is.na(resultBoot$Results[, "BMD_MA"])))
   expect_equal(unname(resultBoot$Results[, "BMD_MA"]), c(16.6498456454902, 30.3165644996408), tolerance = 1e-1)
   expect_equal(resultBoot$Boot.samples.used, 50, tolerance = 1)
   expect_equal(unname(resultBoot$interval[,"BMDL_MA"]), c(2.12333205792999,25.08635184492), tolerance = 1)
   expect_equal(unname(resultBoot$interval[,"BMDU_MA"]), c(27.8840479507664,38.6902619548875), tolerance = 1)
+  expect_equal(resultBoot$modelWeights, c(6.99464226206443e-07, 0.44271627558062, 2.55154979958813e-06, 0.557280473405355), tolerance = 0.05)
+  expect_equal(resultKang$modelWeights, resultBoot$modelWeights)
 })
diff --git a/tests/testthat/test-drmMMRE.R b/tests/testthat/test-drmMMRE.R
index d199803..d691e4f 100644
--- a/tests/testthat/test-drmMMRE.R
+++ b/tests/testthat/test-drmMMRE.R
@@ -36,7 +36,7 @@ test_that("Example usage of drmMMRE function", {
                       ncol = 4, nrow = 4,
                       dimnames = list(c("Coefb", "Coefc", "Coefd", "Coefe"),
                                       c("Coefb", "Coefc", "Coefd", "Coefe"))), 
-               tolerance = 1e-6)
+               tolerance = 1e-4)
   
 })