R/shrinkDSM.R

Defines functions shrinkDSM

Documented in shrinkDSM

#' Markov Chain Monte Carlo (MCMC) for time-varying parameter survival models with shrinkage
#'
#' \code{shrinkDSM} samples from the joint posterior distribution of the parameters of a time-varying
#' parameter survival model with shrinkage and returns the MCMC draws.
#' See also \code{\link{shrinkTVP}} to see more examples of how to modify the prior setup of the time-varying
#' component of the model.
#'
#' @param formula an object of class "formula": a symbolic representation of the model, as in the
#' function \code{lm}. For details, see \code{\link{formula}}.
#' @param data \emph{optional} data frame containing the response variable and the covariates. If not found in \code{data},
#' the variables are taken from \code{environment(formula)}, typically the environment from which \code{shrinkDSM}
#' is called. No \code{NA}s are allowed in the response variable and the covariates.
#' @param mod_type character string that reads either \code{"triple"}, \code{"double"} or \code{"ridge"}.
#' Determines whether the triple gamma, double gamma or ridge prior are used for \code{theta_sr} and \code{beta_mean}.
#' The default is "double".
#' @param delta The status indicator of the last period, 0=censored or 1=event observed.
#' @param S integer vector of time points that start a new interval.
#' Parameters are fixed within an interval and vary across intervals.
#' @param group \emph{optional} grouping indicator for latent factor.
#' @param subset \emph{optional} vector specifying a subset of observations to be used in the fitting process.
#' @param niter positive integer, indicating the number of MCMC iterations
#' to perform, including the burn-in. Has to be larger than or equal to \code{nburn} + 2. The default value is 10000.
#' @param nburn non-negative integer, indicating the number of iterations discarded
#' as burn-in. Has to be smaller than or equal to \code{niter} - 2. The default value is \code{round(niter / 2)}.
#' @param nthin positive integer, indicating the degree of thinning to be performed. Every \code{nthin} draw is kept and returned.
#' The default value is 1, implying that every draw is kept.
#' @param learn_a_xi  logical value indicating whether to learn a_xi, the spike parameter of the state variances.
#' The default value is \code{TRUE}.
#' @param learn_a_tau logical value indicating whether to learn a_tau, the spike parameter of the mean of the
#' initial values of the states. The default value is \code{TRUE}.
#' @param a_xi positive, real number, indicating the (fixed) value for a_xi. Ignored if
#' \code{learn_a_xi} is \code{TRUE} or \code{mod_type} is set to \code{"ridge"}. The default value is 0.1.
#' @param a_tau positive, real number, indicating the (fixed) value for a_tau. Ignored if
#' \code{learn_a_tau} is \code{TRUE} or \code{mod_type} is set to \code{"ridge"}. The default value is 0.1.
#' @param learn_c_xi logical value indicating whether to learn c_xi, the tail parameter of the state variances.
#' Ignored if \code{mod_type} is not set to \code{"triple"} or \code{a_eq_c_xi} is set to \code{TRUE}.
#' The default value is \code{TRUE}.
#' @param learn_c_tau logical value indicating whether to learn c_tau, the tail parameter of the mean of the
#' initial values of the states. Ignored if \code{mod_type} is not set to \code{"triple"} or \code{a_eq_c_tau} is set to \code{TRUE}.
#' The default value is \code{TRUE}.
#' @param c_xi positive, real number, indicating the (fixed) value for c_xi. Ignored if
#' \code{learn_c_xi} is \code{TRUE}, \code{mod_type} is not set to \code{"triple"} or \code{a_eq_c_xi} is set to \code{TRUE}.
#' The default value is 0.1.
#' @param c_tau positive, real number, indicating the (fixed) value for c_tau. Ignored if
#' \code{learn_c_xi} is \code{TRUE}, \code{mod_type} is not set to \code{"triple"}  or \code{a_eq_c_tau} is set to \code{TRUE}.
#' The default value is 0.1.
#' @param a_eq_c_xi logical value indicating whether to force \code{a_xi} and \code{c_xi} to be equal.
#' If set to \code{TRUE}, \code{beta_a_xi} and \code{alpha_a_xi} are used as the hyperparameters and \code{beta_c_xi} and \code{alpha_c_xi} are ignored.
#' Ignored if \code{mod_type} is not set to \code{"triple"}. The default value is \code{FALSE}.
#' @param a_eq_c_tau logical value indicating whether to force \code{a_tau} and \code{c_tau} to be equal.
#' If set to \code{TRUE}, \code{beta_a_tau} and \code{alpha_a_tau} are used as the hyperparameters and \code{beta_c_tau} and \code{alpha_c_tau} are ignored.
#' Ignored if \code{mod_type} is not set to \code{"triple"}. The default value is \code{FALSE}.
#' @param learn_kappa2_B logical value indicating whether to learn kappa2_B, the global level of shrinkage for
#' the state variances. The default value is \code{TRUE}.
#' @param learn_lambda2_B logical value indicating whether to learn the lambda2_B parameter,
#' the global level of shrinkage for the mean of the initial values of the states. The default value is \code{TRUE}.
#' @param kappa2_B positive, real number. Initial value of kappa2_B. The default value is 20.
#' @param lambda2_B positive, real number. Initial value of lambda2_B. The default value is
#' @param hyperprior_param \emph{optional} named list containing hyperparameter values.
#' Not all have to be supplied, with those missing being replaced by the default values.
#' Any list elements that are misnamed will be ignored and a warning will be thrown.
#' All hyperparameter values have to be positive, real numbers. The following hyperparameters can be
#' supplied:
#' \itemize{
#' \item \code{e1}: The default value is 0.001.
#' \item \code{e2}: The default value is 0.001.
#' \item \code{d1}: The default value is 0.001.
#' \item \code{d2}: The default value is 0.001.
#' \item \code{beta_a_xi}: The default value is 10.
#' \item \code{beta_a_tau}: The default value is 10.
#' \item \code{alpha_a_xi}: The default value is 5.
#' \item \code{alpha_a_tau}: The default value is 5.
#' }
#' @param sv_param \emph{optional} named list containing hyperparameter values for the stochastic volatility
#' parameters. Not all have to be supplied, with those missing being replaced by the default values.
#' Any list elements that are misnamed will be ignored and a warning will be thrown. Ignored if
#' \code{group} is missing. The following elements can be supplied:
#' \itemize{
#' \item \code{Bsigma_sv}: positive, real number. The default value is 1.
#' \item \code{a0_sv}: positive, real number. The default value is 5.
#' \item \code{b0_sv}: positive, real number. The default value is 1.5.
#' }
#' @param MH_tuning \emph{optional} named list containing values used to tune the MH steps for \code{a_xi} and \code{a_tau}. Not all have to be supplied, with those missing being replaced by the default values.
#' Any list elements that are misnamed will be ignored and a warning will be thrown.
#' The arguments for \code{a_xi}(\code{a_tau}) are only used if \code{learn_a_xi}(\code{learn_a_tau})
#' is set to \code{TRUE}. Arguments ending in "adaptive" are
#' logical values indicating whether or not to make the MH step for the respective parameter adaptive. Arguments ending in "tuning_par" serve two different purposes.
#' If the respective MH step is not set to be adaptive, it acts as the standard deviation of the proposal distribution. If the respective MH step
#' is set to be adaptive, it acts as the initial standard deviation. Arguments ending in "target_rate" define the acceptance rate the algorithm aims to achieve.
#' Arguments ending in "max_adapt" set the maximum value by which the logarithm of the standard deviation of the proposal distribution is adjusted. Finally,
#' arguments ending in "batch_size" set the batch size after which the standard deviation of the proposal distribution is adjusted.
#' The following elements can be supplied:
#' \itemize{
#' \item \code{a_xi_adaptive}: logical value. The default is \code{TRUE}.
#' \item \code{a_xi_tuning_par}: positive, real number. The default value is 1.
#' \item \code{a_xi_target_rate}: positive, real number, between 0 and 1. The default value is 0.44.
#' \item \code{a_xi_max_adapt}: positive, real number. The default value is 0.01.
#' \item \code{a_xi_batch_size}: positive integer. The default value is 50.
#' \item \code{a_tau_adaptive}: logical value. The default is \code{TRUE}.
#' \item \code{a_tau_tuning_par}: positive, real number. The default value is 1.
#' \item \code{a_tau_target_rate}: positive, real number, between 0 and 1. The default value is 0.44.
#' \item \code{a_tau_max_adapt}: positive, real number. The default value is 0.01.
#' \item \code{a_tau_batch_size}: positive integer. The default value is 50.
#' }
#' @param display_progress logical value indicating whether the progress bar and other informative output should be
#' displayed. The default value is \code{TRUE}.
#'  
#' @return  The value returned is a list object of class \code{shrinkDSM} containing
#' \item{\code{beta}}{\code{list} object containing an \code{mcmc.dsm.tvp} object for the parameter draws from the posterior distribution of the centered
#' states, one for each covariate. In the case that there is only one covariate, this becomes just a single \code{mcmc.dsm.tvp} object.}
#' \item{\code{beta_mean}}{\code{mcmc} object containing the parameter draws from the posterior distribution of beta_mean.}
#' \item{\code{theta_sr}}{\code{mcmc} object containing the parameter draws from the posterior distribution of the square root of theta.}
#' \item{\code{tau2}}{\code{mcmc} object containing the parameter draws from the posterior distribution of tau2.}
#' \item{\code{xi2}}{\code{mcmc} object containing the parameter draws from the posterior distribution of xi2.}
#' \item{\code{lambda2}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of lambda2.
#' Not returned if \code{mod_type} is not \code{"triple"}.}
#' \item{\code{kappa2}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of kappa2.
#' Not returned if \code{mod_type} is not \code{"triple"}.}
#' \item{\code{a_xi}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of a_xi.
#' Not returned if \code{learn_a_xi} is \code{FALSE} or \code{mod_type} is \code{"ridge"}.}
#' \item{\code{a_tau}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of a_tau.
#' Not returned if \code{learn_a_tau} is \code{FALSE} or \code{mod_type} is \code{"ridge"}.}
#' \item{\code{c_xi}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of c_xi.
#' Not returned if \code{learn_c_xi} is \code{FALSE} or \code{mod_type} is not \code{"triple"}.}
#' \item{\code{c_tau}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of c_tau.
#' Not returned if \code{learn_c_tau} is \code{FALSE} or \code{mod_type} is not \code{"triple"}.}
#' \item{\code{lambda2_B}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of lambda2_B.
#' Not returned if \code{learn_lambda2_B} is \code{FALSE} or \code{mod_type} is \code{"ridge"}.}
#' \item{\code{kappa2_B}}{\emph{(optional)} \code{mcmc} object containing the parameter draws from the posterior distribution of kappa2_B.
#' Not returned if \code{learn_kappa2_B} is \code{FALSE} or \code{mod_type} is \code{"ridge"}.}
#' \item{\code{MH_diag}}{\emph{(optional)} named list containing statistics for assessing MH performance. Not returned if no MH steps are required
#' or none of them are specified to be adaptive.}
#' \item{\code{priorvals}}{\code{list} object containing hyperparameter values of the prior distributions, as specified by the user.}
#' \item{\code{model}}{\code{list} object containing the model matrix, model response and formula used.}
#' \item{\code{summaries}}{\code{list} object containing a collection of summary statistics of the posterior draws.}
#' 
#' To display the output, use \code{plot} and \code{summary}. The \code{summary} method displays the specified prior values stored in
#' \code{priorvals} and the posterior summaries stored in \code{summaries}, while the \code{plot} method calls \code{coda}'s \code{plot.mcmc}
#' or the \code{plot.mcmc.dsm.tvp} method. Furthermore, all functions that can be applied to \code{coda::mcmc} objects
#' (e.g. \code{coda::acfplot}) can be applied to all output elements that are \code{coda} compatible.
#' 
#' @examples
#' \donttest{
#' set.seed(123)
#' data("gastric")
#'
#' # Create intervals for piecewise exponential model
#' intervals <- divisionpoints(gastric$time, gastric$status, 2)
#'
#' # Estimate baseline model
#' mod <- shrinkDSM(time ~ radiation, gastric,
#'                  delta = gastric$status, S = intervals)
#'
#' # Estimate model with different prior setup
#' mod2 <- shrinkDSM(time ~ radiation, gastric,
#'                  delta = gastric$status, S = intervals,
#'                  mod_type = "triple")
#'
#' # Change some of the hyperparameters
#' mod3 <- shrinkDSM(time ~ radiation, gastric,
#'                  delta = gastric$status, S = intervals,
#'                  mod_type = "triple",
#'                  hyperprior_param = list(beta_a_xi = 5,
#'                                          alpha_a_xi = 10))
#' }
#' @export
shrinkDSM <- function(formula,
                      data,
                      mod_type = "double",
                      delta,
                      S,
                      group,
                      subset,
                      niter = 10000,
                      nburn = round(niter/2),
                      nthin = 1,
                      learn_a_xi = TRUE,
                      learn_a_tau = TRUE,
                      a_xi = 0.1,
                      a_tau = 0.1,
                      learn_c_xi = TRUE,
                      learn_c_tau = TRUE,
                      c_xi = 0.1,
                      c_tau = 0.1,
                      a_eq_c_xi = FALSE,
                      a_eq_c_tau = FALSE,
                      learn_kappa2_B = TRUE,
                      learn_lambda2_B = TRUE,
                      kappa2_B = 20,
                      lambda2_B = 20,
                      hyperprior_param,
                      sv_param,
                      MH_tuning,
                      display_progress = TRUE){
  assert(missing(group) || length(group) == nrow(data), "Grouping indicator, group, has to be omitted or the same length as data")
  assert(length(delta) == nrow(data), "Status indicator, delta, has to be the same length as data")

  ## Implementation of formula interface:
  mf <- match.call(expand.dots = FALSE)
  m <- match(x = c("formula", "data", "subset"),
             table = names(mf), nomatch = 0L)
  mf <- mf[c(1L, m)]

  mf$drop.unused.levels <- TRUE
  mf[[1L]] <- quote(stats::model.frame)
  mf <- eval(expr = mf, envir = parent.frame())
  for(name in names(mf)){
    assert(!any(class(mf[[name]]) %in% c("POSIXct", "POSIXt","Date")), "No date variables allowed as predictors")
  }
  # Create Vector y
  y <- model.response(mf, "numeric")
  mt <- attr(x = mf, which = "terms")
  # Create Matrix X with dummies and transformations
  z <- model.matrix(object = mt, data = mf)

  colnames(z)[colnames(z) == "(Intercept)"] <- "Intercept"

  # Check that there are no NAs in y and x
  assert(!any(is.na(y)), "No NA values are allowed in survival time")
  assert(all(y>0), "Survival times must be positive. Zeros are not allowed.")
  assert(!any(is.na(z)), "No NA values are allowed in covariates")

  assert(mod_type %in% c("double", "triple", "ridge"), "Allowed model types are: ridge, double, and triple")
  # default hyperparameter values
  default_hyper <- list(c0 = 2.5,
                        g0 = 5,
                        G0 = 5 / (2.5 - 1),
                        e1 = 0.001,
                        e2 = 0.001,
                        d1 = 0.001,
                        d2 = 0.001,
                        beta_a_xi = 10,
                        beta_a_tau = 10,
                        alpha_a_xi = 5,
                        alpha_a_tau = 5,
                        beta_c_xi = 2,
                        beta_c_tau = 2,
                        alpha_c_xi = 5,
                        alpha_c_tau = 5)

  if (!missing(group)) {
    group <- checkvalues(group)
    default_hyper$sigma2_phi <- rep(1, length(unique(group$values)))
  } else {
    group_sort <- c(0)
    default_hyper$sigma2_phi <- c(0)
  }

  # default sv params
  default_hyper_sv <- list(Bsigma_sv = 1,
                           a0_sv = 5,
                           b0_sv = 1.5)

  # default tuning parameters
  default_tuning_par <- list(a_xi_adaptive = TRUE,
                             a_xi_tuning_par = 1,
                             a_xi_target_rate = 0.44,
                             a_xi_max_adapt = 0.01,
                             a_xi_batch_size = 50,
                             a_tau_adaptive = TRUE,
                             a_tau_tuning_par = 1,
                             a_tau_target_rate = 0.44,
                             a_tau_max_adapt = 0.01,
                             a_tau_batch_size = 50,
                             c_xi_adaptive = TRUE,
                             c_xi_tuning_par = 1,
                             c_xi_target_rate = 0.44,
                             c_xi_max_adapt = 0.01,
                             c_xi_batch_size = 50,
                             c_tau_adaptive = TRUE,
                             c_tau_tuning_par = 1,
                             c_tau_target_rate = 0.44,
                             c_tau_max_adapt = 0.01,
                             c_tau_batch_size = 50)

  ## Merge user supplied and default parameters
  if (missing(hyperprior_param)) {
    hyperprior_param <- default_hyper
  } else {
    hyperprior_param <- list_merger(default_hyper, hyperprior_param)
  }

  if (missing(sv_param)) {
    sv_param <- default_hyper_sv
  } else {
    sv_param <- list_merger(default_hyper_sv, sv_param)
  }

  if (missing(MH_tuning)) {
    MH_tuning <- default_tuning_par
  } else {
    MH_tuning <- list_merger(default_tuning_par, MH_tuning)
  }
  # Check if all numeric inputs are correct
  to_test_num <- list(lambda2_B = lambda2_B,
                      kappa2_B = kappa2_B,
                      a_xi = a_xi,
                      a_tau = a_tau,
                      c_xi = c_xi,
                      c_tau = c_tau)
  if(!missing(group)) {
    cond = length(hyperprior_param$sigma2_phi) == length(unique(group$values)) && all(!sapply(hyperprior_param$sigma2_phi, numeric_input_bad))
    assert(cond, "all elements of sigma2_phi have to be positive real numbers")}
  if (missing(hyperprior_param) == FALSE){
    to_test_num <- c(to_test_num, hyperprior_param[names(hyperprior_param) != "sigma2_phi"])
  }

  if (missing(sv_param) == FALSE){
    to_test_num <- c(to_test_num, sv_param)
  }

  if (missing(MH_tuning) == FALSE){
    to_test_num <- c(to_test_num, MH_tuning[!grepl("(batch|adaptive)", names(MH_tuning))])
  }

  bad_inp <- sapply(to_test_num, numeric_input_bad)


  if (any(bad_inp)){
    stand_names <- c(names(default_hyper), names(default_hyper_sv), "lambda2_B", "kappa2_B", "a_xi", "a_tau", "c_xi", "c_tau")
    bad_inp_names <- names(to_test_num)[bad_inp]
    bad_inp_names <- bad_inp_names[bad_inp_names %in% stand_names]
    stop(paste0(paste(bad_inp_names, collapse = ", "),
                ifelse(length(bad_inp_names) == 1, " has", " have"),
                " to be a real, positive number"))
  }


  # Check the adapt_rates seperately
  if (any(0 > MH_tuning[grepl("rate", names(MH_tuning))] | MH_tuning[grepl("rate", names(MH_tuning))] > 1)) {
    stop("all target_rate parameters in MH_tuning have to be > 0 and < 1")
  }

  # Check if all integer inputs are correct
  to_test_int <- c(niter = niter,
                   nburn = nburn,
                   nthin = nthin,
                   MH_tuning[grepl("batch", names(MH_tuning))])

  bad_int_inp <- sapply(to_test_int, int_input_bad)

  if (any(bad_int_inp)){
    bad_inp_names <- names(to_test_int)[bad_int_inp]
    stop(paste0(paste(bad_inp_names, collapse = ", "),
                ifelse(length(bad_inp_names) == 1, " has", " have"),
                " to be a single, positive integer"))

  }

  if ((niter - nburn) < 2){
    stop("niter has to be larger than or equal to nburn + 2")
  }

  if (nthin == 0){
    stop("nthin can not be 0")
  }

  if ((niter - nburn)/2 < nthin){
    stop("nthin can not be larger than (niter - nburn)/2")
  }

  # Check if all boolean inputs are correct
  to_test_bool <- c(learn_lambda2_B = learn_lambda2_B,
                    learn_kappa2_B = learn_kappa2_B,
                    learn_a_xi = learn_a_xi,
                    learn_a_tau = learn_a_tau,
                    display_progress = display_progress,
                    MH_tuning[grepl("adaptive", names(MH_tuning))])

  bad_bool_inp <- sapply(to_test_bool, bool_input_bad)

  if (any(bad_bool_inp)){
    bad_inp_names <- names(to_test_bool)[bad_bool_inp]
    stop(paste0(paste(bad_inp_names, collapse = ", "),
                ifelse(length(bad_inp_names) == 1, " has", " have"),
                " to be a single logical value"))

  }

  # Check if formula is a formula
  if (inherits(formula, "formula") == FALSE){
    stop("formula is not of class formula")
  }



  # Extract colnames

  d <- ncol(z)
  if (!is.null(colnames(z))){
    col_names <- colnames(z)
  } else {
    col_names <- as.character(1:d)
  }

  # Sort objects
  order <- order(y, decreasing = TRUE)
  y_sort <- y[order]
  z_sort <- z[order, ]
  if(!missing(group)){
    # Group is guaranteed to be the values here
    group_sort <- group$values[order]
  }

  if(2 %in% delta){delta <- delta - 1}
  assert(all(delta %in% c(0,1)),
    "delta must contain only 0/1, 1/2, or TRUE/FALSE")

  z_sort <- as.matrix(z_sort)
  delta_sort <- as.matrix(as.integer(delta[order]))


  runtime <- system.time({
  res <- do_shrinkDSM(y_sort,
                      z_sort,
                      mod_type,
                      delta_sort,
                      S,
                      group_sort,
                      niter,
                      nburn,
                      nthin,
                      hyperprior_param$d1,
                      hyperprior_param$d2,
                      hyperprior_param$e1,
                      hyperprior_param$e2,
                      hyperprior_param$sigma2_phi,
                      learn_lambda2_B,
                      learn_kappa2_B,
                      lambda2_B,
                      kappa2_B,
                      learn_a_xi,
                      learn_a_tau,
                      a_xi,
                      a_tau,
                      learn_c_xi,
                      learn_c_tau,
                      c_xi,
                      c_tau,
                      a_eq_c_xi,
                      a_eq_c_tau,
                      MH_tuning$a_xi_tuning_par,
                      MH_tuning$a_tau_tuning_par,
                      MH_tuning$c_xi_tuning_par,
                      MH_tuning$c_tau_tuning_par,
                      hyperprior_param$beta_a_xi,
                      hyperprior_param$beta_a_tau,
                      hyperprior_param$alpha_a_xi,
                      hyperprior_param$alpha_a_tau,
                      hyperprior_param$beta_c_xi,
                      hyperprior_param$beta_c_tau,
                      hyperprior_param$alpha_c_xi,
                      hyperprior_param$alpha_c_tau,
                      sv_param$Bsigma_sv,
                      sv_param$a0_sv,
                      sv_param$b0_sv,
                      display_progress,
                      unlist(MH_tuning[grep("adaptive", names(MH_tuning))]),
                      unlist(MH_tuning[grep("target", names(MH_tuning))]),
                      unlist(MH_tuning[grep("max", names(MH_tuning))]),
                      unlist(MH_tuning[grep("size", names(MH_tuning))]))
  })
  if(display_progress == TRUE){
    cat("Timing (elapsed): ", file=stderr())
    cat(runtime["elapsed"], file=stderr())
    cat(" seconds.\n", file=stderr())
    cat(round((niter + nburn)/runtime[3]), "iterations per second.\n\n", file=stderr())
    cat("Converting results to coda objects and summarizing draws... ", file=stderr())
  }

  # Remove empty storage elements
  res[sapply(res, function(x) 0 %in% dim(x))] <- NULL
  res$MH_diag[sapply(res$MH_diag, function(x) 0 %in% dim(x))] <- NULL

    # Create object to hold prior values
  res$priorvals <- c(hyperprior_param,
                     sv_param,
                     a_xi = a_xi,
                     a_tau = a_tau,
                     lambda2_B = lambda2_B,
                     kappa2_B = kappa2_B)

  # Add data to output
  res[["model"]] <- list()
  res$model$z <- z
  res$model$y <- y
  res$model$formula <- formula
  res$model$xlevels <- .getXlevels(mt, mf)
  res$model$terms <- mt
  res$model$model <- mf


  # add attributes to the individual objects if they are distributions or individual statistics
  nsave <- floor((niter - nburn)/nthin)
  for (i in names(res)){

    attr(res[[i]], "type") <- ifelse(nsave %in% dim(res[[i]]), "sample", "stat")

    # Name each individual sample for plotting frontend
    if (attr(res[[i]], "type") == "sample"){

      if (i == "phi"){

        colnames(res[[i]]) <- paste0(i, unique(group_sort))

      } else if (dim(res[[i]])[2] == d){

        colnames(res[[i]]) <- paste0(i, "_",  col_names)

      } else if (dim(res[[i]])[2] == 2 * d) {

        colnames(res[[i]]) <- paste0(i, "_", rep(col_names, 2))


      } else {
        colnames(res[[i]]) <- i

      }
    }

    # Change objects to be coda compatible
    # Only apply to posterior samples
    if (attr(res[[i]], "type") == "sample"){

      # Differentiate between TVP and non TVP
      if (is.na(dim(res[[i]])[3]) == FALSE){

        # Create a sub list containing an mcmc object for each parameter in TVP case
        dat <- res[[i]]
        res[[i]] <- list()
        for (j in 1:dim(dat)[2]){
          res[[i]][[j]] <- as.mcmc(t(dat[, j, ]), start = niter - nburn, end = niter, thin = nthin)
          colnames(res[[i]][[j]]) <- paste0(i, "_", j, "_", 1:ncol(res[[i]][[j]]))

          # make it of class mcmc.tvp for custom plotting function
          class(res[[i]][[j]]) <- c("mcmc.dsm.tvp", "mcmc")
          attr(res[[i]][[j]], "S") <- S#c(S, max(res$model$y))
          attr(res[[i]][[j]], "lastsurvtime") <- max(res$model$y)
          attr(res[[i]][[j]], "type") <- "sample"
        }

        if (length(res[[i]]) == 1){
          res[[i]] <- res[[i]][[j]]
        }

        # Make it of type 'sample' again
        attr(res[[i]], "type") <- "sample"

        # Rename
        if (dim(dat)[2] > 1){
          names(res[[i]]) <- colnames(dat)
        }


      } else {

        res[[i]] <- as.mcmc(res[[i]], start = niter - nburn, end = niter, thin = nthin)

      }
    }

    # Create summary of posterior
    if (is.list(res[[i]]) == FALSE & attr(res[[i]], "type") == "sample") {
      if (i != "theta_sr" & i != "beta") {
        res$summaries[[i]] <- t(apply(res[[i]], 2, function(x){

          obj <- as.mcmc(x, start = niter - nburn, end = niter, thin = nthin)
          ESS <- tryCatch(coda::effectiveSize(obj),
                          error = function(err) {
                            warning("Calculation of effective sample size failed for one or more variable(s). This can happen if the prior placed on the model induces extreme shrinkage.")
                            return(NA)
                          }, silent = TRUE)

          return(c("mean" = mean(obj),
                   "sd" = sd(obj),
                   "median" = median(obj),
                   "HPD" = HPDinterval(obj)[c(1, 2)],
                   "ESS" = round(ESS)))
        }))
      } else if (i == "theta_sr") {
        res$summaries[[i]] <- t(apply(res[[i]], 2, function(x){

          obj <- as.mcmc(abs(x), start = niter - nburn, end = niter, thin = nthin)
          ESS <- tryCatch(coda::effectiveSize(obj),
                          error = function(err) {
                            warning("Calculation of effective sample size failed for one or more variable(s). This can happen if the prior placed on the model induces extreme shrinkage.")
                            return(NA)
                          }, silent = TRUE)

          return(c("mean" = mean(obj),
                   "sd" = sd(obj),
                   "median" = median(obj),
                   "HPD" = HPDinterval(obj)[c(1, 2)],
                   "ESS" = round(ESS)))
        }))
      }
    }
  }

  if (display_progress == TRUE) {
    cat("Done!\n", file = stderr())
  }

  # add some attributes for the methods and plotting
  attr(res, "class") <- "shrinkDSM"
  attr(res, "S") <- S#c(S, max(res$model$y))
  attr(res, "lastsurvtime") <- max(res$model$y)
  if (!missing(group)) {
    attr(res, "group") <- group
  }
  attr(res, "learn_a_xi") <- learn_a_xi
  attr(res, "learn_a_tau") <- learn_a_tau
  attr(res, "learn_kappa2_B") <- learn_kappa2_B
  attr(res, "learn_lambda2_B") <- learn_lambda2_B
  attr(res, "niter") <- niter
  attr(res, "nburn") <- nburn
  attr(res, "nthin") <- nthin
  attr(res, "colnames") <-  col_names

  return(res)
}

Try the shrinkDSM package in your browser

Any scripts or data that you put into this service are public.

shrinkDSM documentation built on Sept. 6, 2021, 9:07 a.m.