Causal Effects
+Dimitris +Rizopoulos
+ +2024-01-08
+ + Source:vignettes/Causal_Effects.Rmd
+ Causal_Effects.RmdCausal Effects from Joint Models +
+We will illustrate the calculation of causal effects from joint
+models using the PBC dataset for the longitudinal outcome
+serBilir and the composite event transplantation or death.
+We start by fitting a joint model to the data. In the longitudinal
+submodel, we specify nonlinear subject-specific trajectories using
+natural cubic splines. In the fixed-effects part, we also include the
+treatment effect and its interaction with time. In the survival
+submodel, we only include the treatment effect.
+pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
+lmeFit <- lme(log(serBilir) ~ ns(year, 3, B = c(0, 14.4)) * drug,
+ data = pbc2, random = ~ ns(year, 3, B = c(0, 14.4)) | id,
+ control = lmeControl(opt = "optim"))
+CoxFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id)
+jmFit <- jm(CoxFit, lmeFit, time_var = "year")
+summary(jmFit)
+#>
+#> Call:
+#> jm(Surv_object = CoxFit, Mixed_objects = lmeFit, time_var = "year")
+#>
+#> Data Descriptives:
+#> Number of Groups: 312 Number of events: 169 (54.2%)
+#> Number of Observations:
+#> log(serBilir): 1945
+#>
+#> DIC WAIC LPML
+#> marginal 4217.836 4527.502 -2937.053
+#> conditional 6503.006 6206.639 -3454.250
+#>
+#> Random-effects covariance matrix:
+#>
+#> StdDev Corr
+#> (Intr) 0.9958 (Intr) n(,3,B=c(0,14.4))1 n(,3,B=c(0,14.4))2
+#> n(,3,B=c(0,14.4))1 1.5005 0.2219
+#> n(,3,B=c(0,14.4))2 1.6970 0.4318 0.7505
+#> n(,3,B=c(0,14.4))3 1.8995 0.4247 0.2076 0.6798
+#>
+#> Survival Outcome:
+#> Mean StDev 2.5% 97.5% P Rhat
+#> drugD-penicil -0.0138 0.2061 -0.4079 0.3859 0.9533 1.0037
+#> value(log(serBilir)) 1.3001 0.0881 1.1354 1.4741 0.0000 1.0109
+#>
+#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
+#> Mean StDev 2.5% 97.5% P Rhat
+#> (Intercept) 0.5863 0.0808 0.4309 0.7430 0.0000 1.0007
+#> ns(,3,B=c(0,14.4))1 1.1251 0.1719 0.7905 1.4659 0.0000 1.0058
+#> ns(,3,B=c(0,14.4))2 2.1954 0.2329 1.7749 2.6913 0.0000 1.0201
+#> ns(,3,B=c(0,14.4))3 2.4177 0.3416 1.8350 3.1714 0.0000 1.0342
+#> drugD-penicil -0.1061 0.1163 -0.3386 0.1217 0.3633 1.0017
+#> n(,3,B=c(0,14.4))1: 0.1869 0.2369 -0.2802 0.6490 0.4218 1.0055
+#> n(,3,B=c(0,14.4))2: -0.4776 0.3224 -1.1454 0.1111 0.1211 1.0175
+#> n(,3,B=c(0,14.4))3: -0.6994 0.4770 -1.7031 0.1481 0.1233 1.0384
+#> sigma 0.2894 0.0065 0.2776 0.3027 0.0000 1.0061
+#>
+#> MCMC summary:
+#> chains: 3
+#> iterations per chain: 3500
+#> burn-in per chain: 500
+#> thinning: 1
+#> time: 20 secThe coefficient for drugD-penicil for the survival
+outcome in the output produced by the summary() method
+denotes the residual/direct effect of treatment on the risk of the
+composite event. It does not include the effect of treatment that
+follows via the serum bilirubin pathway.
We will illustrate the calculation of causal risk differences for the +group of patients that have the same distribution of serum bilirubin +values as Patient 2:
+
+xyplot(log(serBilir) ~ year, data = pbc2, subset = id == 2, type = "b",
+ xlab = "Follow-up time (years)", ylab = "log{serum bilirubin (mg/dL)}",
+ main = "Patient 2")
We calculate the risk difference for the composite event between the
+active treatment D-penicillamine and placebo at the horizon time
+t_horiz = 6 using the longitudinal data up to year
+t0 = 4. To achieve this, we create a dataset with this
+patient’s data. This patient received the active treatment
+D-penicillamine; hence, we also create a version of her data with the
+drug variable set to placebo:
+t0 <- 4
+t_horiz <- 6
+dataP2_Dpenici <- pbc2[pbc2$id == 2 & pbc2$year <= t0, ]
+dataP2_Dpenici$years <- t0
+dataP2_Dpenici$status2 <- 0
+
+dataP2_placebo <- dataP2_Dpenici
+dataP2_placebo$drug <- factor("placebo", levels = levels(pbc2$drug))Note that in the dataP2_placebo dataset, we need to
+specify that drug is a factor with two levels. We also
+specify that the last time point we know the patient was still
+event-free was t0.
We estimate the cumulative risk for the composite event at
+t_horiz under the active treatment arm using the
+predict() method:
+Pr1 <- predict(jmFit, newdata = dataP2_Dpenici, process = "event",
+ times = t_horiz, return_mcmc = TRUE)We have set the argument return_mcmc to
+TRUE to enable the calculation of a credible interval that
+accounts for the MCMC uncertainty. We produce the same estimate under
+the placebo arm:
+Pr0 <- predict(jmFit, newdata = dataP2_placebo, process = "event",
+ times = t_horiz, return_mcmc = TRUE)The estimated risk difference and its 95% credible interval are
+calculated by the corresponding elements of the Pr1 and
+Pr0 objects, i.e.,
+# estimate
+Pr1$pred[2L] - Pr0$pred[2L]
+#> [1] 0.002916962
+
+# MCMC variability
+quantile(Pr1$mcmc[2L, ] - Pr0$mcmc[2L, ], probs = c(0.025, 0.975))
+#> 2.5% 97.5%
+#> -0.1790865 0.1952537Time-varying treatments +
+An extended example with time-varying treatments / intermediate +events that showcases a calculation of the variance of the causal +effects that includes the sampling variability is available here.
+- + Causal Effects + +
- Competing Risks
-
@@ -132,7 +135,7 @@
Competing Risks
Dimitris Rizopoulos
-2023-10-12
+2024-01-08
Source:vignettes/Competing_Risks.Rmd@@ -293,7 +296,7 @@Competing_Risks.RmdFit models#> iterations per chain: 25000 #> burn-in per chain: 5000 #> thinning: 5 -#> time: 6.4 min +#> time: 6.6 min @@ -73,6 +73,9 @@
- + Causal Effects + +
- Competing Risks
-
@@ -132,7 +135,7 @@
Dynamic Predictions
Dimitris Rizopoulos
-2023-10-12
+2024-01-08
Source:vignettes/Dynamic_Predictions.Rmddiff --git a/docs/articles/JMbayes2.html b/docs/articles/JMbayes2.html index 33371c3..5a97c42 100644 --- a/docs/articles/JMbayes2.html +++ b/docs/articles/JMbayes2.html @@ -53,7 +53,7 @@ @@ -73,6 +73,9 @@Dynamic_Predictions.Rmd- + Causal Effects + +
- Competing Risks
-
@@ -132,7 +135,7 @@
Univariate and Multivariate Joint Models
Dimitris Rizopoulos
-2023-10-12
+2024-01-08
Source:vignettes/JMbayes2.Rmd@@ -212,7 +215,7 @@JMbayes2.RmdUnivariate#> iterations per chain: 3500 #> burn-in per chain: 500 #> thinning: 1 -#> time: 15 sec +#> time: 16 sec
The output of the
summary()method provides some descriptive statistics of the sample at hand, followed by some fit statistics based on the marginal (random effects are integrated out @@ -427,7 +430,7 @@Functional forms#> iterations per chain: 12000 #> burn-in per chain: 2000 #> thinning: 5 -#> time: 2 min +#> time: 1.9 min
As seen above, the
functional_formsargument is a named list with elements corresponding to the longitudinal outcomes. If a longitudinal outcome is not specified in this list, then the default, diff --git a/docs/articles/Multi_State_Processes.html b/docs/articles/Multi_State_Processes.html index 3b2ca8c..5ae9074 100644 --- a/docs/articles/Multi_State_Processes.html +++ b/docs/articles/Multi_State_Processes.html @@ -53,7 +53,7 @@ @@ -73,6 +73,9 @@- + Causal Effects + +
- Competing Risks
-
@@ -132,7 +135,7 @@
Non-Gaussian Mixed Models
Dimitris Rizopoulos
-2023-10-12
+2024-01-08
Source:vignettes/Non_Gaussian_Mixed_Models.Rmd@@ -731,7 +734,7 @@Non_Gaussian_Mixed_Models.RmdNegative binomial mixed models#> iterations per chain: 3500 #> burn-in per chain: 500 #> thinning: 1 -#> time: 26 sec +#> time: 25 sec @@ -873,7 +876,7 @@
Beta-binomial longitudinal outcomes #> iterations per chain: 3500 #> burn-in per chain: 500 #> thinning: 1 -#> time: 48 sec +#> time: 47 sec diff --git a/docs/articles/Recurring_Events.html b/docs/articles/Recurring_Events.html index fe4d9a5..af88f52 100644 --- a/docs/articles/Recurring_Events.html +++ b/docs/articles/Recurring_Events.html @@ -53,7 +53,7 @@ @@ -73,6 +73,9 @@
- + Causal Effects + +
- Competing Risks
-
@@ -132,7 +135,7 @@
Recurrent Events
Pedro Miranda Afonso
-2023-10-12
+2024-01-08
Source:vignettes/Recurring_Events.Rmd@@ -161,29 +164,22 @@Recurring_Events.RmdIntroduction -
Figure 1 Visual representation of an hazard function under the
++Figure 1 Visual representation of an hazard function under the gap or calendar timescale. During the folow-up the subject experienced three events. - +
The model also accommodates discontinuous risk intervals, i.e., periods in which the subject is not at risk of experiencing the recurring event (Figure 2). For example, while a patient is in the hospital, they are not at risk of being hospitalized again.
-
Figure 2 Visual representation of an hazard function under the
++Figure 2 Visual representation of an hazard function under the gap or calendar timescale, while accounting for non-risk periods (gray areas). During the folow-up the subject experienced three events. - +
A joint model with \(p\) normally distributed longitudinal markers, a terminal event process, and a @@ -245,8 +241,9 @@
Data failure-time. We assume that the underlying value of the longitudinal outcome is associated with both risk models and use the gap timescale. The reader can easily extend this example to accommodate -multiple longitudinal markers with other forms of association, exclude -the terminal event, or use a different timescale. +multiple longitudinal markers with other forms of association, include +competing risks, consider only the recurrent events process, or use a +different timescale.
gen_data <- function(){ n <- 500 # desired number of subjects @@ -305,14 +302,14 @@Data rmv_ids <- sample(which(ter_na), sum(ter_na) - n_target) ter_na[rmv_ids] <- FALSE # remove the excess of subjects ter <- data.frame(id = seq_len(n)[ter_na], - time = ter_times[ter_na], - group = group[ter_na], - age = age[ter_na]) + tstop = ter_times[ter_na], + group = group[ter_na], + age = age[ter_na]) frailty <- frailty[ter_na] b <- b[ter_na, , drop = FALSE] cens_times <- tmax - ter$status <- as.numeric(ter$time <= cens_times) # event indicator - ter$time <- pmin(ter$time, cens_times) # add censoring time + ter$status <- as.numeric(ter$tstop <= cens_times) # event indicator + ter$tstop <- pmin(ter$tstop, cens_times) # add censoring time remove(gammas_t, sigma_t, group, W_t, eta_t, alpha_t, invS_t, u_t, i, root, n_target, rmv_ids, ter_times, cens_times, n, alphaF, age, ter_na, sigmaF) @@ -352,7 +349,7 @@
Data j <- 1 for(i in seq_along(ter$id)) { tstart <- 0 - while(!is.na(tstart) & tstart < ter$time[i]) { + while(!is.na(tstart) & tstart < ter$tstop[i]) { u_r <- runif(1) root <- try(uniroot(invS_r, interval = c(1e-05, 250), # arbitrary upper limit u = u_r, i = i, tstart = tstart)$root, TRUE) @@ -371,7 +368,7 @@
Data rec$id <- match(rec$id, unique(rec$id)) # rename IDs rec$group <- ter$group[rec$id] rec$age <- ter$age[rec$id] - rec$Stime <- ter$time[rec$id] + rec$Stime <- ter$tstop[rec$id] rec$status <- as.numeric(!is.na(rec$tstop) & rec$tstop < rec$Stime) # event indicator rec$tstop <- pmin(rec$tstop, rec$Stime, na.rm = TRUE) # add cens time rec$Stime <- NULL @@ -388,7 +385,7 @@
Data data = long) eta_y <- as.vector(X %*% betas + rowSums(Z * b[long$id, ])) long$y <- rnorm(length(eta_y), eta_y, sigma_y) - long_cens <- long$time <= rep(ter$time, times = rle(long$id)$lengths) + long_cens <- long$time <= rep(ter$tstop, times = rle(long$id)$lengths) long <- long[long_cens, , drop = FALSE] # drop censored encounters remove(kn, Bkn, X, betas, Z, b, eta_y, sigma_y, n_i, tmax, long_cens, scale) ############################################################################## @@ -407,10 +404,9 @@
Data in the package the
rc_setup()function:cox_data <- rc_setup(rc_data = recu_data, trm_data = term_data, - rc_idVar = "id", rc_statusVar = "status", - rc_startVar = "tstart", rc_stopVar = "tstop", - trm_idVar = "id", trm_statusVar = "status", - trm_stopVar = "time", + idVar = "id", statusVar = "status", + startVar = "tstart", stopVar = "tstop", + trm_censLevel = 0, nameStrata = "strata", nameStatus = "status")Each subject has as many rows in the new data frame as the number of their recurrent risk periods plus one for the terminal event. The data @@ -425,26 +421,22 @@
Data
+#> 1 1 0.0000000 0.3756627 1 R +#> 2 1 0.4275324 0.7832841 1 R +#> 3 1 0.8724938 1.0863887 1 R +#> 4 1 1.1212953 1.8741434 1 R +#> 5 1 1.9372355 2.7843451 1 R +#> 6 1 2.7906559 3.4618622 1 R +#> 7 1 3.5166929 3.5830169 1 R +#> 8 1 3.6251219 4.0375415 0 R +#> 9 1 0.0000000 4.0375415 1 T1cox_data[cox_data$id == 1, c("id", "tstart", "tstop", "status", "strata")] #> id tstart tstop status strata -#> 1 1 0.0000000 0.3756627 1 Rec -#> 2 1 0.4275324 0.7832841 1 Rec -#> 3 1 0.8724938 1.0863887 1 Rec -#> 4 1 1.1212953 1.8741434 1 Rec -#> 5 1 1.9372355 2.7843451 1 Rec -#> 6 1 2.7906559 3.4618622 1 Rec -#> 7 1 3.5166929 3.5830169 1 Rec -#> 8 1 3.6251219 4.0375415 0 Rec -#> 9 1 0.0000000 4.0375415 1 Ter-@@ -73,6 +73,9 @@
Figure 3 Visual representation of the failure-time data during
++Figure 3 Visual representation of the failure-time data during the follow-up for subject 1. The horizontal black line denotes risk periods, while the blue line non-risk periods. ‘R’ and ‘T’ represent a recurrent and terminal event, respectively. - +
- + Causal Effects + +
- Competing Risks
-
@@ -133,7 +136,7 @@
Combined Dynamic Predictions via Super
Dimitris Rizopoulos
-2023-10-12
+2024-01-08
Source:vignettes/Super_Learning.Rmd@@ -297,8 +300,8 @@Super_Learning.RmdExample jmFit3 <- update(jmFit1, functional_forms = ~ value(log(serBilir)) + area(log(serBilir))) ### - lmeFit2 <- lme(log(serBilir) ~ ns(year, 3, B = c(0, 14.4)) + sex + age, - data = data, random = ~ ns(year, 3, B = c(0, 14.4)) | id, + lmeFit2 <- lme(log(serBilir) ~ ns(year, 2, B = c(0, 14.4)) + sex + age, + data = data, random = ~ ns(year, 2, B = c(0, 14.4)) | id, control = lmeControl(opt = "optim")) CoxFit2 <- coxph(Surv(years, status2) ~ sex + age, data = data_id) jmFit4 <- jm(CoxFit2, lmeFit2, time_var = "year") @@ -352,17 +355,17 @@
Example#> #> Cross-Validated Prediction Error using the Library of Joint Models 'Models_folds' #> -#> Super Learning Estimated Integrated Brier score: 0.056 +#> Super Learning Estimated Integrated Brier score: 0.0551 #> In the time interval: [6, 8) #> For the 166 subjects at risk at time 6 #> Number of subjects with an event in [6, 8): 18 #> Number of subjects with a censored time in [6, 8): 44 #> Accounting for censoring using model-based weights #> -#> Integrated Brier score per model: 0.0596 0.0588 0.0613 0.0629 0.065 -#> Weights per model: 0.1134 0.565 0.3029 0.0086 0.0101 +#> Integrated Brier score per model: 0.0601 0.0575 0.0626 0.0579 0.0611 +#> Weights per model: 0.2945 0.5574 0.0585 0.0669 0.0227 #> Number of folds: 5 -
We observe that the first three models dominate the weights. We also +
We observe that the first two models dominate the weights. We also note that the integrated Brier score based on the combined predictions is smaller than the integrated Brier score of each individual model. To calculate the model weights using the expected predictive cross-entropy @@ -375,16 +378,16 @@
Example#> #> Cross-Validated Expected Predictive Cross-Entropy using the Library of Joint Models 'Models_folds' #> -#> Super Learning Estimated EPCE: 0.3109 +#> Super Learning Estimated EPCE: 0.304 #> In the time interval: [6, 8) #> For the 166 subjects at risk at time 6 #> Number of subjects with an event in [6, 8): 18 #> Number of subjects with a censored time in [6, 8): 44 #> -#> EPCE per model: 0.3568 0.3599 0.3639 0.3992 0.4321 -#> Weights per model: 0.0028 0.5386 0.4586 0 0 +#> EPCE per model: 0.3545 0.348 0.3708 0.3617 0.413 +#> Weights per model: 0.0027 0.5659 0.4306 0 9e-04 #> Number of folds: 5 -
The EPCE results are inline with the ones from the integrated Brier +
The EPCE results are similar with the ones from the integrated Brier score; however, now only models
M2andM3share the weights. Again we see that the EPCE based on the combined cross-validated predictions is smaller than the EPCE based on the diff --git a/docs/articles/Super_Learning_files/figure-html/unnamed-chunk-1-1.png b/docs/articles/Super_Learning_files/figure-html/unnamed-chunk-1-1.png index af20efb..752c433 100644 Binary files a/docs/articles/Super_Learning_files/figure-html/unnamed-chunk-1-1.png and b/docs/articles/Super_Learning_files/figure-html/unnamed-chunk-1-1.png differ diff --git a/docs/articles/Time_Varying_Effects.html b/docs/articles/Time_Varying_Effects.html index b187e73..be5e76b 100644 --- a/docs/articles/Time_Varying_Effects.html +++ b/docs/articles/Time_Varying_Effects.html @@ -53,7 +53,7 @@ @@ -73,6 +73,9 @@- + Causal Effects + +
- Competing Risks
-
@@ -133,7 +136,7 @@
Time-Dependent Predictive Accuracy Measures for Joint Models
parallel = c("snow", "multicore"), cores = parallelly::availableCores(omit = 1L), ...) -create_folds(data, V = 5, id_var = "id", seed = 123L) +create_folds(data, V = 5, id_var = "id", strata = NULL, seed = 123L)@@ -217,6 +220,10 @@@@ -42,6 +42,9 @@Arguments
character string denoting the name of the subject id variable in
data.- strata
+character vector with the names of stratifying variables.
- seed
integer denoting the seed.
- + Causal Effects + +
- Competing Risks
-
diff --git a/docs/reference/coda_methods.html b/docs/reference/coda_methods.html
index 91dd0d1..182c793 100644
--- a/docs/reference/coda_methods.html
+++ b/docs/reference/coda_methods.html
@@ -23,7 +23,7 @@
@@ -41,6 +41,9 @@
- + Causal Effects + +
- Competing Risks
-
diff --git a/docs/reference/cr_setup.html b/docs/reference/cr_setup.html
index 082da6a..7fdbc75 100644
--- a/docs/reference/cr_setup.html
+++ b/docs/reference/cr_setup.html
@@ -25,7 +25,7 @@
@@ -43,6 +43,9 @@
- + Causal Effects + +
- Competing Risks
-
diff --git a/docs/reference/index.html b/docs/reference/index.html
index 66351fd..8cae925 100644
--- a/docs/reference/index.html
+++ b/docs/reference/index.html
@@ -23,7 +23,7 @@
@@ -41,6 +41,9 @@
- + Causal Effects + +
- Competing Risks
-
diff --git a/docs/reference/jm.html b/docs/reference/jm.html
index cbbff24..327b763 100644
--- a/docs/reference/jm.html
+++ b/docs/reference/jm.html
@@ -23,7 +23,7 @@
@@ -41,6 +41,9 @@
- + Causal Effects + +
- Competing Risks
-
@@ -493,7 +496,7 @@
Examples
#> iterations per chain: 11000 #> burn-in per chain: 1000 #> thinning: 1 -#> time: 35 sec +#> time: 34 sec traceplot(joint_model_fit_1)
@@ -611,7 +614,7 @@ Examples
#> iterations per chain: 11000 #> burn-in per chain: 1000 #> thinning: 1 -#> time: 1.4 min +#> time: 1.3 min traceplot(joint_model_fit_2)
diff --git a/docs/reference/methods.html b/docs/reference/methods.html
index c731cc4..9fa4383 100644
--- a/docs/reference/methods.html
+++ b/docs/reference/methods.html
@@ -23,7 +23,7 @@
@@ -41,6 +41,9 @@