Nothing
#' 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'.
#' @returns 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 <- numeric(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()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.