diff --git a/R/cfr_rolling.R b/R/cfr_rolling.R
index 6491c258..e67b848c 100644
--- a/R/cfr_rolling.R
+++ b/R/cfr_rolling.R
@@ -121,6 +121,9 @@ cfr_rolling <- function(data,
     # guaranteed in the early stages of an outbreak due to poor data
     p_mid_values <- cumulative_deaths / round(cumulative_outcomes)
 
+    # Handle edge cases where p_mid might be NA or Inf
+    p_mid_values[is.na(p_mid_values) | is.infinite(p_mid_values)] <- 0
+
     if (any(is.infinite(p_mid_values) | p_mid_values < 1e-4)) {
       message(
         "Some daily ratios of total deaths to total cases with known outcome",
diff --git a/R/estimate_severity.R b/R/estimate_severity.R
index f2bfdbfd..b310876e 100644
--- a/R/estimate_severity.R
+++ b/R/estimate_severity.R
@@ -133,17 +133,16 @@
 .select_func_likelihood <- function(total_cases, poisson_threshold, p_mid) {
   # NOTE: internal function is not input checked
   # switch likelihood function based on total cases and p_mid
-  # Binomial approx
-  if (total_cases < poisson_threshold || (p_mid >= 0.05)) {
-    func_likelihood <- function(total_outcomes, total_deaths, pp) {
-      lchoose(round(total_outcomes), total_deaths) +
-        (total_deaths * log(pp)) +
-        (((total_outcomes) - total_deaths) * log(1.0 - pp))
-    }
+
+  # Default to binomial likelihood
+  func_likelihood <- function(total_outcomes, total_deaths, pp) {
+    lchoose(round(total_outcomes), total_deaths) +
+      (total_deaths * log(pp)) +
+      (((total_outcomes) - total_deaths) * log(1.0 - pp))
   }
 
-  # Poisson approx
-  if ((total_cases >= poisson_threshold) && (p_mid < 0.05)) {
+  # Poisson approx - only switch to Poisson if conditions are met
+  if ((total_cases >= poisson_threshold) && (!is.na(p_mid)) && (p_mid < 0.05)) {
     func_likelihood <- function(total_outcomes, total_deaths, pp) {
       stats::dpois(
         total_deaths, pp * round(total_outcomes),