R/misc.R

In csetraynor/mctte: Compartmental models applied to time-to-event.

## Includes code snippets of the rstanarm package 
# Copyright (C) 2013, 2014, 2015, 2016, 2017 Trustees of Columbia University

# Check and set scale parameters for priors
#
# @param scale Value of scale parameter (can be NULL).
# @param default Default value to use if \code{scale} is NULL.
# @param link String naming the link function or NULL.
# @return If a probit link is being used, \code{scale} (or \code{default} if
#   \code{scale} is NULL) is scaled by \code{dnorm(0) / dlogis(0)}. Otherwise
#   either \code{scale} or \code{default} is returned.
set_prior_scale <- function(scale, default, link) {
  stopifnot(is.numeric(default), is.character(link) || is.null(link))
  if (is.null(scale)) 
    scale <- default
  if (isTRUE(link == "probit"))
    scale <- scale * dnorm(0) / dlogis(0)
  
  return(scale)
}

# Test for a given family
#
# @param x A character vector (probably x = family(fit)$family)
is.binomial <- function(x) x == "binomial"
is.gaussian <- function(x) x == "gaussian"
is.gamma    <- function(x) x == "Gamma"
is.ig       <- function(x) x == "inverse.gaussian"
is.nb       <- function(x) x == "neg_binomial_2"
is.poisson  <- function(x) x == "poisson"
is.beta     <- function(x) x == "beta" || x == "Beta regression"


# Maybe broadcast 
#
# @param x A vector or scalar.
# @param n Number of replications to possibly make. 
# @return If \code{x} has no length the \code{0} replicated \code{n} times is
#   returned. If \code{x} has length 1, the \code{x} replicated \code{n} times
#   is returned. Otherwise \code{x} itself is returned.
maybe_broadcast <- function(x, n) {
  if (!length(x)) {
    rep(0, times = n)
  } else if (length(x) == 1L) {
    rep(x, times = n)
  } else {
    x
  }
}

# Return the default scale parameter for 'prior_aux'.
#
# @param basehaz A list with info about the baseline hazard; see 'handle_basehaz'.
# @return A scalar.
get_default_aux_scale2 <- function(basehaz) {
  nm <- basehaz
  if (nm %in% c("weibull", "weibull-aft", "gompertz")) 2 else 20
}

#' Calculate dt
#' Obtain time points for numerical differentiation in Stan. 
#' The default method is to use the symmetric difference qution.
#' @param time time variable
#' @param h the small number default to 1e-3
#' @param ... currently unused
#' @export
#' @description For internal usage
calculate_dt <- function(time, h = 1e-3){
  N <-  length(time)*2
  time_dt <- rep(NA, N)
  odd_index <- odd(1:N) 
  even_index <- even(1:N)
  
  time_dt[odd_index] <- time - h
  time_dt[even_index] <- time + h
  
  as.numeric(time_dt)
}

### helpers
odd <- function(x) x%%2 != 0 
even <- function(x) x%%2 == 0 

'%!in%' <- function(x,y)!('%in%'(x,y))

nlist <- function(...) {
  m <- match.call()
  out <- list(...)
  no_names <- is.null(names(out))
  has_name <- if (no_names) FALSE else nzchar(names(out))
  if (all(has_name))
    return(out)
  nms <- as.character(m)[-1L]
  if (no_names) {
    names(out) <- nms
  } else {
    names(out)[!has_name] <- nms[!has_name]
  }
  
  return(out)
}


#### helpers 
get_knots <- function(num_knots, tt, spline_degree){
  knots <- unname(quantile(c(0, max(tt)),probs=seq(from=0, to=1,length.out = num_knots)))
  knots
}


am <- function(x, ...) as.matrix(x, ...)
ai <- function(x, ...) as.integer(x, ...)