diff --git a/.jules/bolt.md b/.jules/bolt.md
index f658475..970bc0f 100644
--- a/.jules/bolt.md
+++ b/.jules/bolt.md
@@ -15,3 +15,6 @@
 ## 2024-05-15 - [R Performance: ifelse Overhead]
 **Learning:** In R, ifelse evaluates both true and false branches entirely before subsetting, which is very inefficient for vector operations.
 **Action:** Optimize this by preallocating with res <- Y * 0 to preserve attributes and using vectorized subsetting like if any cond res subset <- ...
+## 2024-05-15 - R Performance: ifelse Overhead in nested subsets
+**Learning:** In R, `ifelse()` evaluates both true and false branches entirely before subsetting, which is inefficient. By preallocating with `res <- Y * 0` (to preserve attributes) and using vectorized subsetting, we eliminate this overhead.
+**Action:** Replace `ifelse()` with preallocation that preserves attributes (e.g., `res <- Y * 0` instead of `numeric(n)`) and vectorized subsetting (e.g., `if (any(cond)) res[cond] <- ...`) to optimize performance in frequently evaluated loops or operations.
diff --git a/R/llcont.R b/R/llcont.R
index d8e496a..9d77ec5 100644
--- a/R/llcont.R
+++ b/R/llcont.R
@@ -53,12 +53,22 @@ 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_res <- y[, 1] * 0
+               cond_n <- n > 0
+               cond_n[is.na(cond_n)] <- FALSE
+               if (any(cond_n)) y_res[cond_n] <- (y[, 1]/n)[cond_n]
+               y <- y_res
              } 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_res <- wt * 0
+             cond_m <- m > 0
+             cond_m[is.na(cond_m)] <- FALSE
+             if (any(cond_m)) wt_res[cond_m] <- (wt/m)[cond_m]
+             wt <- wt_res
              dbinom(round(m * y), round(m), mpreds, log = TRUE) * wt
            },
            quasibinomial = {