diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 6158e7167..fbb81b113 100644 --- a/R/cv_varsel.R +++ b/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, diff --git a/R/extend_family.R b/R/extend_family.R index 89d570bf3..d77a08a31 100644 --- a/R/extend_family.R +++ b/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,10 +167,27 @@ #' 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 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. +#' #' 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. @@ -178,7 +195,8 @@ #' 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: diff --git a/R/formula.R b/R/formula.R index 479745d54..439ec6c09 100644 --- a/R/formula.R +++ b/R/formula.R @@ -742,6 +742,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 diff --git a/R/latent.R b/R/latent.R index 23d4519e4..0d0ee0b90 100644 --- a/R/latent.R +++ b/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,7 +61,8 @@ 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) @@ -67,7 +70,7 @@ latent_ll_oscale_poiss <- function(ilpreds, } 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) { diff --git a/R/methods.R b/R/methods.R index d0758d912..07691dfda 100644 --- a/R/methods.R +++ b/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)) { diff --git a/R/misc.R b/R/misc.R index 356aa59da..28ee8b9c9 100644 --- a/R/misc.R +++ b/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)) } diff --git a/R/projfun.R b/R/projfun.R index 930972ab9..f09170620 100644 --- a/R/projfun.R +++ b/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, diff --git a/R/refmodel.R b/R/refmodel.R index 17286999d..51e4a55f5 100644 --- a/R/refmodel.R +++ b/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)) ) diff --git a/R/summary_funs.R b/R/summary_funs.R index d67d6b3e7..f8b5ad62b 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -3,10 +3,13 @@ # values). # # @param y_wobs_test A `list` (but we encourage to use a `data.frame`), at least -# with elements (columns) `y` (response values) and `wobs` (observation +# with elements (columns) `y` (response values), and `wobs` (observation # 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`.") diff --git a/R/varsel.R b/R/varsel.R index d9cc50a06..6b71cb7e4 100644 --- a/R/varsel.R +++ b/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) ) } diff --git a/man/extend_family.Rd b/man/extend_family.Rd index a376044bb..6cd27ec15 100644 --- a/man/extend_family.Rd +++ b/man/extend_family.Rd @@ -193,8 +193,8 @@ given in \code{family$cats}, after taking \code{latent_y_unqs} into account). The function supplied to argument \code{latent_ll_oscale} needs to have the prototype -\if{html}{\out{