diff --git a/R/llcont.R b/R/llcont.R index d8e496a..7463431 100644 --- a/R/llcont.R +++ b/R/llcont.R @@ -53,12 +53,16 @@ llcont.glm <- function(x, ...){ if(is.matrix(y)) { ## Bolt: replaced apply(..., 1, sum) with optimized rowSums() for performance n <- rowSums(y) - y <- ifelse(n == 0, 0, y[, 1]/n) + ## Bolt: replaced ifelse with preallocation and vectorized subsetting for performance + y <- y[, 1]/n + y[n == 0] <- 0 } else { n <- rep.int(1, length(y)) } m <- if (any(n > 1)) n else wt - wt <- ifelse(m > 0, (wt/m), 0) + ## Bolt: replaced ifelse with preallocation and vectorized subsetting for performance + wt <- wt / m + wt[m <= 0] <- 0 dbinom(round(m * y), round(m), mpreds, log = TRUE) * wt }, quasibinomial = { diff --git a/tests/testthat/test_llcont_glm_binomial_zero_contract.R b/tests/testthat/test_llcont_glm_binomial_zero_contract.R new file mode 100644 index 0000000..d86af8d --- /dev/null +++ b/tests/testthat/test_llcont_glm_binomial_zero_contract.R @@ -0,0 +1,20 @@ +test_that("glm binomial zero totals and zero prior weights preserve likelihood contributions", { + observed <- data.frame( + success = c(0, 1, 2, 3, 1, 4, 2), + failure = c(0, 4, 3, 2, 4, 1, 3), + x = 0:6, + prior_weight = c(1, 1, 1, 0, 1, 1, 1) + ) + + fit <- glm( + cbind(success, failure) ~ x, + data = observed, + family = binomial(), + weights = prior_weight + ) + contributions <- llcont(fit) + + expect_true(all(is.finite(contributions))) + expect_equal(contributions[c(1, 4)], c(0, 0)) + expect_equal(sum(contributions), as.numeric(logLik(fit))) +})