Support censoring (in latent projection)

.Rbuildignore

@@ -14,3 +14,4 @@
 ^CRAN-SUBMISSION$
 ^docs$
 ^_pkgdown\.yml$
+^pkgdown$

DESCRIPTION

@@ -59,6 +59,7 @@ Imports:
     mclogit
 Suggests:
     ggrepel,
+    ggfortify,
     rstanarm,
     brms,
     nlme,
@@ -88,5 +89,5 @@ Additional_repositories:
     https://stan-dev.r-universe.dev/
 LazyData: TRUE
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 VignetteBuilder: knitr, rmarkdown

NEWS.md

@@ -4,11 +4,13 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 # projpred 2.9.1.9000 (development)
 
-
 ## Major changes
 
+* Added support for censored observations when using the latent projection (with response-scale analyses). This makes it possible, e.g., to use the latent projection for time-to-event models (also known as models for survival analysis). Due to this new feature, the function passed to argument `latent_ll_oscale` of `extend_family()` now needs to have an argument `cens`. See `?extend_family` (section "Latent projection") as well as the latent-projection vignette (section ["Censored observations (survival analysis)"](https://mc-stan.org/projpred/articles/latent.html#cens)) for more information and examples. Note that only the performance statistics `"elpd"`, `"mlpd"`, and `"gmpd"` take censoring into account (on response scale). (GitHub: #528)
+
 ## Minor changes
 
+
 ## Bug fixes
 
 

R/cv_varsel.R

@@ -633,8 +633,11 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
     }
     loglik_forPSIS <- refmodel$family$latent_ll_oscale(
       mu_offs_oscale, dis = refmodel$dis, y_oscale = refmodel$y_oscale,
-      wobs = refmodel$wobs, cl_ref = seq_along(refmodel$wdraws_ref),
-      wdraws_ref = refmodel$wdraws_ref
+      wobs = refmodel$wobs,
+      cens = eval_el2_not_null(attr(refmodel$family$latent_ll_oscale,
+                                    "cens_var"),
+                               refmodel$fetch_data()),
+      cl_ref = seq_along(refmodel$wdraws_ref), wdraws_ref = refmodel$wdraws_ref
     )
     if (!is.matrix(loglik_forPSIS)) {
       stop("Unexpected structure for the output of `latent_ll_oscale`.")
@@ -810,6 +813,9 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
       log_lik_ref <- refmodel$family$latent_ll_oscale(
         refdist_eval_mu_offs_oscale, dis = refdist_eval$dis,
         y_oscale = refmodel$y_oscale[inds], wobs = refmodel$wobs[inds],
+        cens = eval_el2_not_null(attr(refmodel$family$latent_ll_oscale,
+                                      "cens_var"),
+                                 refmodel$fetch_data(obs = inds)),
         cl_ref = refdist_eval$cl, wdraws_ref = refdist_eval$wdraws_orig
       )
       if (all(is.na(log_lik_ref))) {
@@ -918,6 +924,9 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
         log_lik_sub_oscale <- refmodel$family$latent_ll_oscale(
           mu_k_oscale, dis = perf_eval_out[["dis_sub"]][[k]],
           y_oscale = refmodel$y_oscale[inds], wobs = refmodel$wobs[inds],
+          cens = eval_el2_not_null(attr(refmodel$family$latent_ll_oscale,
+                                        "cens_var"),
+                                   refmodel$fetch_data(obs = inds)),
           cl_ref = refdist_eval$cl, wdraws_ref = refdist_eval$wdraws_orig
         )
         loo_sub_oscale[[k]][inds] <- apply(log_lik_sub_oscale + lw_sub, 2,
@@ -1264,6 +1273,9 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
       loglik_mlvlRan <- refmodel$family$latent_ll_oscale(
         mu_offs_mlvlRan_oscale_odim, dis = refmodel$dis,
         y_oscale = refmodel$y_oscale, wobs = refmodel$wobs,
+        cens = eval_el2_not_null(attr(refmodel$family$latent_ll_oscale,
+                                      "cens_var"),
+                                 refmodel$fetch_data()),
         cl_ref = seq_along(refmodel$wdraws_ref),
         wdraws_ref = refmodel$wdraws_ref
       )
@@ -1434,6 +1446,7 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws, nclusters,
     mu_test <- fold$refmodel$family$linkinv(eta_test)
     summaries_ref <- weighted_summary_means(
       y_wobs_test = y_wobs_test[fold$omitted, , drop = FALSE],
+      data_aux_test = refmodel$fetch_data(obs = fold$omitted),
       family = fold$refmodel$family,
       wdraws = fold$refmodel$wdraws_ref,
       mu = mu_test,

R/extend_family.R

@@ -156,8 +156,8 @@
 #' The function supplied to argument `latent_ll_oscale` needs to have the
 #' prototype
 #' ```{r, eval = FALSE}
-#' latent_ll_oscale(ilpreds, dis, y_oscale, wobs = rep(1, length(y_oscale)),
-#'                  cl_ref, wdraws_ref = rep(1, length(cl_ref)))
+#' latent_ll_oscale(ilpreds, dis, y_oscale, wobs = rep(1, ncol(ilpreds)),
+#'                  cens, cl_ref, wdraws_ref = rep(1, length(cl_ref)))
 #' ```
 #' where:
 #' * `ilpreds` accepts the return value from `latent_ilink`.
@@ -167,18 +167,38 @@
 #' the original response scale.
 #' * `wobs` accepts a numeric vector of length \eqn{N} containing observation
 #' weights.
+#' * `cens` accepts a vector containing censoring indicators for the
+#' observations for which to calculate the response-scale log-likelihood values
+#' (i.e., for the observations from the second dimension of `ilpreds`). This is
+#' only relevant if attribute `cens_var` of `latent_ll_oscale` is not `NULL`
+#' (see below).
 #' * `cl_ref` accepts the same input as argument `cl_ref` of `latent_ilink`.
 #' * `wdraws_ref` accepts the same input as argument `wdraws_ref` of
 #' `latent_ilink`.
 #'
+#' In case of censoring (in the response values, i.e., survival or time-to-event
+#' analysis), the latent projection (with response-scale analyses) can be used
+#' by setting an attribute `cens_var` of the `latent_ll_oscale` function to a
+#' right-hand side formula with the name of the variable containing the
+#' censoring indicators (e.g., `0` = uncensored, `1` = censored) on its
+#' right-hand side. This variable named in the `cens_var` attribute is then
+#' retrieved (internally, whenever calling the `latent_ll_oscale` function) from
+#' the original dataset (possibly subsetted to the observations corresponding to
+#' the second dimension of `ilpreds`), `newdata`, or element `data` from
+#' [varsel()]'s argument `d_test`, whichever is applicable. The content of the
+#' retrieved variable is passed to argument `cens` of the `latent_ll_oscale`
+#' function. Note that only the performance statistics `"elpd"`, `"mlpd"`, and
+#' `"gmpd"` take censoring into account (on response scale).
+#'
 #' The return value of `latent_ll_oscale` needs to be an \eqn{S \times N}{S x N}
 #' matrix containing the response-scale (not latent-scale) log-likelihood values
 #' for the \eqn{N} observations from its inputs.
 #'
 #' The function supplied to argument `latent_ppd_oscale` needs to have the
 #' prototype
 #' ```{r, eval = FALSE}
-#' latent_ppd_oscale(ilpreds_resamp, dis_resamp, wobs, cl_ref,
+#' latent_ppd_oscale(ilpreds_resamp, dis_resamp,
+#'                   wobs = rep(1, ncol(ilpreds_resamp)), cl_ref,
 #'                   wdraws_ref = rep(1, length(cl_ref)), idxs_prjdraws)
 #' ```
 #' where:
@@ -322,6 +342,16 @@ extend_family <- function(family,
                   "work on response scale (only on latent scale).")
         }
       }
+      # Check for argument `dis` because it was added in v2.9.0 and hence after
+      # the introduction of the latent projection from v2.4.0:
+      if (!any(c("dis", "...") %in% methods::formalArgs(latent_ll_oscale))) {
+        stop("Please add argument `dis` to the `latent_ll_oscale` function.")
+      }
+      # Check for argument `cens` because it was added after v2.9.1 and hence
+      # after the introduction of the latent projection from v2.4.0:
+      if (!any(c("cens", "...") %in% methods::formalArgs(latent_ll_oscale))) {
+        stop("Please add argument `cens` to the `latent_ll_oscale` function.")
+      }
       if (is.null(latent_ppd_oscale)) {
         if (!is.null(family$cats)) {
           latent_ppd_oscale <- latent_ppd_oscale_cats
@@ -337,6 +367,13 @@ extend_family <- function(family,
                   "response scale (only on latent scale).")
         }
       }
+      # Check for argument `dis_resamp` because it was added in v2.9.0 and hence
+      # after the introduction of the latent projection from v2.4.0:
+      if (!any(c("dis_resamp", "...") %in%
+               methods::formalArgs(latent_ppd_oscale))) {
+        stop("Please add argument `dis_resamp` to the `latent_ppd_oscale` ",
+             "function.")
+      }
       family$latent_ilink <- latent_ilink
       family$latent_ll_oscale <- latent_ll_oscale
       family$latent_ppd_oscale <- latent_ppd_oscale

R/formula.R

@@ -760,6 +760,16 @@ eval_el2 <- function(formula, data) {
   eval(formula[[2]], data, environment(formula))
 }
 
+# Same as eval_el2(), but avoids evaluation of argument `data` if `formula` is
+# `NULL`:
+eval_el2_not_null <- function(formula, data) {
+  if (!is.null(formula)) {
+    return(eval_el2(formula = formula, data = data))
+  } else {
+    return(NULL)
+  }
+}
+
 ## Extract left hand side of a formula as a formula itself by removing the right
 ## hand side.
 ## @param x Formula

R/latent.R

@@ -13,14 +13,15 @@
 latent_ll_oscale_cats <- function(ilpreds,
                                   dis = rep(NA, nrow(ilpreds)),
                                   y_oscale,
-                                  wobs = rep(1, length(y_oscale)),
+                                  wobs = rep(1, ncol(ilpreds)),
+                                  cens,
                                   cl_ref,
                                   wdraws_ref = rep(1, length(cl_ref))) {
   return(ll_cats(ilpreds, margin_draws = 1, y = y_oscale, wobs = wobs))
 }
 latent_ppd_oscale_cats <- function(ilpreds_resamp,
                                    dis_resamp = rep(NA, nrow(ilpreds_resamp)),
-                                   wobs,
+                                   wobs = rep(1, ncol(ilpreds_resamp)),
                                    cl_ref,
                                    wdraws_ref = rep(1, length(cl_ref)),
                                    idxs_prjdraws) {
@@ -32,7 +33,8 @@ latent_ppd_oscale_cats <- function(ilpreds_resamp,
 latent_ll_oscale_binom_nocats <- function(ilpreds,
                                           dis = rep(NA, nrow(ilpreds)),
                                           y_oscale,
-                                          wobs = rep(1, length(y_oscale)),
+                                          wobs = rep(1, ncol(ilpreds)),
+                                          cens,
                                           cl_ref,
                                           wdraws_ref = rep(1, length(cl_ref))) {
   # Ensure finite log() values:
@@ -45,7 +47,7 @@ latent_ll_oscale_binom_nocats <- function(ilpreds,
 }
 latent_ppd_oscale_binom_nocats <- function(ilpreds_resamp,
                                            dis_resamp = rep(NA, nrow(ilpreds_resamp)),
-                                           wobs,
+                                           wobs = rep(1, ncol(ilpreds_resamp)),
                                            cl_ref,
                                            wdraws_ref = rep(1, length(cl_ref)),
                                            idxs_prjdraws) {
@@ -59,15 +61,16 @@ latent_ppd_oscale_binom_nocats <- function(ilpreds_resamp,
 latent_ll_oscale_poiss <- function(ilpreds,
                                    dis = rep(NA, nrow(ilpreds)),
                                    y_oscale,
-                                   wobs = rep(1, length(y_oscale)),
+                                   wobs = rep(1, ncol(ilpreds)),
+                                   cens,
                                    cl_ref,
                                    wdraws_ref = rep(1, length(cl_ref))) {
   ll_unw <- dpois(y_oscale, lambda = t(ilpreds), log = TRUE)
   return(t(wobs * ll_unw))
 }
 latent_ppd_oscale_poiss <- function(ilpreds_resamp,
                                     dis_resamp = rep(NA, nrow(ilpreds_resamp)),
-                                    wobs,
+                                    wobs = rep(1, ncol(ilpreds_resamp)),
                                     cl_ref,
                                     wdraws_ref = rep(1, length(cl_ref)),
                                     idxs_prjdraws) {
@@ -82,14 +85,15 @@ latent_ppd_oscale_poiss <- function(ilpreds_resamp,
 latent_ll_oscale_NA <- function(ilpreds,
                                 dis = rep(NA, nrow(ilpreds)),
                                 y_oscale,
-                                wobs = rep(1, length(y_oscale)),
+                                wobs = rep(1, ncol(ilpreds)),
+                                cens,
                                 cl_ref,
                                 wdraws_ref = rep(1, length(cl_ref))) {
   return(array(dim = dim(ilpreds)[1:2]))
 }
 latent_ppd_oscale_NA <- function(ilpreds_resamp,
                                  dis_resamp = rep(NA, nrow(ilpreds_resamp)),
-                                 wobs,
+                                 wobs = rep(1, ncol(ilpreds_resamp)),
                                  cl_ref,
                                  wdraws_ref = rep(1, length(cl_ref)),
                                  idxs_prjdraws) {

R/methods.R

@@ -352,8 +352,9 @@ proj_linpred_aux <- function(proj, newdata, offsetnew, weightsnew,
       ynew <- eval_lhs(formula = proj$refmodel$formula, data = newdata)
     }
   }
-  lpd_out <- compute_lpd(ynew = ynew, pred_sub = pred_sub, proj = proj,
-                         weights = weights, transformed = transform)
+  lpd_out <- compute_lpd(ynew = ynew, newdata = newdata, pred_sub = pred_sub,
+                         proj = proj, weights = weights,
+                         transformed = transform)
   if (integrated) {
     if (proj$refmodel$family$for_latent && transform &&
         length(dim(pred_sub)) == 3) {
@@ -409,7 +410,7 @@ proj_linpred_aux <- function(proj, newdata, offsetnew, weightsnew,
   return(nlist(pred = pred_sub, lpd = lpd_out))
 }
 
-compute_lpd <- function(ynew, pred_sub, proj, weights, transformed) {
+compute_lpd <- function(ynew, newdata, pred_sub, proj, weights, transformed) {
   if (!is.null(ynew)) {
     ## compute also the log-density
     target <- get_standard_y(ynew, weights, proj$refmodel$family)
@@ -446,6 +447,9 @@ compute_lpd <- function(ynew, pred_sub, proj, weights, transformed) {
     if (proj$refmodel$family$for_latent && transformed) {
       ll_oscale_out <- proj$refmodel$family$latent_ll_oscale(
         pred_sub, dis = proj$dis, y_oscale = ynew, wobs = weights,
+        cens = eval_el2_not_null(attr(proj$refmodel$family$latent_ll_oscale,
+                                      "cens_var"),
+                                 newdata %||% proj$refmodel$fetch_data()),
         cl_ref = proj$cl_ref, wdraws_ref = proj$wdraws_ref
       )
       if (!is.matrix(ll_oscale_out)) {

R/misc.R

@@ -193,6 +193,13 @@ validate_vsel_object_stats <- function(object, stats, resp_oscale = TRUE) {
              "latent Gaussian distribution is used there).")
       }
     }
+    if (object$refmodel$family$for_latent &&
+        !is.null(attr(object$refmodel$family$latent_ll_oscale, "cens_var")) &&
+        resp_oscale &&
+        !stat %in% c("elpd", "mlpd", "gmpd")) {
+      warning("Performance statistic `\"", stat, "\"` does not take ",
+              "censoring into account.")
+    }
   }
   return(invisible(TRUE))
 }

R/projfun.R

@@ -129,6 +129,8 @@ perf_eval <- function(search_path,
       sub_summary <- weighted_summary_means(
         y_wobs_test = data.frame(y = y_test, y_oscale = y_oscale_test,
                                  wobs = wobs_test),
+        data_aux_test = newdata_test %||%
+          refmodel_fulldata$fetch_data(obs = indices_test),
         family = refmodel_fulldata$family,
         wdraws = submodl$wdraws_prj,
         mu = refmodel_fulldata$mu_fun(submodl$outdmin,

R/refmodel.R

@@ -586,6 +586,9 @@ predict.refmodel <- function(object, newdata = NULL, ynew = NULL,
       }
       loglik <- refmodel$family$latent_ll_oscale(
         mu_oscale, dis = refmodel$dis, y_oscale = ynew, wobs = weightsnew,
+        cens = eval_el2_not_null(attr(refmodel$family$latent_ll_oscale,
+                                      "cens_var"),
+                                 newdata %||% refmodel$fetch_data()),
         cl_ref = seq_along(refmodel$wdraws_ref),
         wdraws_ref = rep(1, length(refmodel$wdraws_ref))
       )

R/summary_funs.R

@@ -7,6 +7,9 @@
 #   weights). In case of the latent projection, this `list` (or `data.frame`)
 #   also needs to contain `y_oscale` (response values on the original response
 #   scale, i.e., the non-latent response values).
+# @param data_aux_test Used for evaluating `attr(family$latent_ll_oscale,
+#   "cens_var")` when constructing the input for argument `cens` of the
+#   `family$latent_ll_oscale` function.
 # @param family A `family` object.
 # @param wdraws A vector of weights for the parameter draws.
 # @param mu A matrix of expected values for `y`.
@@ -28,7 +31,8 @@
 #
 # @return A `list` with elements `mu` and `lppd` which are both vectors
 #   containing the values for the quantities from the description above.
-weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
+weighted_summary_means <- function(y_wobs_test, data_aux_test, family, wdraws,
+                                   mu, dis, cl_ref,
                                    wdraws_ref = rep(1, length(cl_ref))) {
   if (!is.matrix(mu) || any(dim(mu) == 0)) {
     stop("Unexpected structure for `mu`. Do the return values of ",
@@ -54,7 +58,10 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
     }
     loglik_oscale <- family$latent_ll_oscale(
       mu_oscale, dis = dis, y_oscale = y_wobs_test$y_oscale,
-      wobs = y_wobs_test$wobs, cl_ref = cl_ref, wdraws_ref = wdraws_ref
+      wobs = y_wobs_test$wobs,
+      cens = eval_el2_not_null(attr(family$latent_ll_oscale, "cens_var"),
+                               data_aux_test),
+      cl_ref = cl_ref, wdraws_ref = wdraws_ref
     )
     if (!is.matrix(loglik_oscale)) {
       stop("Unexpected structure for the output of `latent_ll_oscale`.")

R/varsel.R

@@ -118,6 +118,10 @@
 #' If not `NULL`, then `d_test` needs to be a `list` with the following
 #' elements:
 #' * `data`: a `data.frame` containing the predictor variables for the test set.
+#' In case of the latent projection, this `data.frame` is also used for
+#' evaluating attribute `cens_var` of the `latent_ll_oscale` function, so if
+#' `cens_var` is not `NULL`, `data` also needs to contain the data for the
+#' variable from `cens_var`.
 #' * `offset`: a numeric vector containing the offset values for the test set
 #' (if there is no offset, use a vector of zeros).
 #' * `weights`: a numeric vector containing the observation weights for the test
@@ -460,9 +464,10 @@ varsel.refmodel <- function(
       mu_test <- refmodel$family$linkinv(eta_test)
     }
     ref <- weighted_summary_means(
-      y_wobs_test = y_wobs_test, family = refmodel$family,
-      wdraws = refmodel$wdraws_ref, mu = mu_test, dis = refmodel$dis,
-      cl_ref = seq_along(refmodel$wdraws_ref)
+      y_wobs_test = y_wobs_test,
+      data_aux_test = d_test$data %||% refmodel$fetch_data(),
+      family = refmodel$family, wdraws = refmodel$wdraws_ref, mu = mu_test,
+      dis = refmodel$dis, cl_ref = seq_along(refmodel$wdraws_ref)
     )
   }