Fit Bayesian Linear and Generalized Linear Mixed-Effects Models
Maximum a posteriori estimation for linear and generalized
linear mixed-effects models in a Bayesian setting. Built off of
lmer.
blmer(formula, data = NULL, REML = TRUE,
control = lmerControl(), start = NULL, verbose = 0L,
subset, weights, na.action, offset, contrasts = NULL,
devFunOnly = FALSE, cov.prior = wishart,
fixef.prior = NULL, resid.prior = NULL, ...)
bglmer(formula, data = NULL, family = gaussian,
control = glmerControl(), start = NULL, verbose = 0L,
nAGQ = 1L, subset, weights, na.action, offset,
contrasts = NULL, mustart, etastart,
devFunOnly = FALSE, cov.prior = wishart,
fixef.prior = NULL, ...)cov.prior |
a BLME prior or list of priors with allowable
distributions: |
fixef.prior |
a BLME prior of family |
resid.prior |
a BLME prior of family |
start |
like the |
formula, data, REML, family, control, verbose, nAGQ,
mustart, etastart, devFunOnly, ... |
model specification arguments as in |
subset, weights, na.action, offset, contrasts |
further model
specification arguments as in |
The bulk of the usage for blmer and bglmer closely
follows the functions lmer and
glmer. Those help pages provide a good overview of
fitting linear and generalized linear mixed models. The primary
distinction is that blmer and bglmer allow the user to
do Bayesian inference or penalized maximum likelihood, with priors imposed on the different
model components. For the specifics of any distribution listed below,
see the distributions page.
Covariance Prior
The cov.prior argument applies a prior over the
covariance matrix of the random effects/modeled coefficients.
As there is one covariance matrix for every named grouping factor -
that is every element that appears to the right of a vertical bar
("|") in the model formula - it is possible to apply as many
different priors as there are said factors.
The general formats of an argument to blmer or bglmer
for such a prior are of the form:
cov.prior = factor.name ~ covariance.distribution(option1 = value1, ...)
cov.prior = list(fc.nm ~ dist1, fc.nm ~ dist2, ..., default.distribution)
If the “factor.name ~” construct is ommitted, the prior
is interpretted as a default and applied to all factors that
lack specific priors of their own. Options are not required,
but permit fine-tuning of the model.
Supported distributions are gamma, invgamma, wishart,
invwishart, NULL, and custom.
The common.scale option, a logical, determines whether or
not the prior applies to in the absolute-real world
sense (value = FALSE), or if the prior is applied to the random effect
covariance divided by the estimated residual variance (TRUE). As a practical matter,
when false computation can be slower as the profiled common scale may
no longer have a closed-form solution. As such, the default for all
cases is TRUE.
Other options are specified along with the specific distributions and defaults are explained in the blme distributions page.
Fixed Effects Prior
Priors on the fixed effects, or unmodeled coefficients, are specified in a fashion similar to that of covariance priors. The general format is
fixef.prior = multivariate.distribution(options1 = value1, ...)
At present, the implemented multivariate distributions are normal, t,
horseshoe, and NULL. t and horseshoe priors cannot be used
when REML is TRUE, as that integral does not have a closed form solution.
Residual Variance Prior
The general format for a residual variance prior is the same as for a fixed
effect prior. The supported distributions are point, gamma,
invgamma.
An object of class "bmerMod", for which many methods
are available. See there for details.
data("sleepstudy", package = "lme4")
### Examples using a covariance prior ##
# Here we are ignoring convergence warnings just to illustate how the package
# is used: this is not a good idea in practice..
control <- lmerControl(check.conv.grad = "ignore")
(fm1 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
control = control,
cov.prior = gamma))
(fm2 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
control = control,
cov.prior = gamma(shape = 2, rate = 0.5, posterior.scale = 'sd')))
(fm3 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
control = control,
cov.prior = wishart))
(fm4 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
control = control,
cov.prior = invwishart(df = 5, scale = diag(0.5, 2))))
# Custom prior
penaltyFn <- function(sigma)
dcauchy(sigma, 0, 10, log = TRUE)
(fm5 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
cov.prior = custom(penaltyFn, chol = TRUE, scale = "log")))
### Examples using a fixed effect prior ###
(fm6 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
cov.prior = NULL,
fixef.prior = normal))
(fm7 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
cov.prior = NULL,
fixef.prior = normal(cov = diag(0.5, 2), common.scale = FALSE)))
### Example using a residual variance prior ###
# This is the "eight schools" data set; the mode should be at the boundary
# of the space.
control <- lmerControl(check.conv.singular = "ignore",
check.nobs.vs.nRE = "ignore",
check.nobs.vs.nlev = "ignore")
y <- c(28, 8, -3, 7, -1, 1, 18, 12)
sigma <- c(15, 10, 16, 11, 9, 11, 10, 18)
g <- 1:8
(schools <- blmer(y ~ 1 + (1 | g), control = control, REML = FALSE,
resid.prior = point, cov.prior = NULL,
weights = 1 / sigma^2))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.