Bug fix for bhazard with no observed event times

chjackson · Nov 22, 2023 · f98312c · f98312c
1 parent ea20567
commit f98312c
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 17 deletions.
diff --git a/R/deriv.R b/R/deriv.R
@@ -68,7 +68,7 @@ Dminusloglik.flexsurv <- function(optpars, Y, X=0, weights, bhazard, rtrunc, dli
 }
 
 dderiv <- function(ddfn, ddcall, X, mx, dlist){
-    if (length(ddcall$t) == 0) matrix(0) else { 
+    if (length(ddcall$t) == 0) array(dim=c(0,length(dlist$pars))) else { 
         res.base <- do.call(ddfn, ddcall)
         res.beta <- Dcov(res.base, X, mx, dlist)
         cbind(res.base, res.beta)

diff --git a/R/deriv2.R b/R/deriv2.R
@@ -268,21 +268,23 @@ D2minusloglik.flexsurv <- function(optpars, Y, X=0, weights, bhazard, rtrunc, dl
         d1LS <- dderiv(dfns$DLS, ddcall, X[dead,,drop=FALSE], mx, dlist)
         d2LS <- d2deriv(dfns$D2LS, ddcall, X[dead,,drop=FALSE], mx, dlist)
         dhazinv <- surv/dens*(d1LS - d1Ld)
-
-        ## given two matrices with dims (n,p), 
-        ## return array of dim(n,p,p) with outer product of each pair of matrix rows
-        vouter <- function(y,z){
-          rows <- seq_len(nrow(y))
-          res <- mapply(outer, split(y, rows), split(z, rows), SIMPLIFY=FALSE)
-          res <- simplify2array(res, except = NULL)
-          aperm(res, c(3,1,2))
-        }
-
-        doff <- - offseti*bw*(offseti*vouter(dhazinv, dhazinv)*bw +
-                              (d2Ld - d2LS)*surv/dens +
-                              vouter(d1Ld - d1LS, dhazinv))
 
-        res <- res - apply(doff*weights[dead],2:3,sum)
+        if (any(dead)){
+            ## given two matrices with dims (n,p), 
+            ## return array of dim(n,p,p) with outer product of each pair of matrix rows
+            vouter <- function(y,z){
+                rows <- seq_len(nrow(y))
+                res <- mapply(outer, split(y, rows), split(z, rows), SIMPLIFY=FALSE)
+                res <- simplify2array(res, except = NULL)
+                aperm(res, c(3,1,2))
+            }
+
+            doff <- - offseti*bw*(offseti*vouter(dhazinv, dhazinv)*bw +
+                                      (d2Ld - d2LS)*surv/dens +
+                                      vouter(d1Ld - d1LS, dhazinv))
+
+            res <- res - apply(doff*weights[dead],2:3,sum)
+        }
     }
     ## currently wastefully calculates derivs for fixed pars then discards them
     optpars <- setdiff(1:npars, fixedpars)

diff --git a/R/flexsurvreg.R b/R/flexsurvreg.R
@@ -51,7 +51,7 @@ logLikFactory <- function(Y, X=0, weights, bhazard, rtrunc, dlist,
     nbpars <- length(dlist$pars)
     insert.locations <- setdiff(seq_len(npars), fixedpars)
 
-    ## which are the subjects with events
+    ## which are the subjects with known event times
     event <- Y[,"status"] == 1
     event.times <- Y[event, "time1"]
     left.censor <- Y[!event, "time2"]

diff --git a/tests/testthat/test_flexsurvreg.R b/tests/testthat/test_flexsurvreg.R
@@ -477,8 +477,12 @@ test_that("No events in the data",{
     set.seed(1)
     tmin <- rexp(100, 1) 
     tmax <- tmin + 0.1
+    bhazard <- rep(0.0001, 100)
     mod <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1, dist="exponential")
     expect_equal(mod$loglik, -337.9815, tolerance=1e-03)
+    modb <- flexsurvreg(Surv(tmin,tmax,type='interval2')~1,
+                       bhazard=bhazard,dist="exponential")
+    expect_equal(mod$loglik, modb$loglik, tolerance=1e-03)
 })
 
 test_that("No censoring in the data",{
@@ -520,4 +524,4 @@ test_that("With and without analytic Hessian", {
   fln <- flexsurvreg(formula = Surv(ovarian$futime, ovarian$fustat) ~ 1, dist="gengamma",
                      hess.control = list(numeric=TRUE))
   expect_equivalent(fla$cov, fln$cov) # analytic derivatives not available
-})
+})