/
survstan.Rmd
107 lines (78 loc) · 3.67 KB
/
survstan.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
---
title: "Introduction to the R package survstan"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Introduction to the R package survstan}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
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, `yp` for Yang \& Prentice (YP) models, or `eh` for extended hazard (EH) 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 an 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:
```{r setup, message=FALSE}
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:
```{r}
summary(survreg)
summary(survstan)
```
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()`.
```{r}
survstan <- phreg(
Surv(futime, fustat) ~ ecog.ps + rx,
data = ovarian, dist = "weibull"
)
cox <- coxph(
Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian
)
coef(survstan)
coef(cox)
```
<!-- Now we and the semiparametric proportional hazards model using the `survival::coxph()` -->
<!-- The usege of other functions implemented in the __survstan__ package are illustrated in the code presented below: -->
<!-- ```{r} -->
<!-- # fitting the model: -->
<!-- fit <- aftreg( -->
<!-- Surv(futime, fustat) ~ ., -->
<!-- dist = "weibull", data = ovarian -->
<!-- ) -->
<!-- # investigating the fitted model: -->
<!-- estimates(fit) -->
<!-- coef(fit) -->
<!-- confint(fit) -->
<!-- summary(fit) -->
<!-- tidy(fit) -->
<!-- vcov(fit) -->
<!-- # residual plots: -->
<!-- ggresiduals(fit, type = "coxsnell") -->
<!-- ggresiduals(fit, type = "martingale") -->
<!-- ggresiduals(fit, type = "deviance") -->
<!-- # Deviance analysis: -->
<!-- fit1 <- aftreg(Surv(futime, fustat) ~ 1, baseline = "weibull", data = ovarian, init = 0) -->
<!-- fit2 <- aftreg(Surv(futime, fustat) ~ rx, baseline = "weibull", data = ovarian, init = 0) -->
<!-- fit3 <- aftreg(Surv(futime, fustat) ~ rx + ecog.ps, baseline = "weibull", data = ovarian, init = 0) -->
<!-- anova(fit1, fit2, fit3) -->
<!-- # Model comparison using AIC: -->
<!-- AIC(fit1, fit2, fit3) -->
<!-- ``` -->