From 0d6b4459532869d399eb0874dd32bfce06c9ee90 Mon Sep 17 00:00:00 2001 From: Dimitris Rizopoulos Date: Fri, 2 Feb 2024 14:32:24 +0100 Subject: [PATCH] updates --- DESCRIPTION | 4 ++-- R/basic_methods.R | 9 +++++---- R/predict_funs.R | 5 +++-- man/predict.Rd | 4 +++- src/mcmc_fit.cpp | 15 +++++++++++---- 5 files changed, 24 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a8842cf..79ea7b5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ 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 = "d.rizopoulos@erasmusmc.nl", role = c("aut", "cre"), comment = c(ORCID = '0000-0001-9397-0900')), person("Grigorios", "Papageorgiou", email = "g.papageorgiou@erasmusmc.nl", @@ -9,7 +9,7 @@ Authors@R: c(person("Dimitris", "Rizopoulos", email = "d.rizopoulos@erasmusmc.nl person("Pedro", "Miranda Afonso", email = "p.mirandaafonso@erasmusmc.nl", role = "aut")) Maintainer: Dimitris Rizopoulos -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 index ce505b9..7edc911 100644 --- 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 index 38254ca..64e03e3 100644 --- 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 index 23d0a6d..64c2546 100644 --- 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 index ab3c2c5..0991358 100644 --- 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(control["use_Y"]); uword n_samples = as(control["n_samples"]); uword n_iter = as(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;