updates

drizopoulos · Feb 2, 2024 · 0d6b445 · 0d6b445
1 parent d4a6226
commit 0d6b445
Show file tree

Hide file tree

Showing 5 changed files with 24 additions and 13 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,15 +1,15 @@
 Package: JMbayes2
 Type: Package
 Title: Extended Joint Models for Longitudinal and Time-to-Event Data
-Version: 0.4-8
+Version: 0.4-9
 Authors@R: c(person("Dimitris", "Rizopoulos", email = "[email protected]",
                   role = c("aut", "cre"), comment = c(ORCID = '0000-0001-9397-0900')),
            person("Grigorios", "Papageorgiou", email = "[email protected]",
                   role = "aut"),
            person("Pedro", "Miranda Afonso", email = "[email protected]",
                   role = "aut"))
 Maintainer: Dimitris Rizopoulos <[email protected]>
-Date: 2024-01-02
+Date: 2024-02-02
 BugReports: https://github.com/drizopoulos/JMbayes2/issues
 Description: Fit joint models for longitudinal and time-to-event data under the Bayesian approach. Multiple longitudinal outcomes of mixed type (continuous/categorical) and multiple event times (competing risks and multi-state processes) are accommodated. Rizopoulos (2012, ISBN:9781439872864).
 Suggests: lattice, knitr, rmarkdown, pkgdown

diff --git a/R/basic_methods.R b/R/basic_methods.R
@@ -600,7 +600,7 @@ predict.jm <- function (object, newdata = NULL, newdata2 = NULL,
                         process = c("longitudinal", "event"),
                         type_pred = c("response", "link"),
                         type = c("subject_specific", "mean_subject"),
-                        level = 0.95, return_newdata = FALSE,
+                        level = 0.95, return_newdata = FALSE, use_Y = TRUE,
                         return_mcmc = FALSE, n_samples = 200L, n_mcmc = 55L,
                         parallel = c("snow", "multicore"),
                         cores = NULL, seed = 123L, ...) {
@@ -735,8 +735,9 @@ predict.jm <- function (object, newdata = NULL, newdata2 = NULL,
         else length(unique(newdata[[id_var]]))
         cores <- if (n > 20) 4L else 1L
     }
-    components_newdata <- get_components_newdata(object, newdata, n_samples,
-                                                 n_mcmc, parallel, cores, seed)
+    components_newdata <-
+        get_components_newdata(object, newdata, n_samples,
+                               n_mcmc, parallel, cores, seed, use_Y)
     if (process == "longitudinal") {
         predict_Long(object, components_newdata, newdata, newdata2, times,
                      times_per_id, type, type_pred, level, return_newdata,
@@ -1034,7 +1035,7 @@ rc_setup <- function(rc_data, trm_data,
   trm_data <- trm_data[order(trm_data[[idVar]]), ]
   if(any(rc_data[[stopVar]] > trm_data[[stopVar]][rc_data[[idVar]]])) {
     stop(paste0("'", stopVar, "' in the recurring event data cannot be larger than '", stopVar," in the terminal event data.'"))
-  } 
+  }
   # create new dataset
   ## CR dataset
   n <- nrow(trm_data)

diff --git a/R/predict_funs.R b/R/predict_funs.R
@@ -533,7 +533,7 @@ prepare_DataE_preds <- function (object, newdataL, newdataE,
 }
 
 get_components_newdata <- function (object, newdata, n_samples, n_mcmc,
-                                    parallel, cores, seed) {
+                                    parallel, cores, seed, use_Y = TRUE) {
     # prepare the data for calculations
     newdataL <- if (!is.data.frame(newdata)) newdata[["newdataL"]] else newdata
     newdataE <- if (!is.data.frame(newdata)) newdata[["newdataE"]] else newdata
@@ -577,7 +577,8 @@ get_components_newdata <- function (object, newdata, n_samples, n_mcmc,
                 "MCMC iterations in the fitted model.")
         n_samples <- M
     }
-    control <- list(GK_k = object$control$GK_k, n_samples = n_samples, n_iter = n_mcmc)
+    control <- list(GK_k = object$control$GK_k, n_samples = n_samples,
+                    n_iter = n_mcmc, use_Y = use_Y)
     id_samples <-
         split(seq_len(control$n_samples),
               rep(seq_len(cores), each = ceiling(control$n_samples / cores),

diff --git a/man/predict.Rd b/man/predict.Rd
@@ -18,7 +18,7 @@ Predict method for object of class \code{"jm"}.
   times_per_id = FALSE, process = c("longitudinal", "event"),
   type_pred = c("response", "link"),
   type = c("subject_specific", "mean_subject"),
-  level = 0.95, return_newdata = FALSE, return_mcmc = FALSE,
+  level = 0.95, return_newdata = FALSE, use_Y = TRUE, return_mcmc = FALSE,
   n_samples = 200L, n_mcmc = 55L, parallel = c("snow", "multicore"),
   cores = NULL, seed = 123L,
   \dots)
@@ -67,6 +67,8 @@ subjects in \code{newdata}.}
 
 \item{return_newdata}{logical; should \code{predict()} return the predictions as extra columns in \code{newdata} and \code{newdata2}.}
 
+\item{use_Y}{logical; should the longitudinal measurements be used in the posterior of the random effects.}
+
 \item{return_mcmc}{logical; if \code{TRUE} the mcmc sample for the predictions is returned. It can be \code{TRUE} only in conjuction with \code{return_newdata} being \code{FALSE}.}
 
 \item{n_samples}{the number of samples to use from the original MCMC sample of \code{object}.}

diff --git a/src/mcmc_fit.cpp b/src/mcmc_fit.cpp
@@ -285,7 +285,7 @@ List mcmc_cpp (List model_data, List model_info, List initial_values,
   if(any_terminal) {
     for (uword j = 0; j < n_strata - 1; ++j) {
       alphaF_H.rows(which_term_H.at(j)).fill(alphaF.at(j));
-      alphaF_h.rows(which_term_h.at(j)).fill(alphaF.at(j)); 
+      alphaF_h.rows(which_term_h.at(j)).fill(alphaF.at(j));
     }
   }
   vec frailty_H(WH_gammas.n_rows, fill::zeros);
@@ -528,9 +528,9 @@ List mcmc_cpp (List model_data, List model_info, List initial_values,
                  W0H_bs_gammas, W0h_bs_gammas, W0H2_bs_gammas,
                  WH_gammas, Wh_gammas, WH2_gammas,
                  log_Pwk, log_Pwk2, id_H_fast, id_h_fast, which_event,
-                 which_right_event, which_left, which_interval, unq_idL, 
-                 n_burnin, recurrent, frailtyH_sigmaF_alphaF, 
-                 frailtyh_sigmaF_alphaF, save_random_effects, res_b, res_b_last, 
+                 which_right_event, which_left, which_interval, unq_idL,
+                 n_burnin, recurrent, frailtyH_sigmaF_alphaF,
+                 frailtyh_sigmaF_alphaF, save_random_effects, res_b, res_b_last,
                  cumsum_b, outprod_b, n_iter);
 
         // update intercepts
@@ -908,6 +908,7 @@ arma::cube simulate_REs (List Data, List MCMC, List control) {
   //////////////////////////////
   // Sampling Random Effects //
   /////////////////////////////
+  bool use_Y = as<bool>(control["use_Y"]);
   uword n_samples = as<uword>(control["n_samples"]);
   uword n_iter = as<uword>(control["n_iter"]);
   uword n_b = b_mat.n_rows;
@@ -980,6 +981,9 @@ arma::cube simulate_REs (List Data, List MCMC, List control) {
     ///
     vec logLik_re = log_re(b_mat, L_it, sds_it);
     // calculate the denominator
+    if (!use_Y) {
+        logLik_long = 0.0 * logLik_long;
+    }
     vec denominator_b =
       logLik_long + logLik_surv + logLik_re;
     for (uword i = 0; i < n_iter; ++i) {
@@ -1029,6 +1033,9 @@ arma::cube simulate_REs (List Data, List MCMC, List control) {
         //
         vec logLik_re_proposed = log_re(proposed_b_mat, L_it, sds_it);
         //
+        if (!use_Y) {
+            logLik_long_proposed = 0.0 * logLik_long_proposed;
+        }
         vec numerator_b =
           logLik_long_proposed + logLik_surv_proposed + logLik_re_proposed;
         vec log_ratio = numerator_b - denominator_b;