R/mc_gl.R

In mcmcsae: Markov Chain Monte Carlo Small Area Estimation

#' Create a model object for group-level regression effects within a
#' generic random effects component.
#'
#' This function is intended to be used to specify the \code{formula.gl} argument to
#' the \code{\link{gen}} model component specification function.
#' Group-level predictors and hierarchical centering are
#' not used by default, and they currently cannot be used in a model component that is sampled
#' together with another model component in the same Gibbs block.
#'
#' @export
#' @param formula a formula specifying the group-level predictors to be used within a model
#'  component. If no \code{data} is supplied the group-level predictors are derived as
#'  group-level means from the unit-level data passed as \code{data} argument to
#'  \code{\link{create_sampler}} or \code{\link{generate_data}}.
#' @param remove.redundant whether redundant columns should be removed from the design matrix.
#'  Default is \code{FALSE}.
#' @param prior prior specification for the group-level effects. Currently only
#'  normal priors with mean 0 can be specified, using function \code{\link{pr_normal}}.
#' @param Q0 prior precision matrix for the group-level 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(precision = Q0.value)}.
#' @param data group-level data frame in which the group-level variables specified in
#'  \code{formula} are looked up.
#' @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 this name will be
#'  the name of the corresponding generic random effects component appended by '_gl'.
#' @return An object with precomputed quantities for sampling from
#'  prior or conditional posterior distributions for this model component. Only intended
#'  for internal use by other package functions.
glreg <- function(formula=NULL, remove.redundant=FALSE, prior=NULL, Q0=NULL,
                  data=NULL, name="") {

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

  if (is.null(data)) {
    # derive group-level design matrix as group means of unit-level data
    if (is.null(formula)) stop("no group-level formula: cannot derive group-level design matrix")
    if (intercept_only(formula)) {
      X <- matrix(rep.int(1, e$l), ncol=1L, dimnames=list(NULL, "(Intercept)"))
    } else {
      X <- model_matrix(formula, e$e$data, sparse=FALSE)
      if (remove.redundant) X <- remove_redundancy(X)
      if (is.null(e$factor)) stop("cannot derive group-level design matrix")
      factor.info <- get_factor_info(e$factor, e$e$data)
      if (any("spline" == factor.info$types)) stop("unsupported combination: splines and group-level covariates")
      fac <- combine_factors(factor.info$variables, e$e$data)
      X <- economizeMatrix(crossprod(aggrMatrix(fac, mean=TRUE), X), strip.names=FALSE, check=TRUE)
      rm(fac)
    }
  } else {  # formula + group-level data provided
    X <- model_matrix(formula, data, sparse=FALSE)
    if (remove.redundant) X <- remove_redundancy(X)
    # TODO if formula is used, match levels of factor to glp$data
  }
  if (nrow(X) != e$l) stop("wrong number of rows of group-level design matrix")
  if (!is.null(colnames(X))) {
    e$e$coef.names[[name]] <- colnames(X)
  }
  X <- unname(X)  # l x p0

  p0 <- ncol(X)
  # if p0 > 1, each varying effect has its own coefficient and the same design matrix glp$X
  q <- p0 * e$q0  # number of gl coefficients

  if (!is.null(Q0)) {
    warn("argument 'Q0' is deprecated; please use argument 'prior' instead")
    if (is.null(prior)) prior <- pr_normal(precision = if (is.null(Q0)) 0 else Q0)
  } else {
    if (is.null(prior)) prior <- pr_normal(mean=0, precision=0)
  }

  if (prior$type != "normal" || any(prior$mean != 0)) stop("only a normal prior with mean zero is currently supported for group-level effects")
  # if modeled.Q need sparse Q0 and XX in order to combine them easily using a block-diagonal sparse template
  prior$init(q, e$e$coef.names[[name]], sparse=e$e$modeled.Q, sigma=!e$e$sigma.fixed)
  informative.prior <- prior$informative
  Q0 <- prior$precision
  Q0b0 <- rep.int(0, q)  # need this in draw function, even though only b0=0 is supported

  # for modeled Q, the XX block of XX.ext is updated -> choose sparse
  XX.ext <- economizeMatrix(bdiag(e$XX, Q0), symmetric=TRUE, sparse=if (e$e$modeled.Q) TRUE else NULL)
  IU0 <- economizeMatrix(cbind(Diagonal(e$l), -X), drop.zeros=TRUE)
  if (e$Leroux_update) {
    # TODO crossprod_sym may introduce 0 fill-in --> kronecker template may fail for "unstructured" or "diagonal" var
    #      may drop (prune) 0s, but need to do that in each next crossprod_sym call too! (unavoidable?)
    w <- runif(1L, 0.25, 0.75)
    QA.ext <- crossprod_sym(IU0, e$mat_sum_Leroux(e$QA, e$idL, w, 1-w))
    rm(w)
  } else {
    QA.ext <- economizeMatrix(crossprod_sym(IU0, e$QA), symmetric=TRUE, drop.zeros=TRUE)
  }
  if (!is.null(e$R)) {
    R <- economizeMatrix(crossprod(kronecker(IU0, Diagonal(e$q0)), e$R), allow.tabMatrix=FALSE)
  }

  rm(data, e)
  environment()
}