To fit a survival model in 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 StdErr z.value p.value
#> ecog.ps -0.3851 0.5270 -0.7307 0.4649
#> rx2 0.5287 0.5292 0.9991 0.3178
#>
#> Baseline parameters:
#> estimate se 2.5% 97.5%
#> alpha 1.13141 0.25221 0.63709 1.6257
#> gamma 1678.06655 1139.62654 -555.56042 3911.6935
#> ---
#> loglik = -97.08449 AIC = 202.169Next, 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().