vignettes/Time_Varying_Effects.Rmd
Time_Varying_Effects.Rmd
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
-CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)Our aim is to assess the strength of the association between the risk
of the composite event and the level of serum bilirubin. We will
describe the patient-specific profiles over time for this biomarker
using a linear mixed-effects model, where both in the fixed- and
random-effects we include an intercept, the linear and quadratic time
effect. In the fixed-effects we also include the interaction of the time
-effect and sex. The syntax to fit this model with lme()
+effect and sex. The syntax to fit this model with lme()
is:
-fm <- lme(log(serBilir) ~ poly(year, 2) * sex, data = pbc2,
- random = ~ poly(year, 2) | id, control = lmeControl(opt = 'optim'))fm <- lme(log(serBilir) ~ poly(year, 2) * sex, data = pbc2,
+ random = ~ poly(year, 2) | id, control = lmeControl(opt = 'optim'))
The default call to jm() adds the subject-specific
linear predictor of the mixed model as a time-varying covariate in the
survival relative risk model:
To specify that the association of serum bilirubin may change over
time we include an interaction of this time-varying covariate with a
-natural cubic spline of time using function ns() from the
+natural cubic spline of time using function ns() from the
splines package. Important Note: for
-this to work correctly we need to explicitly specify the knots and
+this to work correctly we need to explicitly specify the internal and
boundary knots for the B-splines basis, i.e., in the following example
we set the internal knots at 3, 6 and 9 years, and the boundary knots at
0 and 14.5 years:
-form_splines <- ~ value(log(serBilir)) * ns(year, k = c(3, 6, 9), B = c(0, 14.5))
+form_splines <- ~ value(log(serBilir)) * ns(year, k = c(3, 6, 9), B = c(0, 14.5))
jointFit2 <- update(jointFit1, functional_forms = form_splines,
n_iter = 6500L, n_burnin = 2500L)
summary(jointFit2)
@@ -289,11 +289,12 @@ Non Proportional Hazards#> thinning: 1
#> time: 33 sec
The spline coefficients do not have a straightforward interpretation. -We therefore visualize the time-varying association of serum bilirubin -using the following piece of code:
+We therefore visualize the time-varying association of log serum +bilirubin with the hazard of the composite event using the following +piece of code:
x_times <- seq(0.001, 12, length = 501)
-X <- cbind(1, ns(x_times, knots = c(3, 6, 9), B = c(0, 14.5)))
+X <- cbind(1, ns(x_times, knots = c(3, 6, 9), B = c(0, 14.5)))
mcmc_alphas <- do.call('rbind', jointFit2$mcmc$alphas)
log_hr <- X %*% t(mcmc_alphas)
log_hr_mean <- rowMeans(log_hr)
@@ -302,7 +303,7 @@ Non Proportional Hazards
matplot(x_times, cbind(exp(log_hr_mean), exp(log_hr_low), exp(log_hr_upp)),
type = "l", col = c("red", "black", "black"), lty = c(1, 2, 2), lwd = 2,
- xlab = "Follow-up Time (years)", ylab = "Hazard Ratio serum Bilirubin",
+ xlab = "Follow-up Time (years)", ylab = "Hazard Ratio log serum Bilirubin",
ylim = c(0.5, 6.4))
abline(h = exp(coef(jointFit1)$association), lty = 2, col = "red")
abline(h = 1, lty = 2)
@@ -322,8 +323,8 @@ Non Proportional Hazards#>
#> The criteria are calculated based on the marginal log-likelihood.
Both the WAIC and LPML indicate that jointFit1 is a
-better model than jointFit2. The DIC is marginally smaller
-for jointFit2.
jointFit2. The DIC has the same magnitude
+for both models.
diff --git a/docs/articles/Time_Varying_Effects_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/Time_Varying_Effects_files/figure-html/unnamed-chunk-5-1.png
index 2d0372c..a792a25 100644
Binary files a/docs/articles/Time_Varying_Effects_files/figure-html/unnamed-chunk-5-1.png and b/docs/articles/Time_Varying_Effects_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/vignettes/Time_Varying_Effects.Rmd b/vignettes/Time_Varying_Effects.Rmd
index e5f799b..6fd34ef 100644
--- a/vignettes/Time_Varying_Effects.Rmd
+++ b/vignettes/Time_Varying_Effects.Rmd
@@ -40,7 +40,7 @@ jointFit1 <- jm(CoxFit, fm, time_var = "year")
summary(jointFit1)
```
-To specify that the association of serum bilirubin may change over time we include an interaction of this time-varying covariate with a natural cubic spline of time using function `ns()` from the **splines** package. **Important Note:** for this to work correctly we need to explicitly specify the knots and boundary knots for the B-splines basis, i.e., in the following example we set the internal knots at 3, 6 and 9 years, and the boundary knots at 0 and 14.5 years:
+To specify that the association of serum bilirubin may change over time we include an interaction of this time-varying covariate with a natural cubic spline of time using function `ns()` from the **splines** package. **Important Note:** for this to work correctly we need to explicitly specify the internal and boundary knots for the B-splines basis, i.e., in the following example we set the internal knots at 3, 6 and 9 years, and the boundary knots at 0 and 14.5 years:
```{r}
form_splines <- ~ value(log(serBilir)) * ns(year, k = c(3, 6, 9), B = c(0, 14.5))
jointFit2 <- update(jointFit1, functional_forms = form_splines,
@@ -48,7 +48,7 @@ jointFit2 <- update(jointFit1, functional_forms = form_splines,
summary(jointFit2)
```
-The spline coefficients do not have a straightforward interpretation. We therefore visualize the time-varying association of serum bilirubin using the following piece of code:
+The spline coefficients do not have a straightforward interpretation. We therefore visualize the time-varying association of log serum bilirubin with the hazard of the composite event using the following piece of code:
```{r, fig.align = "center", fig.width = 8.5, fig.height = 7.5}
x_times <- seq(0.001, 12, length = 501)
X <- cbind(1, ns(x_times, knots = c(3, 6, 9), B = c(0, 14.5)))
@@ -60,7 +60,7 @@ log_hr_upp <- apply(log_hr, 1, quantile, probs = 0.975)
matplot(x_times, cbind(exp(log_hr_mean), exp(log_hr_low), exp(log_hr_upp)),
type = "l", col = c("red", "black", "black"), lty = c(1, 2, 2), lwd = 2,
- xlab = "Follow-up Time (years)", ylab = "Hazard Ratio serum Bilirubin",
+ xlab = "Follow-up Time (years)", ylab = "Hazard Ratio log serum Bilirubin",
ylim = c(0.5, 6.4))
abline(h = exp(coef(jointFit1)$association), lty = 2, col = "red")
abline(h = 1, lty = 2)
@@ -73,4 +73,4 @@ We observe that the 95% credible interval for the time-varying coefficient inclu
compare_jm(jointFit1, jointFit2)
```
-Both the WAIC and LPML indicate that `jointFit1` is a better model than `jointFit2`. The DIC is marginally smaller for `jointFit2`.
+Both the WAIC and LPML indicate that `jointFit1` is a better model than `jointFit2`. The DIC has the same magnitude for both models.