Introduction to the R package survstan

To fit a survival model with the survstan package, the user must choose among one of the available fitting functions *reg(), where * stands for the type of regression model, that is, aft for accelerated failure time (AFT) models, ah for accelerates hazards (AH) models, ph for proportional hazards (PO) models, po for proportional (PO) models, or yp for Yang & Prentice (YP) models.

The specification of the survival formula passed to the chosen *reg() function follows the same syntax adopted in the survival package, so that transition to the survstan package can be smoothly as possible for those familiar with the survival package.

The code below shows the model fitting of a AFT model with Weibull baseline distribution using the survstan::aftreg() function. For comparison purposes, we also fit to the same data the Weibull regression model using survival::survreg() function:

library(survstan)
library(dplyr)

ovarian <- ovarian %>%
  mutate(
    across(c("rx", "resid.ds"), as.factor)
  )

survreg <- survreg(
  Surv(futime, fustat) ~ ecog.ps + rx, 
  dist = "weibull", data = ovarian
)

survstan <- aftreg(
  Surv(futime, fustat) ~ ecog.ps + rx, 
  dist = "weibull", data = ovarian
)

Although the model specification is quite similar, there are some important differences that the user should be aware of. While the model fitted using the survival::survreg() function uses the log scale representation of the AFT model with the presence of an intercept term in the linear predictor, the survstan::aftreg() considers the original time scale for model fitting without the presence of the intercept term in the linear predictor.

To see that, let us summarize the fitted models:

summary(survreg)
#> 
#> Call:
#> survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, 
#>     dist = "weibull")
#>              Value Std. Error     z       p
#> (Intercept)  7.425      0.929  7.99 1.3e-15
#> ecog.ps     -0.385      0.527 -0.73    0.47
#> rx2          0.529      0.529  1.00    0.32
#> Log(scale)  -0.123      0.252 -0.49    0.62
#> 
#> Scale= 0.884 
#> 
#> Weibull distribution
#> Loglik(model)= -97.1   Loglik(intercept only)= -98
#>  Chisq= 1.74 on 2 degrees of freedom, p= 0.42 
#> Number of Newton-Raphson Iterations: 5 
#> n= 26
summary(survstan)
#> Call:
#> aftreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, 
#>     dist = "weibull")
#> 
#> Accelerated failure time model fit with weibull baseline distribution 
#> 
#> Regression coefficients:
#>         Estimate Std. Error z value Pr(>|z|)
#> ecog.ps  -0.3851     0.5270 -0.7307   0.4649
#> rx2       0.5287     0.5292  0.9991   0.3178
#> 
#> Baseline parameters:
#>         Estimate Std. Error       2.5%      97.5%
#> alpha    1.13141    0.28535    0.69014     1.8548
#> gamma 1678.06655 1558.57309  271.78211 10360.9005
#> --- 
#> loglik = -97.08449   AIC = 202.169

Next, we show how to fit a PH model using the survstan::phreg() function. For comparison purposes, the semiparametric Cox model is also fitted to the same data using the function survival::coxph().

survstan <- phreg(
  Surv(futime, fustat) ~ ecog.ps + rx, 
  data = ovarian, dist = "weibull"
)
cox <- coxph(
  Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian
)
coef(survstan)
#>    ecog.ps        rx2 
#>  0.4355449 -0.5981575
coef(cox)
#>    ecog.ps        rx2 
#>  0.3697972 -0.5782271