From e9798763d10bc8aab41de12e52e45791494c4104 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 20 May 2024 17:28:41 +0300 Subject: [PATCH] Initial roxygen pass --- .gitignore | 5 + DESCRIPTION | 2 +- NAMESPACE | 38 ++++++++ R/doc-rstanarm-package.R | 2 +- R/jm_make_assoc_terms.R | 52 ++++++++--- R/log_lik.R | 195 +++++++++++++++++++++++++++++---------- R/loo-prediction.R | 64 +++++++------ R/misc.R | 128 +++++++++++++++++++------ R/posterior_survfit.R | 2 - R/posterior_traj.R | 2 - R/pp_data.R | 81 ++++++++++++---- R/predictive_error.R | 24 ++--- R/ps_check.R | 1 - R/stan_glm.fit.R | 26 ++++-- R/stanreg_list.R | 30 ++++-- 15 files changed, 480 insertions(+), 172 deletions(-) diff --git a/.gitignore b/.gitignore index 62f913f22..e61a23194 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,8 @@ vignettes/*.html vignettes/*_files .DS_Store revdep/* +*.o +*.cc +*.h +src/Makevars +R/stanmodels.R diff --git a/DESCRIPTION b/DESCRIPTION index 7313116f0..b67604600 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -79,4 +79,4 @@ UseLTO: true NeedsCompilation: yes URL: https://mc-stan.org/rstanarm/, https://discourse.mc-stan.org BugReports: https://github.com/stan-dev/rstanarm/issues -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index baf720ff3..c2c1fea8a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,7 +14,13 @@ S3method(as_draws_rvars,stanreg) S3method(bayes_R2,stanreg) S3method(coef,stanmvreg) S3method(coef,stanreg) +S3method(collapse_within_groups,default) +S3method(collapse_within_groups,matrix) S3method(confint,stanreg) +S3method(evaluate_log_survival,default) +S3method(evaluate_log_survival,matrix) +S3method(extract_pars,stanmvreg) +S3method(extract_pars,stansurv) S3method(family,stanmvreg) S3method(family,stanreg) S3method(fitted,stanmvreg) @@ -23,6 +29,8 @@ S3method(fixef,stanmvreg) S3method(fixef,stanreg) S3method(formula,stanmvreg) S3method(formula,stanreg) +S3method(get_model_data,stanmvreg) +S3method(get_model_data,stansurv) S3method(get_surv,stanjm) S3method(get_surv,stansurv) S3method(get_x,default) @@ -35,6 +43,15 @@ S3method(get_z,lmerMod) S3method(get_z,stanmvreg) S3method(kfold,stanreg) S3method(launch_shinystan,stanreg) +S3method(linear_predictor,default) +S3method(linear_predictor,matrix) +S3method(linkinv,character) +S3method(linkinv,family) +S3method(linkinv,stanmvreg) +S3method(linkinv,stanreg) +S3method(ll_args,stanjm) +S3method(ll_args,stanreg) +S3method(ll_args,stansurv) S3method(log_lik,stanjm) S3method(log_lik,stanmvreg) S3method(log_lik,stanreg) @@ -73,6 +90,7 @@ S3method(pp_check,stanreg) S3method(predict,stanreg) S3method(predictive_error,matrix) S3method(predictive_error,ppd) +S3method(predictive_error,stanmvreg) S3method(predictive_error,stanreg) S3method(predictive_interval,matrix) S3method(predictive_interval,ppd) @@ -87,18 +105,28 @@ S3method(print,summary.stanreg) S3method(print,survfit.stanjm) S3method(print,survfit.stansurv) S3method(prior_summary,stanreg) +S3method(psis,stanreg) +S3method(quadrature_sum,default) +S3method(quadrature_sum,matrix) S3method(ranef,stanmvreg) S3method(ranef,stanreg) +S3method(rename_loos,stanreg) +S3method(rename_loos,stanreg_list) S3method(residuals,stanmvreg) S3method(residuals,stanreg) S3method(se,stanmvreg) S3method(se,stanreg) S3method(sigma,stanmvreg) S3method(sigma,stanreg) +S3method(split2,matrix) +S3method(split2,vector) S3method(summary,stanmvreg) S3method(summary,stanreg) S3method(terms,stanmvreg) S3method(terms,stanreg) +S3method(truncate,numeric) +S3method(unpad_reTrms,array) +S3method(unpad_reTrms,default) S3method(update,stanjm) S3method(update,stanmvreg) S3method(update,stanreg) @@ -115,13 +143,17 @@ export(as_draws_matrix) export(as_draws_rvars) export(bayes_R2) export(cauchy) +export(collapse_within_groups) export(compare_models) export(decov) export(default_prior_coef) export(default_prior_intercept) export(dirichlet) +export(evaluate_log_survival) export(exponential) +export(extract_pars) export(fixef) +export(get_model_data) export(get_surv) export(get_x) export(get_y) @@ -133,7 +165,10 @@ export(kfold) export(laplace) export(lasso) export(launch_shinystan) +export(linear_predictor) +export(linkinv) export(lkj) +export(ll_args) export(log_lik) export(logit) export(loo) @@ -166,7 +201,9 @@ export(prior_options) export(prior_summary) export(product_normal) export(ps_check) +export(quadrature_sum) export(ranef) +export(rename_loos) export(se) export(sigma) export(stan_aov) @@ -196,6 +233,7 @@ export(stanmvreg_list) export(stanreg_list) export(student_t) export(tve) +export(unpad_reTrms) export(waic) if(getRversion()>='3.3.0') importFrom(stats, sigma) else importFrom(lme4,sigma) diff --git a/R/doc-rstanarm-package.R b/R/doc-rstanarm-package.R index a80773c9c..5dcc60d3f 100644 --- a/R/doc-rstanarm-package.R +++ b/R/doc-rstanarm-package.R @@ -104,4 +104,4 @@ #' @template reference-bayesvis #' @template reference-muth #' -NULL +"_PACKAGE" diff --git a/R/jm_make_assoc_terms.R b/R/jm_make_assoc_terms.R index be823b2f6..90955e698 100644 --- a/R/jm_make_assoc_terms.R +++ b/R/jm_make_assoc_terms.R @@ -376,21 +376,36 @@ get_element <- function(parts, m = 1, which = "eta", ...) { } } -# Collapse the linear predictor across the lower level units -# clustered an individual, using the function specified in the -# 'grp_assoc' argument -# -# @param eta The linear predictor evaluated for all lower level groups -# at the quadrature points. -# @param grp_idx An N*2 array providing the indices of the first (col 1) -# and last (col 2) observations in eta that correspond to individuals -# i = 1,...,N. -# @param grp_assoc Character string, the function to use to collapse -# across the lower level units clustered within individuals. -# @return A vector or matrix, depending on the method called. +#' Collapse the linear predictor across the lower level units +#' clustered an individual, using the function specified in the +#' 'grp_assoc' argument +#' +#' @param eta The linear predictor evaluated for all lower level groups +#' at the quadrature points. +#' @param grp_idx An N*2 array providing the indices of the first (col 1) +#' and last (col 2) observations in eta that correspond to individuals +#' i = 1,...,N. +#' @param grp_assoc Character string, the function to use to collapse +#' across the lower level units clustered within individuals. +#' @return A vector or matrix, depending on the method called. +#' @export collapse_within_groups <- function(eta, grp_idx, grp_assoc = "sum") { UseMethod("collapse_within_groups") } + +#' Collapse the linear predictor across the lower level units +#' clustered an individual, using the function specified in the +#' 'grp_assoc' argument +#' +#' @param eta The linear predictor evaluated for all lower level groups +#' at the quadrature points. +#' @param grp_idx An N*2 array providing the indices of the first (col 1) +#' and last (col 2) observations in eta that correspond to individuals +#' i = 1,...,N. +#' @param grp_assoc Character string, the function to use to collapse +#' across the lower level units clustered within individuals. +#' @return A vector or matrix, depending on the method called. +#' @export collapse_within_groups.default <- function(eta, grp_idx, grp_assoc) { N <- nrow(grp_idx) val <- rep(NA, N) @@ -400,6 +415,19 @@ collapse_within_groups.default <- function(eta, grp_idx, grp_assoc) { } val } +#' Collapse the linear predictor across the lower level units +#' clustered an individual, using the function specified in the +#' 'grp_assoc' argument +#' +#' @param eta The linear predictor evaluated for all lower level groups +#' at the quadrature points. +#' @param grp_idx An N*2 array providing the indices of the first (col 1) +#' and last (col 2) observations in eta that correspond to individuals +#' i = 1,...,N. +#' @param grp_assoc Character string, the function to use to collapse +#' across the lower level units clustered within individuals. +#' @return A vector or matrix, depending on the method called. +#' @export collapse_within_groups.matrix <- function(eta, grp_idx, grp_assoc) { N <- nrow(grp_idx) val <- matrix(NA, nrow = nrow(eta), ncol = N) diff --git a/R/log_lik.R b/R/log_lik.R index 960e902d6..049528cf5 100644 --- a/R/log_lik.R +++ b/R/log_lik.R @@ -201,25 +201,42 @@ ll_fun <- function(x, m = NULL) { get(fun, mode = "function") } -# get arguments needed for ll_fun -# @param object stanreg object -# @param newdata same as posterior predict -# @param offset vector of offsets (only required if model has offset term and +#' get arguments needed for ll_fun +#' @param object stanreg object +#' @param newdata same as posterior predict +#' @param offset vector of offsets (only required if model has offset term and # newdata is specified) -# @param m Integer specifying which submodel for stanmvreg objects -# @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold -# @param ... For models without group-specific terms (i.e., not stan_[g]lmer), -# if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to -# pass 'newx' and 'stanmat' from reloo or kfold (bypassing pp_data). This is a -# workaround in case there are issues with newdata containing factors with -# only a single level. Or for stanmvreg objects, then ... can be used to pass -# 'stanmat', which may be a matrix with a reduced number of draws (potentially -# just a single MCMC draw). -# @return a named list with elements data, draws, S (posterior sample size) and -# N = number of observations +#' @param m Integer specifying which submodel for stanmvreg objects +#' @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold +#' @param ... For models without group-specific terms (i.e., not stan_[g]lmer), +#' if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to +#' pass 'newx' and 'stanmat' from reloo or kfold (bypassing pp_data). This is a +#' workaround in case there are issues with newdata containing factors with +#' only a single level. Or for stanmvreg objects, then ... can be used to pass +#' 'stanmat', which may be a matrix with a reduced number of draws (potentially +#' just a single MCMC draw). +#' @return a named list with elements data, draws, S (posterior sample size) and +#' N = number of observations +#' @export ll_args <- function(object, ...) UseMethod("ll_args") -#--- ll_args for stanreg models +#' get arguments needed for ll_fun +#' @param object stanreg object +#' @param newdata same as posterior predict +#' @param offset vector of offsets (only required if model has offset term and +# newdata is specified) +#' @param m Integer specifying which submodel for stanmvreg objects +#' @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold +#' @param ... For models without group-specific terms (i.e., not stan_[g]lmer), +#' if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to +#' pass 'newx' and 'stanmat' from reloo or kfold (bypassing pp_data). This is a +#' workaround in case there are issues with newdata containing factors with +#' only a single level. Or for stanmvreg objects, then ... can be used to pass +#' 'stanmat', which may be a matrix with a reduced number of draws (potentially +#' just a single MCMC draw). +#' @return a named list with elements data, draws, S (posterior sample size) and +#' N = number of observations +#' @export ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, reloo_or_kfold = FALSE, ...) { validate_stanreg_object(object) @@ -391,7 +408,19 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, return(out) } -#--- ll_args for stansurv models +#' get arguments needed for ll_fun +#' @param object stansurv object +#' @param newdata same as posterior predict +#' @param ... For models without group-specific terms (i.e., not stan_[g]lmer), +#' if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to +#' pass 'newx' and 'stanmat' from reloo or kfold (bypassing pp_data). This is a +#' workaround in case there are issues with newdata containing factors with +#' only a single level. Or for stanmvreg objects, then ... can be used to pass +#' 'stanmat', which may be a matrix with a reduced number of draws (potentially +#' just a single MCMC draw). +#' @return a named list with elements data, draws, S (posterior sample size) and +#' N = number of observations +#' @export ll_args.stansurv <- function(object, newdata = NULL, ...) { validate_stansurv_object(object) @@ -754,16 +783,17 @@ ll_args.stansurv <- function(object, newdata = NULL, ...) { # log-likelihood functions for stanjm objects only ---------------------- -# Alternative ll_args method for stanjm objects that allows data and pars to be -# passed directly, rather than constructed using pp_data within the ll_args -# method. This can be much faster when used in the MH algorithm within -# posterior_survfit, since it doesn't require repeated calls to pp_data. -# -# @param object A stanmvreg object -# @param data Output from .pp_data_jm -# @param pars Output from extract_pars -# @param m Integer specifying which submodel -# @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold +#' Alternative ll_args method for stanjm objects that allows data and pars to be +#' passed directly, rather than constructed using pp_data within the ll_args +#' method. This can be much faster when used in the MH algorithm within +#' posterior_survfit, since it doesn't require repeated calls to pp_data. +#' +#' @param object A stanmvreg object +#' @param data Output from .pp_data_jm +#' @param pars Output from extract_pars +#' @param m Integer specifying which submodel +#' @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold +#' @export ll_args.stanjm <- function(object, data, pars, m = 1, reloo_or_kfold = FALSE, ...) { validate_stanjm_object(object) @@ -1102,24 +1132,39 @@ evaluate_log_basehaz2 <- function(times, basehaz, coefs) { log_basehaz } -# Evaluate the log baseline hazard at the specified times -# given the vector or matrix of MCMC draws for the baseline -# hazard coeffients / parameters -# -# @param log_haz A vector containing the log hazard for each -# individual, evaluated at each of the quadrature points. The -# vector should be ordered such that the first N elements contain -# the log_haz evaluated for each individual at quadrature point 1, -# then the next N elements are the log_haz evaluated for each -# individual at quadrature point 2, and so on. -# @param qnodes Integer specifying the number of quadrature nodes -# at which the log hazard was evaluated for each individual. -# @param qwts A vector of unstandardised GK quadrature weights. -# @return A vector or matrix of log survival probabilities. +#' Evaluate the log baseline hazard at the specified times +#' given the vector or matrix of MCMC draws for the baseline +#' hazard coeffients / parameters +#' +#' @param log_haz A vector containing the log hazard for each +#' individual, evaluated at each of the quadrature points. The +#' vector should be ordered such that the first N elements contain +#' the log_haz evaluated for each individual at quadrature point 1, +#' then the next N elements are the log_haz evaluated for each +#' individual at quadrature point 2, and so on. +#' @param qnodes Integer specifying the number of quadrature nodes +#' at which the log hazard was evaluated for each individual. +#' @param qwts A vector of unstandardised GK quadrature weights. +#' @return A vector or matrix of log survival probabilities. +#' @export evaluate_log_survival <- function(log_haz, qnodes, qwts) { UseMethod("evaluate_log_survival") } - +#' Evaluate the log baseline hazard at the specified times +#' given the vector or matrix of MCMC draws for the baseline +#' hazard coeffients / parameters +#' +#' @param log_haz A vector containing the log hazard for each +#' individual, evaluated at each of the quadrature points. The +#' vector should be ordered such that the first N elements contain +#' the log_haz evaluated for each individual at quadrature point 1, +#' then the next N elements are the log_haz evaluated for each +#' individual at quadrature point 2, and so on. +#' @param qnodes Integer specifying the number of quadrature nodes +#' at which the log hazard was evaluated for each individual. +#' @param qwts A vector of unstandardised GK quadrature weights. +#' @return A vector or matrix of log survival probabilities. +#' @export evaluate_log_survival.default <- function(log_haz, qnodes, qwts) { # convert log hazard to hazard haz <- exp(log_haz) @@ -1132,6 +1177,21 @@ evaluate_log_survival.default <- function(log_haz, qnodes, qwts) { -cum_haz } +#' Evaluate the log baseline hazard at the specified times +#' given the vector or matrix of MCMC draws for the baseline +#' hazard coeffients / parameters +#' +#' @param log_haz A vector containing the log hazard for each +#' individual, evaluated at each of the quadrature points. The +#' vector should be ordered such that the first N elements contain +#' the log_haz evaluated for each individual at quadrature point 1, +#' then the next N elements are the log_haz evaluated for each +#' individual at quadrature point 2, and so on. +#' @param qnodes Integer specifying the number of quadrature nodes +#' at which the log hazard was evaluated for each individual. +#' @param qwts A vector of unstandardised GK quadrature weights. +#' @return A vector or matrix of log survival probabilities. +#' @export evaluate_log_survival.matrix <- function(log_haz, qnodes, qwts) { # convert log hazard to hazard haz <- exp(log_haz) @@ -1273,34 +1333,60 @@ evaluate_log_surv <- function(times, basehaz, betas, b = NULL, aux, #--------------- +#' Apply quadrature weights and sum over nodes +#' +#' @param x Inputs +#' @param qnodes Nodes +#' @param qwts Quadrature weights +#' @export quadrature_sum <- function(x, qnodes, qwts) { UseMethod("quadrature_sum") } +#' Apply quadrature weights and sum over nodes +#' +#' @param x Inputs +#' @param qnodes Nodes +#' @param qwts Quadrature weights +#' @export quadrature_sum.default <- function(x, qnodes, qwts) { weighted_x <- qwts * x # apply quadrature weights splitted_x <- split_vector(x, n_segments = qnodes) # split at each quad node Reduce('+', splitted_x) # sum over the quad nodes } +#' Apply quadrature weights and sum over nodes +#' +#' @param x Inputs +#' @param qnodes Nodes +#' @param qwts Quadrature weights +#' @export quadrature_sum.matrix <- function(x, qnodes, qwts) { weighted_x <- sweep_multiply(x, qwts, margin = 2L) # apply quadrature weights splitted_x <- array2list(weighted_x, nsplits = qnodes) # split at each quad node Reduce('+', splitted_x) # sum over the quad nodes } -# Split a vector or matrix into a specified number of segments and return -# each segment as an element of a list. The matrix method allows splitting -# across the column (bycol = TRUE) or row margin (bycol = FALSE). -# -# @param x A vector or matrix. -# @param n_segments Integer specifying the number of segments. -# @param bycol Logical, should a matrix be split along the column or row margin? -# @return A list with n_segments elements. +#' Split a vector or matrix into a specified number of segments and return +#' each segment as an element of a list. The matrix method allows splitting +#' across the column (bycol = TRUE) or row margin (bycol = FALSE). +#' +#' @param x A vector or matrix. +#' @param n_segments Integer specifying the number of segments. +#' @param bycol Logical, should a matrix be split along the column or row margin? +#' @return A list with n_segments elements. split2 <- function(x, n_segments = 1, ...) { UseMethod("split2") } +#' Split a vector or matrix into a specified number of segments and return +#' each segment as an element of a list. The matrix method allows splitting +#' across the column (bycol = TRUE) or row margin (bycol = FALSE). +#' +#' @param x A vector +#' @param n_segments Integer specifying the number of segments. +#' @param bycol Logical, should a matrix be split along the column or row margin? +#' @export split2.vector <- function(x, n_segments = 1, ...) { len <- length(x) segment_length <- len %/% n_segments @@ -1308,7 +1394,14 @@ split2.vector <- function(x, n_segments = 1, ...) { stop("Dividing x by n_segments does not result in an integer.") split(x, rep(1:n_segments, each = segment_length)) } - +#' Split a vector or matrix into a specified number of segments and return +#' each segment as an element of a list. The matrix method allows splitting +#' across the column (bycol = TRUE) or row margin (bycol = FALSE). +#' +#' @param x A matrix +#' @param n_segments Integer specifying the number of segments. +#' @param bycol Logical, should a matrix be split along the column or row margin? +#' @export split2.matrix <- function(x, n_segments = 1, bycol = TRUE) { len <- if (bycol) ncol(x) else nrow(x) segment_length <- len %/% n_segments diff --git a/R/loo-prediction.R b/R/loo-prediction.R index e2c19d38c..ea840c846 100644 --- a/R/loo-prediction.R +++ b/R/loo-prediction.R @@ -1,68 +1,68 @@ #' Compute weighted expectations using LOO -#' +#' #' These functions are wrappers around the \code{\link[loo]{E_loo}} function #' (\pkg{loo} package) that provide compatibility for \pkg{rstanarm} models. -#' +#' #' @export #' @aliases loo_predict loo_linpred loo_predictive_interval -#' +#' #' @template reference-loo #' @template reference-bayesvis #' @templateVar stanregArg object #' @template args-stanreg-object -#' @param psis_object An object returned by \code{\link[loo]{psis}}. If missing +#' @param psis_object An object returned by \code{\link[loo]{psis}}. If missing #' then \code{psis} will be run internally, which may be time consuming #' for models fit to very large datasets. #' @param ... Currently unused. #' @inheritParams loo::E_loo -#' -#' @return A list with elements \code{value} and \code{pareto_k}. -#' -#' For \code{loo_predict} and \code{loo_linpred} the value component is a -#' vector with one element per observation. -#' +#' +#' @return A list with elements \code{value} and \code{pareto_k}. +#' +#' For \code{loo_predict} and \code{loo_linpred} the value component is a +#' vector with one element per observation. +#' #' For \code{loo_predictive_interval} the \code{value} component is a matrix #' with one row per observation and two columns (like #' \code{\link{predictive_interval}}). \code{loo_predictive_interval(..., prob #' = p)} is equivalent to \code{loo_predict(..., type = "quantile", probs = #' c(a, 1-a))} with \code{a = (1 - p)/2}, except it transposes the result and #' adds informative column names. -#' +#' #' See \code{\link[loo]{E_loo}} and \code{\link[loo]{pareto-k-diagnostic}} for #' details on the \code{pareto_k} diagnostic. -#' +#' #' @examples #' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \dontrun{ #' if (!exists("example_model")) example(example_model) -#' +#' #' # optionally, log-weights can be pre-computed and reused #' psis_result <- loo::psis(log_ratios = -log_lik(example_model)) -#' +#' #' loo_probs <- loo_linpred(example_model, type = "mean", transform = TRUE, psis_object = psis_result) #' str(loo_probs) -#' +#' #' loo_pred_var <- loo_predict(example_model, type = "var", psis_object = psis_result) #' str(loo_pred_var) -#' +#' #' loo_pred_ints <- loo_predictive_interval(example_model, prob = 0.8, psis_object = psis_result) #' str(loo_pred_ints) #' } #' } loo_predict.stanreg <- - function(object, - type = c("mean", "var", "quantile"), + function(object, + type = c("mean", "var", "quantile"), probs = 0.5, ..., psis_object = NULL) { - + if ("lw" %in% names(list(...))) { stop( - "Due to changes in the 'loo' package, the 'lw' argument ", + "Due to changes in the 'loo' package, the 'lw' argument ", "is no longer supported. Use the 'psis_object' argument instead." ) } - + type <- match.arg(type) log_ratios <- -log_lik(object) if (is.null(psis_object)) { @@ -70,12 +70,12 @@ loo_predict.stanreg <- r_eff <- loo::relative_eff(exp(-log_ratios), chain_id = chain_id_for_loo(object)) psis_object <- loo::psis(log_ratios, r_eff = r_eff) } - + preds <- posterior_predict(object) if (is_polr(object) && !is_scobit(object)) { preds <- polr_yrep_to_numeric(preds) } - + loo::E_loo( x = preds, psis_object = psis_object, @@ -88,22 +88,22 @@ loo_predict.stanreg <- #' @rdname loo_predict.stanreg #' @export #' @param transform Passed to \code{\link{posterior_linpred}}. -#' +#' loo_linpred.stanreg <- function(object, type = c("mean", "var", "quantile"), probs = 0.5, transform = FALSE, - ..., + ..., psis_object = NULL) { - + if ("lw" %in% names(list(...))) { stop( - "Due to changes in the 'loo' package, the 'lw' argument ", + "Due to changes in the 'loo' package, the 'lw' argument ", "is no longer supported. Use the 'psis_object' argument instead." ) } - + type <- match.arg(type) log_ratios <- -log_lik(object) if (is.null(psis_object)) { @@ -111,10 +111,10 @@ loo_linpred.stanreg <- r_eff <- loo::relative_eff(exp(-log_ratios), chain_id = chain_id_for_loo(object)) psis_object <- loo::psis(log_ratios, r_eff = r_eff) } - + type <- match.arg(type) linpreds <- posterior_linpred(object, transform = transform) - + loo::E_loo( x = linpreds, psis_object = psis_object, @@ -155,6 +155,10 @@ loo_predictive_interval.stanreg <- # internal ---------------------------------------------------------------- +#' Perform importance sampling with stanreg objects +#' +#' @export +#' @param log_ratios log_ratios for PSIS computation psis.stanreg <- function(log_ratios, ...) { object <- log_ratios message("Running PSIS to compute weights...") diff --git a/R/misc.R b/R/misc.R index 5e61fb7f5..b64c8c627 100644 --- a/R/misc.R +++ b/R/misc.R @@ -574,18 +574,30 @@ set_prior_scale <- function(scale, default, link) { } -# Methods for creating linear predictor -# -# Make linear predictor vector from x and point estimates for beta, or linear -# predictor matrix from x and full posterior sample of beta. -# -# @param beta A vector or matrix or parameter estimates. -# @param x Predictor matrix. -# @param offset Optional offset vector. -# @return A vector or matrix. +#' Methods for creating linear predictor +#' +#' Make linear predictor vector from x and point estimates for beta, or linear +#' predictor matrix from x and full posterior sample of beta. +#' +#' @param beta A vector or matrix or parameter estimates. +#' @param x Predictor matrix. +#' @param offset Optional offset vector. +#' @return A vector or matrix. +#' @export linear_predictor <- function(beta, x, offset = NULL) { UseMethod("linear_predictor") } + +#' Methods for creating linear predictor +#' +#' Make linear predictor vector from x and point estimates for beta, or linear +#' predictor matrix from x and full posterior sample of beta. +#' +#' @param beta A vector or matrix or parameter estimates. +#' @param x Predictor matrix. +#' @param offset Optional offset vector. +#' @return A vector or matrix. +#' @export linear_predictor.default <- function(beta, x, offset = NULL) { eta <- as.vector(if (NCOL(x) == 1L) x * beta else x %*% beta) if (length(offset)) @@ -593,6 +605,17 @@ linear_predictor.default <- function(beta, x, offset = NULL) { return(eta) } + +#' Methods for creating linear predictor +#' +#' Make linear predictor vector from x and point estimates for beta, or linear +#' predictor matrix from x and full posterior sample of beta. +#' +#' @param beta A vector or matrix or parameter estimates. +#' @param x Predictor matrix. +#' @param offset Optional offset vector. +#' @return A vector or matrix. +#' @export linear_predictor.matrix <- function(beta, x, offset = NULL) { if (NCOL(beta) == 1L) beta <- as.matrix(beta) @@ -680,24 +703,57 @@ get_surv.stanjm <- function(object, ...) { object$survmod$mod$y %ORifNULL% stop("response not found") } -# Get inverse link function -# -# @param x A stanreg object, family object, or string. -# @param ... Other arguments passed to methods. For a \code{stanmvreg} object -# this can be an integer \code{m} specifying the submodel. -# @return The inverse link function associated with x. +#' Get inverse link function +#' +#' @param x A stanreg object, family object, or string. +#' @param ... Other arguments passed to methods. For a \code{stanmvreg} object +#' this can be an integer \code{m} specifying the submodel. +#' @return The inverse link function associated with x. +#' @export linkinv <- function(x, ...) UseMethod("linkinv") + +#' Get inverse link function +#' +#' @param x A stanreg object +#' @param ... Other arguments passed to methods. For a \code{stanmvreg} object +#' this can be an integer \code{m} specifying the submodel. +#' @return The inverse link function associated with x. +#' @export linkinv.stanreg <- function(x, ...) { if (is(x, "polr")) polr_linkinv(x) else family(x)$linkinv } + +#' Get inverse link function +#' +#' @param x A stanmvreg object +#' @param ... Other arguments passed to methods. For a \code{stanmvreg} object +#' this can be an integer \code{m} specifying the submodel. +#' @return The inverse link function associated with x. +#' @export linkinv.stanmvreg <- function(x, m = NULL, ...) { ret <- lapply(family(x), `[[`, "linkinv") stub <- get_stub(x) if (!is.null(m)) ret[[m]] else list_nms(ret, stub = stub) } + +#' Get inverse link function +#' +#' @param x A family object +#' @param ... Other arguments passed to methods. For a \code{stanmvreg} object +#' this can be an integer \code{m} specifying the submodel. +#' @return The inverse link function associated with x. +#' @export linkinv.family <- function(x, ...) { x$linkinv } + +#' Get inverse link function +#' +#' @param x A character string +#' @param ... Other arguments passed to methods. For a \code{stanmvreg} object +#' this can be an integer \code{m} specifying the submodel. +#' @return The inverse link function associated with x. +#' @export linkinv.character <- function(x, ...) { stopifnot(length(x) == 1) polr_linkinv(x) @@ -1549,16 +1605,24 @@ get_time_seq <- function(increments, t0, t1, simplify = TRUE) { return(val) } -# Extract parameters from stanmat and return as a list -# -# @param object A stanmvreg or stansurv object -# @param stanmat A matrix of posterior draws, may be provided if the desired -# stanmat is only a subset of the draws from as.matrix(object$stanfit) -# @return A named list +#' Extract parameters from stanmat and return as a list +#' +#' @param object A stanmvreg or stansurv object +#' @param stanmat A matrix of posterior draws, may be provided if the desired +#' stanmat is only a subset of the draws from as.matrix(object$stanfit) +#' @return A named list +#' @export extract_pars <- function(object, ...) { UseMethod("extract_pars") } +#' Extract parameters from stanmat and return as a list +#' +#' @param object A stanmvreg or stansurv object +#' @param stanmat A matrix of posterior draws, may be provided if the desired +#' stanmat is only a subset of the draws from as.matrix(object$stanfit) +#' @return A named list +#' @export extract_pars.stansurv <- function(object, stanmat = NULL, means = FALSE) { validate_stansurv_object(object) if (is.null(stanmat)) @@ -1580,6 +1644,14 @@ extract_pars.stansurv <- function(object, stanmat = NULL, means = FALSE) { nlist(alpha, beta, beta_tve, aux, smooth, b, stanmat) } + +#' Extract parameters from stanmat and return as a list +#' +#' @param object A stanmvreg or stansurv object +#' @param stanmat A matrix of posterior draws, may be provided if the desired +#' stanmat is only a subset of the draws from as.matrix(object$stanfit) +#' @return A named list +#' @export extract_pars.stanmvreg <- function(object, stanmat = NULL, means = FALSE) { validate_stanmvreg_object(object) M <- get_M(object) @@ -2131,12 +2203,14 @@ sample_stanmat <- function(object, draws = NULL, default_draws = NA) { stanmat } -# Method to truncate a numeric vector at defined limits -# -# @param con A numeric vector. -# @param lower Scalar, the lower limit for the returned vector. -# @param upper Scalar, the upper limit for the returned vector. -# @return A numeric vector. +#' Method to truncate a numeric vector at defined limits +#' +#' @param con A numeric vector. +#' @param lower Scalar, the lower limit for the returned vector. +#' @param upper Scalar, the upper limit for the returned vector. +#' +#' @return A numeric vector. +#' @export truncate.numeric <- function(con, lower = NULL, upper = NULL) { if (!is.null(lower)) con[con < lower] <- lower if (!is.null(upper)) con[con > upper] <- upper diff --git a/R/posterior_survfit.R b/R/posterior_survfit.R index a5cb10eef..f0511062e 100644 --- a/R/posterior_survfit.R +++ b/R/posterior_survfit.R @@ -297,7 +297,6 @@ #' \strong{67}, 819. #' #' @examples -#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ #' # Run example model if not already loaded #' if (!exists("example_jm")) example(example_jm) @@ -1129,7 +1128,6 @@ print.survfit.stanjm <- function(x, digits = 4, ...) { #' \code{\link{posterior_traj}}, \code{\link{plot.predict.stanjm}} #' #' @examples -#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ #' # Run example model if not already loaded #' if (!exists("example_jm")) example(example_jm) diff --git a/R/posterior_traj.R b/R/posterior_traj.R index 0c89d9887..71d0a3847 100644 --- a/R/posterior_traj.R +++ b/R/posterior_traj.R @@ -181,7 +181,6 @@ #' \code{\link{posterior_survfit}} #' #' @examples -#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ #' # Run example model if not already loaded #' if (!exists("example_jm")) example(example_jm) @@ -553,7 +552,6 @@ posterior_traj <- function(object, #' \code{\link{posterior_survfit}}, \code{\link{plot.survfit.stanjm}} #' #' @examples -#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ #' # Run example model if not already loaded #' if (!exists("example_jm")) example(example_jm) diff --git a/R/pp_data.R b/R/pp_data.R index 126c5357a..36400c31c 100644 --- a/R/pp_data.R +++ b/R/pp_data.R @@ -560,28 +560,48 @@ pp_data <- return(res) } -# Return a data frame for each submodel that: -# (1) only includes variables used in the model formula -# (2) only includes rows contained in the glmod/coxmod model frames -# (3) ensures that additional variables that are required -# such as the ID variable or variables used in the -# interaction-type association structures, are included. -# -# It is necessary to drop unneeded variables though so that -# errors are not encountered if the original data contained -# NA values for variables unrelated to the model formula. -# We generate a data frame here for in-sample predictions -# rather than using a model frame, since some quantities will -# need to be recalculated at quadrature points etc, for example -# in posterior_survfit. -# -# @param object A stansurv, stanmvreg or stanjm object. -# @param m Integer specifying which submodel to get the -# prediction data frame for (for stanmvreg or stanjm objects). -# @return A data frame or list of data frames with all the -# (unevaluated) variables required for predictions. +#' Return a data frame for each submodel that: +#' (1) only includes variables used in the model formula +#' (2) only includes rows contained in the glmod/coxmod model frames +#' (3) ensures that additional variables that are required +#' such as the ID variable or variables used in the +#' interaction-type association structures, are included. +#' +#' It is necessary to drop unneeded variables though so that +#' errors are not encountered if the original data contained +#' NA values for variables unrelated to the model formula. +#' We generate a data frame here for in-sample predictions +#' rather than using a model frame, since some quantities will +#' need to be recalculated at quadrature points etc, for example +#' in posterior_survfit. +#' +#' @param object A stansurv, stanmvreg or stanjm object. +#' @param m Integer specifying which submodel to get the +#' prediction data frame for (for stanmvreg or stanjm objects). +#' @return A data frame or list of data frames with all the +#' (unevaluated) variables required for predictions. +#' @export get_model_data <- function(object, ...) UseMethod("get_model_data") +#' Return a data frame for each submodel that: +#' (1) only includes variables used in the model formula +#' (2) only includes rows contained in the glmod/coxmod model frames +#' (3) ensures that additional variables that are required +#' such as the ID variable or variables used in the +#' interaction-type association structures, are included. +#' +#' It is necessary to drop unneeded variables though so that +#' errors are not encountered if the original data contained +#' NA values for variables unrelated to the model formula. +#' We generate a data frame here for in-sample predictions +#' rather than using a model frame, since some quantities will +#' need to be recalculated at quadrature points etc, for example +#' in posterior_survfit. +#' +#' @param object A stansurv object. +#' @return A data frame or list of data frames with all the +#' (unevaluated) variables required for predictions. +#' @export get_model_data.stansurv <- function(object, ...) { validate_stansurv_object(object) terms <- terms(object, fixed.only = FALSE) @@ -589,6 +609,27 @@ get_model_data.stansurv <- function(object, ...) { get_all_vars(terms, object$data)[row_nms, , drop = FALSE] } +#' Return a data frame for each submodel that: +#' (1) only includes variables used in the model formula +#' (2) only includes rows contained in the glmod/coxmod model frames +#' (3) ensures that additional variables that are required +#' such as the ID variable or variables used in the +#' interaction-type association structures, are included. +#' +#' It is necessary to drop unneeded variables though so that +#' errors are not encountered if the original data contained +#' NA values for variables unrelated to the model formula. +#' We generate a data frame here for in-sample predictions +#' rather than using a model frame, since some quantities will +#' need to be recalculated at quadrature points etc, for example +#' in posterior_survfit. +#' +#' @param object A stanmvreg. +#' @param m Integer specifying which submodel to get the +#' prediction data frame for (for stanmvreg or stanjm objects). +#' @return A data frame or list of data frames with all the +#' (unevaluated) variables required for predictions. +#' @export get_model_data.stanmvreg <- function(object, m = NULL, ...) { validate_stanmvreg_object(object) diff --git a/R/predictive_error.R b/R/predictive_error.R index ba9c34364..6cada3ce9 100644 --- a/R/predictive_error.R +++ b/R/predictive_error.R @@ -129,18 +129,18 @@ predictive_error.ppd <- function(object, y, ...) { predictive_error(unclass(object), y = y, ...) } -# @rdname predictive_error.stanreg -# @export -# @param m For \code{stanmvreg} models, the submodel for which to calculate -# the prediction error. Can be an integer, or for \code{\link{stan_mvmer}} -# models it can be \code{"y1"}, \code{"y2"}, etc, or for \code{\link{stan_jm}} -# models it can be \code{"Event"}, \code{"Long1"}, \code{"Long2"}, etc. -# @param t,u Only relevant for \code{\link{stan_jm}} models and when \code{m = "Event"}. -# The argument \code{t} specifies the time up to which individuals must have survived -# as well as being the time up to which the longitudinal data in \code{newdata} -# is available. The argument \code{u} specifies the time at which the -# prediction error should be calculated (i.e. the time horizon). -# +#' @rdname predictive_error.stanreg +#' @export +#' @param m For \code{stanmvreg} models, the submodel for which to calculate +#' the prediction error. Can be an integer, or for \code{\link{stan_mvmer}} +#' models it can be \code{"y1"}, \code{"y2"}, etc, or for \code{\link{stan_jm}} +#' models it can be \code{"Event"}, \code{"Long1"}, \code{"Long2"}, etc. +#' @param t,u Only relevant for \code{\link{stan_jm}} models and when \code{m = "Event"}. +#' The argument \code{t} specifies the time up to which individuals must have survived +#' as well as being the time up to which the longitudinal data in \code{newdata} +#' is available. The argument \code{u} specifies the time at which the +#' prediction error should be calculated (i.e. the time horizon). +#' predictive_error.stanmvreg <- function(object, newdataLong = NULL, diff --git a/R/ps_check.R b/R/ps_check.R index d3d8362f9..2ce6448a2 100644 --- a/R/ps_check.R +++ b/R/ps_check.R @@ -69,7 +69,6 @@ #' (for \code{\link{stan_jm}} models only) #' #' @examples -#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ #' if (!exists("example_jm")) example(example_jm) #' # Compare estimated survival function to Kaplan-Meier curve diff --git a/R/stan_glm.fit.R b/R/stan_glm.fit.R index 4747a2e3b..282461dff 100644 --- a/R/stan_glm.fit.R +++ b/R/stan_glm.fit.R @@ -896,13 +896,20 @@ pad_reTrms <- function(Ztlist, cnms, flist) { return(nlist(Z, cnms, flist)) } -# Drop the extra reTrms from a matrix x -# -# @param x A matrix or array (e.g. the posterior sample or matrix of summary -# stats) -# @param columns Do the columns (TRUE) or rows (FALSE) correspond to the -# variables? +#' Drop the extra reTrms from a matrix x +#' +#' @param x A matrix or array (e.g. the posterior sample or matrix of summary +#' stats) +#' @param columns Do the columns (TRUE) or rows (FALSE) correspond to the +#' variables? +#' @export unpad_reTrms <- function(x, ...) UseMethod("unpad_reTrms") + +#' Drop the extra reTrms from a matrix x +#' +#' @param x A matrix (e.g. the posterior sample or matrix of summary +#' stats) +#' @export unpad_reTrms.default <- function(x, ...) { if (is.matrix(x) || is.array(x)) return(unpad_reTrms.array(x, ...)) @@ -910,6 +917,13 @@ unpad_reTrms.default <- function(x, ...) { x[keep] } +#' Drop the extra reTrms from an array x +#' +#' @param x An array (e.g. the posterior sample or matrix of summary +#' stats) +#' @param columns Do the columns (TRUE) or rows (FALSE) correspond to the +#' variables? +#' @export unpad_reTrms.array <- function(x, columns = TRUE, ...) { ndim <- length(dim(x)) if (ndim > 3) diff --git a/R/stanreg_list.R b/R/stanreg_list.R index b4044c36a..b1b9c5944 100644 --- a/R/stanreg_list.R +++ b/R/stanreg_list.R @@ -273,13 +273,24 @@ stanreg_list_families <- function(mods) { -# loo/waic/kfold objects created by rstanarm have a model_name attribute. -# when a stanreg_list is created those attributes should be changed to match -# the names of the models used for the stanreg_list in case user has specified -# the model_names argument +#' loo/waic/kfold objects created by rstanarm have a model_name attribute. +#' when a stanreg_list is created those attributes should be changed to match +#' the names of the models used for the stanreg_list in case user has specified +#' the model_names argument +#' +#' @param loo/waic/kfold object +#' @return Renamed object +#' +#' @export rename_loos <- function(x,...) UseMethod("rename_loos") -# Change model_name attributes of a loo/waic/kfold object stored in a stanreg object, +#' Change model_name attributes of a loo/waic/kfold object stored in a stanreg object, +#' +#' @param loo/waic/kfold object +#' @param new_model_name new name for object +#' @return Renamed object +#' +#' @export rename_loos.stanreg <- function(x, new_model_name,...) { for (criterion in c("loo", "waic", "kfold")) { if (!is.null(x[[criterion]])) { @@ -289,8 +300,13 @@ rename_loos.stanreg <- function(x, new_model_name,...) { return(x) } -# Change model_name attributes of loo/waic/kfold objects to correspond to -# model names used for stanreg_list +#' Change model_name attributes of loo/waic/kfold objects to correspond to +#' model names used for stanreg_list +#' +#' @param loo/waic/kfold object +#' @return Renamed object +#' +#' @export rename_loos.stanreg_list <- function(x, ...) { for (j in seq_along(x)) { x[[j]] <- rename_loos.stanreg(x[[j]], new_model_name = names(x)[j])