-
Notifications
You must be signed in to change notification settings - Fork 3
/
README.Rmd
139 lines (111 loc) · 6.05 KB
/
README.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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
---
output:
md_document:
variant: markdown_github
bibliography: README.bib
csl: vignettes/bib_style.csl
---
```{r setup, echo = FALSE, cache=FALSE}
knitr::opts_chunk$set(
collapse = TRUE, comment = "#>", fig.path = "README-", dpi = 124,
error = FALSE)
```
# dynamichazard
[](https://github.com/boennecd/dynamichazard/actions)
[](https://www.r-pkg.org/badges/version/dynamichazard)
[](http://cran.rstudio.com/web/packages/dynamichazard/index.html)
The goal of dynamichazard is to estimate time-varying effects in survival analysis. The time-varying effects are estimated with state space models where the coefficients follow a given order random walk. The advantageous of using state space models is that you can extrapolate beyond the last observed time period. For more details, see the ddhazard vignette at https://cran.r-project.org/web/packages/dynamichazard/vignettes/ddhazard.pdf.
The particle filter and smoother methods can estimate more general models then the random walk model. See the [/examples](/examples) directory for some examples.
## Installation
You can install dynamichazard from github with:
```R
# install.packages("remotes")
remotes::install_github("boennecd/dynamichazard")
```
You can also download the package from CRAN by calling:
```R
install.packages("dynamichazard")
```
## Example - ddhazard
I will use the `aids` data set from the `JMbayes` package. The data set is from a a randomized clinical trial between two drugs for HIV or aids patients. The event is when the patient die. Patient are right censored at the end of the study. The data set is longitudinal/panel format with rows for patient. Though, the only longitudinal variable is the `CD4` count (T-cells count) which is presumably affected by the drug. Thus, I will not use it in the model. The other the other columns of interest are:
* `AZT` is one of the two enrollment criteria. It indicates whether the patient was enrolled due to intolerance to the drug zidovudine or whether the drug failed prior to the study start.
* `prevOI` is the other enrollment criteria. Patients are enrolled either due AIDS diagnosis or two CD4 counts of 300 or fewer. The variable indicates which is the case.
* `drug` is either `ddC` or `ddI` depending on which of the two drugs the patient is randomly assigned to.
* `gender`.
The analysis is given below with comments:
```{r ddhazard_fit, cache = TRUE}
library(dynamichazard)
library(JMbayes) # Contain the aids data set
# We remove the data we dont neeed
aids <- aids[aids$Time == aids$stop, ]
aids <- aids[, !colnames(aids) %in% c("Time", "death", "obstime", "CD4")]
# A look at head of data
head(aids)
max(aids$stop) # Last observation time
max(aids$stop[aids$event == 1]) # Last person with event
# Fit model with extended Kalman filter
fit <- ddhazard(
Surv(stop, event) ~ ddFixed_intercept() + ddFixed(AZT) + gender +
ddFixed(drug) + ddFixed(prevOI),
aids,
model = "exponential", # piecewise constant exponentially distributed
# arrivals times
by = .5, # Length of time intervals in state space
# model
max_T = 19, # Last period we observe when modeling
Q = .1^2, # Covariance matrix for state equation in
# first iteration
Q_0 = 1, # Covariance matrix for the prior
control = ddhazard_control(
eps = .001, # tolerance for EM-algorithm
LR = .9, # Learning rate
n_max = 100 # Max number iterations in EM
))
# Plot the estimates. Dashed lines are 95% confidence bounds
plot(fit)
# Bootstrap the estimates
set.seed(87754771)
boot_out <- ddhazard_boot(fit, R = 1000) # R is number of bootstrap samples
# Plot bootstrap estimates. Dashed lines are 2.5% and 97.5% quantiles of the
# bootstrap estimates. Transparent lines are bootstrap estimates
plot(fit, ddhazard_boot = boot_out)
```
Bootstrapping only slightly changes the confidence bounds. An example of a paper analyzing the CD4 count can be found in [@guo2004]. They also fit a static model (time-invariant coefficients) of the survival times with an exponential model. The estimates are comparable with those above as expected.
## Example - particle filter and smoother
A particle filter and smoother is also included in the package. The computational complexity of these methods match those of the extended Kalman filter but with a much larger constant. Below, I fit a model for the `aids` data. We only use a time-varying effect for gender this time.
<!--
knitr::opts_knit$set(output.dir = ".")
knitr::load_cache("pf_fit", path = "README_cache/markdown_github/")
-->
```{r pf_fit, cache = 1}
set.seed(20170907)
pf_fit <- PF_EM(
Surv(stop, event) ~ ddFixed_intercept() + ddFixed(AZT) + gender +
ddFixed(drug) + ddFixed(prevOI),
aids, model = "exponential",
by = .5, max_T = 19, Q = .5^2, Q_0 = 1,
control = PF_control(
n_threads = 4, # I use a quad-core machine
# set number of particles
N_fw_n_bw = 1000, N_first = 1000,
N_smooth = 1, # Does not matter with Brier_O_N_square
smoother = "Brier_O_N_square", # Select smoother
eps = 1e-5, n_max = 1000, est_a_0 = FALSE, nu = 6L,
covar_fac = 1.2, averaging_start = 150L)
#, trace = 1 # comment back to get feedback during estimation
)
```
```{r pf_plots}
# Compare estimates of Q
pf_fit$Q / .5
fit$Q
plot(pf_fit)
plot(pf_fit$log_likes) # log-likelihoods
logLik(pf_fit)
```
```{r pf_better_log_like, cache = 1, dependson = "pf_fit"}
# better estimate of final log-likelihood
logLik(PF_forward_filter(pf_fit, N_fw = 10000, N_first = 10000))
```
For more details, see the "Particle filters in the dynamichazard package" vignette at https://cran.r-project.org/web/packages/dynamichazard/vignettes/Particle_filtering.pdf
# References