updates

DESCRIPTION

@@ -8,7 +8,7 @@ Authors@R: c(
     person("Pedro", "Miranda Afonso", email = "[email protected]", role = "aut")
     )
 Maintainer: Dimitris Rizopoulos <[email protected]>
-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

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

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"

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
 

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 = "_")

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
 }
 

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.}
     }

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{