From a574e17209c9cfcdf0d4def1583ae02c203af60c Mon Sep 17 00:00:00 2001 From: Dimitris Rizopoulos Date: Wed, 29 May 2024 20:13:14 +0200 Subject: [PATCH] updates --- DESCRIPTION | 2 +- Development/CI/causal_effects.R | 4 ++- Development/predict/predict.R | 59 +++++++++++++++++++++++++-------- NEWS.md | 4 ++- R/predict_funs.R | 14 ++++---- man/JMbayes2.Rd | 2 +- man/jm.Rd | 2 +- man/predict.Rd | 4 +-- 8 files changed, 64 insertions(+), 27 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 91b9c21..e010e64 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,7 @@ Authors@R: c( person("Pedro", "Miranda Afonso", email = "p.mirandaafonso@erasmusmc.nl", role = "aut") ) Maintainer: Dimitris Rizopoulos -Date: 2024-05-23 +Date: 2024-05-30 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/Development/CI/causal_effects.R b/Development/CI/causal_effects.R index 64aa0f4..040cbf4 100644 --- a/Development/CI/causal_effects.R +++ b/Development/CI/causal_effects.R @@ -39,6 +39,7 @@ jointFit <- jm(CoxFit, lmeFit, time_var = "year", functional_forms = fForms, summary(jointFit) + #----------------------------- # Conditional Causal Effects - #----------------------------- @@ -51,6 +52,7 @@ newdataL <- pbc2[pbc2$id %in% c(81), ] newdataL$status2 <- 0 # The data.frame 'newdataE_withIE' contains the event information in the # presence of the IE. +pbc2_CR$serBilir <- 0.1 newdataE_withIE <- pbc2_CR[pbc2_CR$id == 81, ] newdataE_withIE$event <- 0 # The data.frame 'newdataE_withoutIE' contains the event information in the @@ -77,7 +79,7 @@ for (i in seq_along(t0)) { newdataE_withoutIE_i$stop <- t0[i] # In the event data.frame with the IE we set that the IE occurs a bit after t0 # and a bit afterward is the last time the patient was available - newdataE_withIE_i <- newdataE_withoutIE_i + newdataE_withIE_i <- newdataE_withIE newdataE_withIE_i$IE <- 1 # We calculate the predictions using the two datasets diff --git a/Development/predict/predict.R b/Development/predict/predict.R index 4a33be1..12253d4 100644 --- a/Development/predict/predict.R +++ b/Development/predict/predict.R @@ -32,8 +32,10 @@ if (FALSE) { file = "C:/Users/drizo/OneDrive/Desktop/predict_simpleExample.RData") load("C:/Users/drizo/OneDrive/Desktop/predict_simpleExample.RData") + load("C:/Users/drizo/OneDrive/Desktop/predict_CRExample.RData") data("pbc2", package = "JM") data("pbc2.id", package = "JM") + pbc2.idCR <- JMbayes2::crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", nameStrata = "CR") @@ -53,22 +55,53 @@ if (FALSE) { source("./R/help_functions.R") source("./R/predict_funs.R") Rcpp::sourceCpp('src/mcmc_fit.cpp') + + source("./Development/CI/prepare_data.R") + source("./Development/CI/causal_effects_fun.R") + newdataL <- pbc2[pbc2$id %in% c(81), ] + newdataL$status2 <- 0 + # The data.frame 'newdataE_withIE' contains the event information in the + # presence of the IE. + pbc2_CR$serBilir <- 0.1 + newdataE_withIE <- pbc2_CR[pbc2_CR$id == 81, ] + newdataE_withIE$event <- 0 + # The data.frame 'newdataE_withoutIE' contains the event information in the + # absence of the IE. + newdataE_withoutIE <- newdataE_withIE[c(1, 3), ] + + t0 <- 5 + + newdataL_i <- newdataL[newdataL$year <= t0, ] + # In the event data.frame without the IE we set the stop time at t0 + newdataE_withoutIE_i <- newdataE_withoutIE + newdataE_withoutIE_i$stop <- t0 + # In the event data.frame with the IE we set that the IE occurs a bit after t0 + # and a bit afterward is the last time the patient was available + newdataE_withIE_i <- newdataE_withIE + newdataE_withIE_i$stop[c(2, 4)] <- t0 + + # We calculate the predictions using the two datasets + newdata_withIE_i <- list(newdataL = newdataL_i, newdataE = newdataE_withIE_i) + newdata_withoutIE_i <- list(newdataL = newdataL_i, newdataE = newdataE_withoutIE_i) + } -object <- jointFit2 -ND <- pbc2[pbc2$id %in% c(2, 3), ] -ND$id <- factor(ND$id) -ND <- ND[ND$year < 1, ] -ND$status2 <- 0 -ND$years <- 1 #with(ND, ave(year, id, FUN = function (x) max(x, na.rm = T))) -ND. <- pbc2.idCR[pbc2.idCR$id %in% c(2, 3), ] -ND.$id <- factor(ND.$id) -ND.$status2 <- 0 -ND.$years <- 1 -newdata = list(newdataL = ND, newdataE = ND.) -newdata2 = NULL -times = c(2, 3) +object <- jointFit +control = NULL +#ND <- pbc2[pbc2$id %in% c(2, 3), ] +#ND$id <- factor(ND$id) +#ND <- ND[ND$year < 1, ] +#ND$status2 <- 0 +#ND$years <- 1 #with(ND, ave(year, id, FUN = function (x) max(x, na.rm = T))) +#ND. <- pbc2.idCR[pbc2.idCR$id %in% c(2, 3), ] +#ND.$id <- factor(ND.$id) +#ND.$status2 <- 0 +#ND.$years <- 1 +#newdata = list(newdataL = ND, newdataE = ND.) +newdata = newdata_withoutIE_i +newdata2 = newdata_withIE_i +times = c(6, 7, 8) process = "event" type_pred = "response" type = "subject_specific" diff --git a/NEWS.md b/NEWS.md index 009240c..7f1b18c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,9 +3,11 @@ ## Major * `jm()` now allows for zero-correlations constraints in the covariance matrix of the random effects. When the mixed models provided in the `Mixed_objects` argument have been fitted assuming a diagonal matrix for the random effects, this will also be assumed in the joint model (in previous versions, this was ignored). In addition, the new argument `which_independent` can be used to specify which longitudinal outcomes are to be assumed independent. +* `jm()` can fit joint models with a combination of interval-censored data and competing risks (e.g., one of the the competing events is interval-censored and the other(s) not). + * A bug in the `predict()` method causing low AUC values has been corrected. -* The time-varying ROC and AUC now allow to correct for censoring in the interval `Tstart` to `Thoriz` using inverse probability of censoring weighting. The default remains model-based weights +* The time-varying ROC and AUC now allow to correct for censoring in the interval `Tstart` to `Thoriz` using inverse probability of censoring weighting. The default remains model-based weights. # JMbayes2 0.4.1 diff --git a/R/predict_funs.R b/R/predict_funs.R index 5441244..4f4771d 100644 --- a/R/predict_funs.R +++ b/R/predict_funs.R @@ -968,14 +968,15 @@ predict_Event <- function (object, components_newdata, newdata, newdata2, # last_times to times numer_intgr <- function (low, upp, - type = c("7-Gauss-Kronrod", - "15-Gauss-Kronrod", "Simpson-3/8")) { + type = c("7-Gauss-Kronrod", "15-Gauss-Kronrod", + "Simpson-3/8")) { type <- match.arg(type) if (type == "Simpson-3/8") { n <- length(upp) - list(qpoints = c(rep(low, n), (2 * low + upp) / 3, + if (n > 1) low <- rep(low, length.out = n) + list(qpoints = cbind(low, (2 * low + upp) / 3, (low + 2 * upp) / 3, upp), - log_weights = log(c(1, 3, 3, 1) * (upp - low) / 8)) + log_weights = log(c(c(1, 3, 3, 1) %o% (upp - low)) / 8)) } else { GK <- if (type == "7-Gauss-Kronrod") gaussKronrod(7L) else gaussKronrod(15L) sk <- GK$sk @@ -984,7 +985,6 @@ predict_Event <- function (object, components_newdata, newdata, newdata2, log_Pwk <- unname(rep(log(P), each = length(sk)) + rep_len(log(GK$wk), length.out = length(st))) list(qpoints = st, log_weights = log_Pwk) - } } # extract strata @@ -1009,12 +1009,12 @@ predict_Event <- function (object, components_newdata, newdata, newdata2, lastRow_strt <- tapply(row.names(d), factor2(strt_i), tail, n = 1L) out <- vector("list", length(qpoints)) for (j in seq_along(qpoints)) { - d[[Time_var]] <- qpoints[j] + d[lastRow_strt, Time_var] <- qpoints[j] d[[".log_weights"]] <- log_w[j] d[[".time"]] <- tt[j] d[[".qpoint"]] <- j r <- vector("list", nstrt_i) - for (i in seq_along(strt_i)) { + for (i in seq_len(nstrt_i)) { dd <- d dd[lastRow_strt[i], event_var] <- 1 dd[[".id_h"]] <- paste(dd[[id_var]][1], j, i, sep = "_") diff --git a/man/JMbayes2.Rd b/man/JMbayes2.Rd index b63b503..b3274cc 100644 --- a/man/JMbayes2.Rd +++ b/man/JMbayes2.Rd @@ -16,7 +16,7 @@ Fit joint models for longitudinal and time-to-event data under the Bayesian appr Package: \tab JMbayes2\cr Type: \tab Package\cr Version: \tab 0.5-0\cr -Date: \tab 2024-04-24\cr +Date: \tab 2024-05-30\cr License: \tab GPL (>=3)\cr } diff --git a/man/jm.Rd b/man/jm.Rd index c892972..03d3568 100644 --- a/man/jm.Rd +++ b/man/jm.Rd @@ -162,7 +162,7 @@ tv(x, knots = NULL, ord = 2L) \item{\code{cores}}{an integer specifying the number of cores to use for running the chains in parallel; no point of setting this greater than \code{n_chains}.} \item{\code{parallel}}{a character string indicating how the parallel sampling of the chains will - be performed. Options are \code{"snow"} (default) and \code{multicore}.} + be performed. Options are \code{"snow"} (default) and \code{"multicore"}.} \item{\code{knots}}{a numeric vector with the position of the knots for the B-spline approximation of the log baseline hazard function.} } diff --git a/man/predict.Rd b/man/predict.Rd index f07e205..fc7a12c 100644 --- a/man/predict.Rd +++ b/man/predict.Rd @@ -76,7 +76,7 @@ subjects in \code{newdata}.} \item{n_mcmc}{the number of Metropolis-Hastings iterations for sampling the random effects per iteration of \code{n_samples}; only the last iteration is retained.} -\item{parallel}{character string; what type of parallel computing to use.} +\item{parallel}{character string; what type of parallel computing to use. Options are \code{"snow"} (default) and \code{"multicore"}.} \item{cores}{how many number of cores to use. If there more than 20 subjects in \code{newdata}, parallel computing is invoked with four cores by default. If \code{cores = 1}, no parallel computing is used.} @@ -105,7 +105,7 @@ longitudinal outcomes are plotted.} \item{ylim_long_outcome_range}{logical; if \code{TRUE}, the range of the y-axis spans across the range of the outcome in the data used to fit the model; not only the range of values of the specific subject being plotted.} -\item{\dots}{extra aguments passed to control.} +\item{\dots}{aguments passed to control.} } \details{