From 99a3b85b6e74715cce5f24be1e53d8b33a6217a8 Mon Sep 17 00:00:00 2001 From: Chris Jackson Date: Fri, 1 Dec 2023 20:46:16 +0000 Subject: [PATCH] Changed bhazard to take death rather than survival probability for interval censoring --- R/flexsurvreg.R | 9 +++------ man/flexsurvreg.Rd | 2 +- tests/testthat/test_flexsurvreg.R | 8 ++++---- 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/R/flexsurvreg.R b/R/flexsurvreg.R index 5ac30b9..f6e6088 100644 --- a/R/flexsurvreg.R +++ b/R/flexsurvreg.R @@ -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))) @@ -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 @@ -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) diff --git a/man/flexsurvreg.Rd b/man/flexsurvreg.Rd index 9cbad84..8bda33e 100644 --- a/man/flexsurvreg.Rd +++ b/man/flexsurvreg.Rd @@ -102,7 +102,7 @@ maximising the weighted sum of the individual-specific log-likelihoods.} \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 diff --git a/tests/testthat/test_flexsurvreg.R b/tests/testthat/test_flexsurvreg.R index 88a0206..35c50de 100644 --- a/tests/testthat/test_flexsurvreg.R +++ b/tests/testthat/test_flexsurvreg.R @@ -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) }) @@ -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) })