diff --git a/R/flexsurvreg.R b/R/flexsurvreg.R index 1f56ca1..5ac30b9 100644 --- a/R/flexsurvreg.R +++ b/R/flexsurvreg.R @@ -343,21 +343,24 @@ check.flexsurv.response <- function(Y){ if (!inherits(Y, "Surv")) stop("Response must be a survival object") ### convert Y from Surv object to numeric matrix + type <- attr(Y, "type") ### though "time" only used for initial values, printed time at risk, empirical hazard - if (attr(Y, "type") == "counting") + if (type == "counting") Y <- cbind(Y, time=Y[,"stop"] - Y[,"start"], time1=Y[,"stop"], time2=Inf) - else if (attr(Y, "type") == "interval"){ + else if (type == "interval"){ # Surv() converts interval2 to interval Y[,"time2"][Y[,"status"]==0] <- Inf # upper bound with right censoring + Y[,"time2"][Y[,"status"]==1] <- Y[,"time1"][Y[,"status"]==1] # event time Y[,"time2"][Y[,"status"]==2] <- Y[,"time1"][Y[,"status"]==2] Y[,"time1"][Y[,"status"]==2] <- 0 # Y <- cbind(Y, start=0, stop=Y[,"time1"], time=Y[,"time1"]) } - else if (attr(Y, "type") == "right") + else if (type == "right") Y <- cbind(Y, start=0, stop=Y[,"time"], time1=Y[,"time"], time2=Inf) else stop("Survival object type \"", attr(Y, "type"), "\"", " not supported") if (any(Y[,"time1"]<0)){ stop("Negative survival times in the data") } + attr(Y, "type") <- type Y } diff --git a/R/plot.flexsurvreg.R b/R/plot.flexsurvreg.R index 5ede4cd..b51a32b 100644 --- a/R/plot.flexsurvreg.R +++ b/R/plot.flexsurvreg.R @@ -96,8 +96,8 @@ plot.flexsurvreg <- function(x, newdata=NULL, X=NULL, type="survival", fn=NULL, add=FALSE,...) { ## don't calculate or plot CIs by default if all covs are categorical -> multiple curves - if (!add && is.null(x[["data"]])) - stop("Goodness-of-fit plots are not available if the data have been removed from the model object") + if (!add && is.null(x[["data"]])) + stop("Goodness-of-fit plots are not available if the data have been removed from the model object") if (is.null(ci)) ci <- ((x$ncovs == 0) || (!all(x$covdata$isfac))) if (!ci) B <- 0 @@ -112,7 +112,10 @@ plot.flexsurvreg <- function(x, newdata=NULL, X=NULL, type="survival", fn=NULL, mf <- model.frame(x) Xraw <- mf[,attr(mf, "covnames.orig"), drop=FALSE] mm <- as.data.frame(model.matrix(x)) - form <- "Surv(dat$Y[,\"start\"],dat$Y[,\"stop\"],dat$Y[,\"status\"]) ~ " + if (isTRUE(attr(dat$Y, "type")=="interval")){ # guard against NULL from older versions + form <- "Surv(dat$Y[,\"time1\"],dat$Y[,\"time2\"],type=\"interval2\") ~ " + } else + form <- "Surv(dat$Y[,\"start\"],dat$Y[,\"stop\"],dat$Y[,\"status\"]) ~ " form <- paste(form, if (x$ncovs > 0 && all(x$covdata$isfac)) paste("mm[,",1:x$ncoveffs,"]", collapse=" + ") else 1) form <- as.formula(form) ## If any continuous covariates, it is hard to define subgroups @@ -139,7 +142,7 @@ plot.flexsurvreg <- function(x, newdata=NULL, X=NULL, type="survival", fn=NULL, group <- if(x$ncovs>0) do.call("interaction", mm) else factor(rep(0,nrow(dat$Y))) Xgroup <- factor(do.call("interaction", as.data.frame(X)), levels=levels(group)) haz <- list() - for (i in 1:nrow(X)) { + for (i in 1:nrow(X)) { mha <- muhaz.args if (is.null(mha$max.time)) mha$max.time <- with(as.data.frame(dat$Y[group==Xgroup[i],,drop=FALSE]), max(time[status==1])) haz[[i]] <- do.call("muhaz", c(list(times=dat$Y[,"time"], delta=dat$Y[,"status"], subset=(group==Xgroup[i])), mha)) diff --git a/tests/testthat/test_flexsurvreg.R b/tests/testthat/test_flexsurvreg.R index d71a837..88a0206 100644 --- a/tests/testthat/test_flexsurvreg.R +++ b/tests/testthat/test_flexsurvreg.R @@ -367,10 +367,12 @@ test_that("Interval censoring",{ sr2 <- survreg(Surv(tmin, status) ~ 1, dist="weibull") expect_equal(sr1$loglik[2], sr2$loglik[2]) - fs1 <- flexsurvreg(Surv(tmin, tmax.sr, type="interval2") ~ 1, dist="weibull") + fs1_inf <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1, dist="weibull") + fs1_na <- flexsurvreg(Surv(tmin, tmax.sr, type="interval2") ~ 1, dist="weibull") fs2 <- flexsurvreg(Surv(simt, status) ~ 1, dist="weibull") - expect_equal(fs1$loglik, fs2$loglik) - expect_equal(fs1$loglik, sr1$loglik[2]) + expect_equal(fs1_inf$loglik, fs2$loglik) + expect_equal(fs1_na$loglik, fs2$loglik) + expect_equal(fs1_inf$loglik, sr1$loglik[2]) ## put an upper bound on censored times tmax <- ifelse(status==1, simt, 0.7) @@ -382,7 +384,14 @@ test_that("Interval censoring",{ status[status==0] <- 3 fs3 <- flexsurvreg(Surv(tmin, tmax, status, type="interval") ~ 1, dist="weibull") expect_equal(fs1$loglik, fs3$loglik) - + + ## interval censoring with zero width interval + set.seed(1) + tmin <- tmax <- rweibull(100, 1.1, 1.5) + fs1 <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1, dist="weibull") + fs2 <- flexsurvreg(Surv(tmin) ~ 1, dist="weibull") + expect_equal(fs1$loglik, fs2$loglik) + ## interval censoring close around the event set.seed(1) tev <- rweibull(100, 1.1, 1.5) @@ -391,7 +400,7 @@ test_that("Interval censoring",{ fs1 <- flexsurvreg(Surv(tev) ~ 1, dist="weibull") fs2 <- flexsurvreg(Surv(tmin, tmax, type="interval2") ~ 1, dist="weibull") expect_equal(fs2$res["shape","est"], fs1$res["shape","est"], tol=1e-03) - + ## 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))