Merge pull request #72 from OHDSI/infinite_profile_fix

Fix for edge case when likelihood is infinite when profiling
OHDSI · Nov 9, 2023 · 862ff7e · 862ff7e
2 parents a2b3a95 + 6527620
commit 862ff7e
Show file tree

Hide file tree

Showing 2 changed files with 35 additions and 6 deletions.
diff --git a/R/ModelFit.R b/R/ModelFit.R
@@ -937,21 +937,22 @@ getCyclopsProfileLogLikelihood <- function(object,
         while (length(grid) != 0) {
             ll <- fixedGridProfileLogLikelihood(object, parm, grid, includePenalty)
             profile <- bind_rows(profile, ll) %>% arrange(.data$point)
-
-            if (any(is.nan(profile$value))) {
-                if (all(is.nan(profile$value))) {
+            invalid <- is.nan(profile$value) | is.infinite(profile$value)
+            if (any(invalid)) {
+                if (all(invalid)) {
                     warning("Failing to compute likelihood at entire initial grid.")
                     return(NULL)
                 }
 
-                start <- min(which(!is.nan(profile$value)))
-                end <- max(which(!is.nan(profile$value)))
+                start <- min(which(!invalid))
+                end <- max(which(!invalid))
                 if (start == end) {
                     warning("Failing to compute likelihood at entire grid except one. Giving up")
                     return(NULL)
                 }
                 profile <- profile[start:end, ]
-                if (any(is.nan(profile$value))) {
+                invalid <- invalid[start:end]
+                if (any(invalid)) {
                     warning("Failing to compute likelihood in non-extreme regions. Giving up.")
                     return(NULL)
                 }

diff --git a/tests/testthat/test-profileLikelihood.R b/tests/testthat/test-profileLikelihood.R
@@ -97,3 +97,31 @@ start, length, event, x1, x2
     expect_equal(coef(fit)["x1"], argMax)
 })
 
+test_that("Check adapative profiling likelihood, other covariate is perfect predictor", {
+
+    test <- read.table(header=T, sep = ",", text = "
+start, length, event, x1, x2
+0, 4,  1,0,1
+0, 4,1,2,1
+0, 4,  0,0,0
+0, 4,1,0,1
+0, 4,  1,1,1
+0, 4,0,1,0
+0, 4,  1,1,1
+")
+    data <- createCyclopsData(event ~ x1 + x2, data = test,
+                              modelType = "cpr")
+
+    # Recreating combination of circumstances that lead to infinite log likelihood
+    # Note: infinite variance would have been selected by cross-validation if we had enough data:
+    fit <- fitCyclopsModel(data,
+                           prior = createPrior("laplace", variance = Inf, exclude = "x1"))
+    expect_warning(getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0, 2), initialGridSize = 10),
+                   "Failing to compute likelihood at entire initial grid")
+
+    # This just produces NaN likelihood:
+    fit <- fitCyclopsModel(data)
+    expect_warning(getCyclopsProfileLogLikelihood(fit, "x1", bounds = c(0, 2), initialGridSize = 10),
+                   "Failing to compute likelihood at entire initial grid")
+})
+