Parametric models with {icpack}"
In icpack: Survival Analysis of Interval-Censored Data

knitr::opts_chunk$set(fig.width = 6, fig.height = 4,
                      echo = TRUE, warning = FALSE, message = FALSE)

Introduction

The purpose of this vignette is to show how {icpack} can be used to fit certain parametric models, namely with constant hazards, Gompertz models and Weibull models. We will use the ovarian cancer data, because with right-censored data we can easily check the results of {icpack} with other methods. Those other methods do not extend to interval-censored data, while {icpack} does.

For information on the ovarian cancer data we refer to the vignette "Smooth modeling of right-censored survival data with {icpack}", or the help file. First the package needs to be loaded.

library(icpack)
library(ggplot2)

Constant hazards

Using icfit it is possible to fit proportional hazards models to right-censored and interval-censored data with constant baseline hazards. It is of interest to start with a model without covariates and right-censored data, since then we know what the maximum likelihood of the underlying constant hazard is, and we also have a straightforward expression for its confidence interval. Here is how the model is fit using icfit.

ic_fit_cst <- icfit(Surv(time, d) ~ 1, data = Ova,
                    bdeg=1, pord=1, lambda=1e+06, update_lambda=FALSE)
ic_fit_cst

We use bdeg=1 to use B-splines of degree 1 (constants), pord=1 to penalize first order differences of coefficients of the B-splines, lambda=1e+06 to use an extremely high penalty, and update_lambda=FALSE to make sure that lambda is not updated in the process. Note that the effective dimension of the log-hazard is effectively 1, so the log-hazard and hazard are constant. A plot of the hazard shows a straight line.

plot(ic_fit_cst)

The value fo the hazard and its 95% confidence interval can be seen by printing the result of predict.

pic_fit_cst <- predict(ic_fit_cst)[, c("time", "haz", "sehaz", "lowerhaz", "upperhaz")]
head(pic_fit_cst)
tail(pic_fit_cst)

We see that the hazard is estimated as $\hat{\lambda} = 0.000747$, with 95\% confidence interval about 0.000662 -- 0.000842. The estimate and confidence interval are obtained using the general theory outlined in Putter, Gampe \& Eilers (2023). Nevertheless, they coincide closely with standard theory, which says that $\hat{\lambda}{\textrm{MLE}} = D/T$, where $D = \sum{i=1}^n d_i$ is the total number of events and $T = \sum_{i=1}^n t_i$ is the total follow-up time in the data. (Here $(t_i, d_i)$ for subject $i=1,\ldots,n$ are the usual time and status indicators in right-censored time-to-event data.) Moreover, the variance of the log of $\hat{\lambda}{\textrm{MLE}}$ is given by $1/D$, so a 95\% confidence interval for $\lambda$ is given by $\exp( \log(\hat{\lambda}{\textrm{MLE}}) \pm 1.96/\sqrt{D} )$. In our data this would yield

D <- sum(Ova$d)
T <- sum(Ova$time)
rate <- D / T
lograte <- log(rate)
selograte <- 1 / sqrt(D)
lowerlograte <- lograte - qnorm(0.975) * selograte
upperlograte <- lograte + qnorm(0.975) * selograte
data.frame(D=D, T=T, rate=rate, lower=exp(lowerlograte), upper=exp(upperlograte))

We see that this is quite close to what icfit has given us.

We can go one step further with icfit and add the covariate effects to the constant baseline hazard.

ic_fit_cst <- icfit(Surv(time, d) ~ Diameter + FIGO + Karnofsky, data = Ova,
                    bdeg=1, pord=1, lambda=1e+06, update_lambda=FALSE)
ic_fit_cst

Also this model can be fitted in an alternative way, namely using direct Poisson regression. To achieve this d is considered the outcome variable and the logarithm of the follow-up time is added as an offset. Here is the code to run it.

Ova$logtime <- log(Ova$time)
pois_fit <- glm(d ~ Diameter + FIGO + Karnofsky + offset(logtime), family=poisson(), data = Ova)
spois_fit <- summary(pois_fit)
spois_fit

We see that the estimates of the covariate effects are quite similar to the estimates from ic_fit_cst. Moreover the (Intercept) term is the logarithm of the (constant) baseline hazard. If we exponeniate it the corresponding baseline hazard and 95\% confidence interval are given by

lograte <- spois_fit$coefficients[1, 1]
selograte <- spois_fit$coefficients[1, 2]
lowerlograte <- lograte - qnorm(0.975) * selograte
upperlograte <- lograte + qnorm(0.975) * selograte
data.frame(rate=exp(lograte), lower=exp(lowerlograte), upper=exp(upperlograte))

The result from icfit is visualised as

plot(ic_fit_cst, xlab="Days since randomisation", ylab="Hazard",
     title="Constant baseline hazard") + ylim(0, 0.0052)

Naturally, direct Poisson regression using glm is preferable to icfit in case of constant hazards both in terms of speed and accuracy, but it is comforting to see that icfit does give rapid and reliable results also here. Moreover, Poisson regression cannot be used directly with interval-censored data.

Precision in icfit can be increased (at the cost of extra computation time) by increasing nt. If we remove the covariates again, we set nt to 500, and look at the hazard and its 95% confidence interval, which again can be seen by printing the result of predict, we do see that the results are quite a bit closer to that of the occurrence over exposure theoretical ones.

ic_fit_cst <- icfit(Surv(time, d) ~ 1, data = Ova, nt=500,
                    bdeg=1, pord=1, lambda=1e+06, update_lambda=FALSE)
pic_fit_cst <- predict(ic_fit_cst)[, c("time", "haz", "sehaz", "lowerhaz", "upperhaz")]
head(pic_fit_cst)
tail(pic_fit_cst)

Note that only with right-censored data we can use the occurrence over exposure theory and Poisson GLM. Importantly, {icpack} can also be used to fit models (with and without covariates) with constant (baseline) hazards for interval-censored data, while there is no standard alternative for this situation.

"Gompertz" models

When we set the degree of the B-splines to 1 and the order of the penalty to 2, and take a very large value of lambda, we get a model for which the log-hazard is a straight line. The corresponding hazard is of the form $h(t) = a e^{bt}$. This corresponds to the hazard of a Gompertz model, provided $b>0$. Gompertz rates are often used to represent rates of mortality; when $b<0$ then this implies that very large, or indeed infinite values are not uncommon, since then $\int_0^\infty h(t) dt$ is finite. For modeling human mortality this obviously is a problem, but for many other practical purposes we may ignore this issue. But implicitly allowing for negative $b$ by icfit is the reason for referring to "Gompertz" rather than Gompertz models.

Here is the fit of the "Gompertz" model using icfit.

ic_fit_gompertz <- icfit(Surv(time, d) ~ 1, data = Ova,
                         bdeg=1, pord=2, lambda=1e+06, update_lambda=FALSE)
ic_fit_gompertz

We see that the effective dimension of the log-hazard is 2, so a straight line. Below is a plot of the hazard.

plot(ic_fit_gompertz) + ylim(0, 0.00125)

The fact that the log-hazard is a straight line is easier to see after transforming the y-axis.

plot(ic_fit_gompertz) + scale_y_continuous(trans='log10')

Note that the "Gompertz" hazard rate is implicitly obtained by setting the order of the penalty to 2. From the result of icfit the values of $a$ and $b$ and their standard errors can only be derived indirectly.

Weibull models

Also Weibull models can be fitted with icfit in the same indirect way, by setting bdeg, pord, and lambda, with one additional trick, namely time warping. The hazard of a Weibull distribution is of the form $h(t) = \alpha \lambda t^{\alpha-1}$, with $\lambda > 0$. If we take pord=2, lambda = 1e+06, and work with $\log(t)$ instead of $t$ itself, the form of the log-hazard would be $\beta_0 + \beta_1 \log(t)$, which corresponds to a Weibull hazard with $\alpha = \beta_1 + 1$, $\lambda = e^{\beta_0} / (\beta_1 + 1)$.

Here is the fit of the Weibull model with icfit. This works fine with our present way of using time in days; time=1 corresponds to 1 day after randomisation. When time is in years, the logarithm of time would be negative for all events and censorings before one year. In that case adding a constant, or alternatively, defining logtime below as log(Ova$time * 365.25) would solve the issue.

Ova$logtime <- log(Ova$time)
ic_fit_weibull <- icfit(Surv(logtime, d) ~ 1, data = Ova,
                         bdeg=1, pord=2, lambda=1e+06, update_lambda=FALSE)
ic_fit_weibull

Here is a plot of the hazard with log-transformed y-axis. Time on the x-axis actually refers to the logarithm (base 10) of the original time, so 3 on the x-axis corresponds to $t=1000$.

plot(ic_fit_weibull) + scale_y_continuous(trans='log')

To see everything on the original time scale and to facilitate future comparison with a direct parametric fit using survreg, let us extract the data elements of the plot, and explicitly transform time back into the original time.

We need to be careful though. Note that the estimated hazard with transformed time measures the expected number of events per transformed time unit, not in the original time unit. So we need not only back transform time, but also the hazard itself. If $\tilde{h}(u)$ denotes the hazard in transformed time $u = g(t)$, then in the original time scale $t = f(u) = g^{-1}(u)$ we have that $h(t) = \tilde{h}(u) / f^\prime(u)$. With $u = g(t) = \log(t)$, $t = f(u) = \exp(u)$, we get $h(t) = \tilde{h}(u) / \exp(u) = \tilde{h}(\log(t)) / t$.

pw <- plot(ic_fit_weibull)$data
head(pw)[, 1:5]
tail(pw)[, 1:5]
pw$u <- pw$time
pw$time <- exp(pw$u)
pw$haz <- pw$haz / pw$time
pw$lowerhaz <- pw$lowerhaz / pw$time
pw$upperhaz <- pw$upperhaz / pw$time
head(pw)[, 1:5]
tail(pw)[, 1:5]

Here is a plot of the result.

plt <- ggplot(pw) +
  geom_ribbon(aes(x = pw$time, ymin = pw$lowerhaz, ymax = pw$upperhaz),
              fill = 'lightgrey') +
  geom_line(aes(x = pw$time, y = pw$lowerhaz), colour = 'blue') +
  geom_line(aes(x = pw$time, y = pw$upperhaz), colour = 'blue') +
  geom_line(aes(x = pw$time, y = pw$haz), size = 1, colour = 'blue', lty = 1) +
  ylim(0, 0.002) + 
  ggtitle("Weibull hazard fit") + xlab("Days since randomisation") +
  ylab("Hazard") + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
plt

Let us compare this now with survreg. It is most convenient to use the wrapper provided in flexsurvreg, because it uses a parametrisation close to ours. The parametrisation in flexsurvreg is of the form $S(t) = \exp(-(t/\mu^\prime)^{\alpha^\prime})$, compared to our $S(t) = \exp(-\lambda t^{\alpha})$, from which we see that $\alpha = \alpha^\prime$ and $\lambda = 1 / (\mu^\prime)^{\alpha^\prime}$. We add the fit from survreg with a dashed line to the original plot.

library(flexsurv)
fsr <- flexsurvreg(Surv(time, d) ~ 1, data = Ova, dist = "weibull")
fsr

alp <- exp(fsr$coef[1])
mu <- exp(fsr$coef[2])
lam <- 1 / mu^alp
hw <- function(t, alp, lam) alp * lam * t^(alp - 1)
maxt <- max(Ova$time)
pw$haz_survreg <- hw(pw$time, alp, lam)

plt <- plt +
  geom_line(aes(x = pw$time, y = pw$haz_survreg), size = 1, lty = 2)
plt