R/mc_vreg.R

In mcmcsae: Markov Chain Monte Carlo Small Area Estimation

#' Create a model component object for a regression component in the variance function of a
#' gaussian sampling distribution
#'
#' This function is intended to be used on the right hand side of the \code{formula.V} argument to
#' \code{\link{create_sampler}} or \code{\link{generate_data}}.
#'
#' @export
#' @param formula a formula for the regression effects explaining the log-variance.
#'  Variable names are looked up in the data frame passed as \code{data} argument to
#'  \code{\link{create_sampler}} or \code{\link{generate_data}}, or in \code{environment(formula)}.
#' @param remove.redundant whether redundant columns should be removed from the design matrix.
#'  Default is \code{FALSE}.
#' @param sparse whether the model matrix associated with \code{formula} should be sparse.
#'  The default is determined by a simple heuristic based on storage size.
#' @param X a (possibly sparse) design matrix can be specified directly, as an alternative to the
#'  creation of one based on \code{formula}. If \code{X} is specified \code{formula} is ignored.
#' @param prior prior specification for the coefficients. Currently only
#'  normal priors are supported, specified using function \code{\link{pr_normal}}.
#' @param Q0 prior precision matrix for the regression effects. The default is a
#'  zero matrix corresponding to a noninformative improper prior.
#'  DEPRECATED, please use argument \code{prior} instead, i.e.
#'  \code{prior = pr_normal(mean = b0.value, precision = Q0.value)}.
#' @param b0 prior mean for the regression effect. Defaults to a zero vector.
#'  DEPRECATED, please use argument \code{prior} instead, i.e.
#'  \code{prior = pr_normal(mean = b0.value, precision = Q0.value)}.
#' @param name the name of the model component. This name is used in the output of the
#'  MCMC simulation function \code{\link{MCMCsim}}. By default the name will be 'vreg'
#'  with the number of the variance model term attached.
#' @return An object with precomputed quantities and functions for sampling from
#'  prior or conditional posterior distributions for this model component. Intended
#'  for internal use by other package functions.
#' @references
#'  E. Cepeda and D. Gamerman (2000).
#'    Bayesian modeling of variance heterogeneity in normal regression models.
#'    Brazilian Journal of Probability and Statistics, 207-221.
#'
#'  T.I. Lin and W.L. Wang (2011).
#'    Bayesian inference in joint modelling of location and scale parameters
#'    of the t distribution for longitudinal data.
#'    Journal of Statistical Planning and Inference 141(4), 1543-1553.
vreg <- function(formula=NULL, remove.redundant=FALSE, sparse=NULL, X=NULL,
                 prior=NULL, Q0=NULL, b0=NULL, name="") {

  e <- sys.frame(-2L)
  type <- "vreg"
  if (name == "") stop("missing model component name")

  if (e$Q0.type == "symm") stop("TBI: vreg component with (compatible) non-diagonal sampling variance matrix")

  if (is.null(X)) {
    X <- model_matrix(formula, e$data, sparse=sparse)
    if (remove.redundant) X <- remove_redundancy(X)
  } else
    X <- economizeMatrix(X, strip.names=FALSE, check=TRUE)
  if (nrow(X) != e$n) stop("design matrix with incompatible number of rows")
  e$coef.names[[name]] <- colnames(X)
  X <- unname(X)
  q <- ncol(X)

  if (!is.null(b0) || !is.null(Q0)) {
    warn("arguments 'b0' and 'Q0' are deprecated; please use argument 'prior' instead")
    if (is.null(prior)) prior <- pr_normal(mean = if (is.null(b0)) 0 else b0, precision = if (is.null(Q0)) 0 else Q0)
  } else {
    if (is.null(prior)) prior <- pr_normal(mean=0, precision=0)
  }
  if (prior$type != "normal") stop("only a normal prior is currently supported for 'vreg' model component effects")
  prior$init(q, e$coef.names[[name]])
  informative.prior <- prior$informative
  Q0 <- prior$precision
  zero.mean <- !informative.prior || all(prior$mean == 0)
  if (zero.mean) {
    prior$mean <- 0
    Q0b0 <- 0
  } else {
    if (length(prior$mean) == 1L)
      Q0b0 <- Q0 %m*v% rep.int(prior$mean, q)
    else
      Q0b0 <- Q0 %m*v% prior$mean
  }
  rprior <- function(p) prior$rprior(p)

  sumX <- 0.5 * colSums(X) - Q0b0

  MVNsampler <- create_TMVN_sampler(Q=0.5*crossprod(X) + Q0, name=name)

  compute_Qfactor <- function(p) exp(X %m*v% (- p[[name]]))

  # function that creates an oos prediction function (closure) based on new data
  make_predict_Vfactor <- function(newdata) {
    Xnew <- model_matrix(formula, newdata)
    # TODO remove columns not corresponding to columns in X
    Xnew <- unname(Xnew)
    rm(newdata)
    function(p) exp(Xnew %m*v% p[[name]])
  }

  if (e$prior.only) return(environment())

  sigma.fixed <- e$sigma.fixed

  proposal_scale <- 2.4/sqrt(q)
  draw <- function(p) {
    # TODO consider case e$Q0.type="symm"
    # random walk MH proposal; see Lin and Wang (2011)
    gamma <- p[[name]]
    gamma.star <- gamma + MVNsampler$draw(p, proposal_scale)[[name]]
    if (sigma.fixed)
      Qres <- p[["Q_"]] * p[["e_"]]^2
    else
      Qres <- p[["Q_"]] * (1 / p[["sigma_"]]^2) * p[["e_"]]^2
    Qratio <- exp(X %m*v% (gamma - gamma.star))
    log.ar <- 0.5 * (dotprodC(gamma, Q0 %m*v% gamma) - dotprodC(gamma.star, Q0 %m*v% gamma.star)) +
              dotprodC(gamma - gamma.star, sumX) + 0.5 * sum((1 - Qratio) * Qres)
    if (log(runif(1L)) < log.ar) {
      p[[name]] <- gamma.star
      p$Q_ <- p[["Q_"]] * Qratio
    }
    p
  }

  start <- function(p) {
    if (is.null(p[[name]])) p[[name]] <- runif(q, -1e-3, 1e-3)
    if (length(p[[name]]) != q) stop("wrong length for start value '", name, "'")
    p
  }

  rm(e)
  environment()
}