Skip to content

Commit

Permalink
Changed bhazard to take death rather than survival probability for in…
Browse files Browse the repository at this point in the history
…terval censoring
  • Loading branch information
chjackson committed Dec 1, 2023
1 parent 2f43e46 commit 99a3b85
Show file tree
Hide file tree
Showing 3 changed files with 8 additions and 11 deletions.
9 changes: 3 additions & 6 deletions R/flexsurvreg.R
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ logLikFactory <- function(Y, X=0, weights, bhazard, rtrunc, dlist,
haz_excess <- exp(loghaz_excess)
logdens_offset <- log(1 + bhazard[event] / haz_excess) # = log(haz_allcause / haz_excess)
if (!all(event)) { # background survival S* and left or interval censoring
b_condsurv <- bhazard[!event] # this is S*(end) / S*(start)
b_condsurv <- 1 - bhazard[!event] # this is S*(end) / S*(start)
b_condsurv[lcens.times==Inf] <- 0 # when end=Inf, i.e. right censoring
}
if (any(is.finite(rtrunc)))
Expand Down Expand Up @@ -464,7 +464,7 @@ compress.model.matrices <- function(mml){
##'
##' \item For people whose event time is
##' left-censored or interval-censored, \code{bhazard} should contain the
##' probability of surviving until the end of the corresponding interval,
##' probability of dying by the end of the corresponding interval,
##' conditionally on being alive at the start.
##'
##' \item For people whose event time
Expand Down Expand Up @@ -856,10 +856,7 @@ flexsurvreg <- function(formula, anc=NULL, data, weights, bhazard, rtrunc, subse
## a) calling model.frame() directly doesn't work. it only looks in
## "data" or the environment of "formula" for extra variables like
## "weights". needs to also look in environment of flexsurvreg.
## m <- model.frame(formula=, data=data, weights=weights, subset=subset, na.action=na.action)
## b) putting block below in a function doesn't work when calling
## flexsurvreg within a function
## m <- make.model.frame(call, formula, data, weights, subset, na.action, ancnames)
## should handle calling flexsurvreg within a function

## Make model frame
indx <- match(c("formula", "data", "weights", "bhazard", "rtrunc", "subset", "na.action"), names(call), nomatch = 0)
Expand Down
2 changes: 1 addition & 1 deletion man/flexsurvreg.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

8 changes: 4 additions & 4 deletions tests/testthat/test_flexsurvreg.R
Original file line number Diff line number Diff line change
Expand Up @@ -404,15 +404,15 @@ test_that("Interval censoring",{
## relative survival with left and interval censoring
## at left cens times, bhazard contains the background cond prob of surviving interval
bh <- rep(0.01, length(tmax))
back_csurv <- rep(exp(-0.002*0.01), length(tmax))
back_cdeath <- 1 - rep(exp(-0.002*0.01), length(tmax))
fs1 <- flexsurvreg(Surv(tev) ~ 1, dist="weibull", bhazard=bh)
fs2 <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1,
dist="weibull", bhazard=back_csurv)
dist="weibull", bhazard=back_cdeath)

fs1 <- flexsurvreg(Surv(tev) ~ 1,
dist="weibull", bhazard=bh)
fs2 <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1,
dist="weibull", bhazard=back_csurv, inits=c(1.24, 1.4))
dist="weibull", bhazard=back_cdeath, inits=c(1.24, 1.4))
expect_equal(fs2$res["shape","est"], fs1$res["shape","est"], tol=1e-03)
})

Expand Down Expand Up @@ -508,7 +508,7 @@ test_that("No events in the data",{
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=exp(-bhazard*(tmax-tmin)),dist="exponential")
bhazard = 1 - exp(-bhazard*(tmax-tmin)), dist="exponential")
expect_equal(mod$res["rate","est"], modb$res["rate","est"], tolerance=1e-02)
})

Expand Down

0 comments on commit 99a3b85

Please sign in to comment.