View source: R/vglmer_regression.R
vglmer | R Documentation |
This function estimates hierarchical models using mean-field variational
inference. vglmer
accepts standard syntax used for lme4
, e.g.,
y ~ x + (x | g)
. Options are described below. Goplerud (2022a; 2022b)
provides details on the variational algorithms.
vglmer(formula, data, family, control = vglmer_control())
formula |
|
data |
|
family |
Options are "binomial", "linear", or "negbin" (experimental).
If "binomial", outcome must be either binary ( |
control |
Adjust internal options for estimation. Must use an object created by vglmer_control. |
Estimation Syntax: The formula
argument takes syntax designed
to be a similar as possible to lme4
. That is, one can specify models
using y ~ x + (1 | g)
where (1 | g)
indicates a random intercept. While
not tested extensively, terms of (1 | g / f)
should work as expected. Terms
of (1 + x || g)
may work, although will raise a warning about duplicated
names of random effects. (1 + x || g)
terms may not work with spline
estimation. To get around this, one can might copy the column g
to
g_copy
and then write (1 | g) + (0 + x | g_copy)
.
Splines: Splines can be added using the term v_s(x)
for a
spline on the variable x
. These are transformed into hierarchical
terms in a standard fashion (e.g. Ruppert et al. 2003) and then estimated
using the variational algorithms. At the present, only truncated linear
functions (type = "tpf"
; the default) and O'Sullivan splines (Wand and
Ormerod 2008) are included. The options are described in more detail at
v_s.
It is possible to have the spline vary across some categorical predictor by
specifying the "by"
argument such as v_s(x, by = g)
. In effect,
this adds additional hierarchical terms for the group-level deviations from
the "global" spline. Note: In contrast to the typical presentation of
these splines interacted with categorical variables (e.g., Ruppert et al.
2003), the default use of "by"
includes the lower order interactions
that are regularized, i.e. (1 + x | g)
, versus their unregularized
version (e.g., x * g
); this can be changed using the by_re
argument described in v_s. Further, all group-level deviations from
the global spline share the same smoothing parameter (same prior
distribution).
Default Settings: By default, the model is estimated using the
"strong" (i.e. fully factorized) variational assumption. Setting
vglmer_control(factorization_method = "weak")
will improve the quality
of the variance approximation but may take considerably more time to
estimate. See Goplerud (2022a) for discussion.
By default, the prior on each random effect variance (\Sigma_j
) uses a Huang-Wand prior (Huang
and Wand 2013) with hyper-parameters \nu_j = 2
and A_{j,k} = 5
.
This is designed to be proper but weakly informative. Other options are
discussed in vglmer_control under the prior_variance
argument.
By default, estimation is accelerated using SQUAREM (Varadhan and Roland
2008) and (one-step-late) parameter expansion for variational Bayes. Under
the default "strong"
factorization, a "translation" expansion is used;
under other factorizations a "mean" expansion is used. These can be adjusted
using vglmer_control. See Goplerud (2022b) for more discussion of
these methods.
This returns an object of class vglmer
. The available methods
(e.g. coef
) can be found using methods(class="vglmer")
.
Contains the estimated distribution of the fixed effects
(\beta
). It is multivariate normal. mean
contains the means;
var
contains the variance matrix; decomp_var
contains a matrix
L
such that L^T L
equals the full variance matrix.
Contains the estimated distribution of the random effects
(\alpha
). They are all multivariate normal. mean
contains the
means; dia.var
contains the variance of each random effect. var
contains the variance matrix of each random effect (j,g). decomp_var
contains a matrix L
such that L^T L
equals the full variance of
the entire set of random effects.
If factorization_method="weak"
, this is a list with one
element (decomp_var
) that contains a matrix L
such that L^T
L
equals the full variance matrix between the fixed and random effects
q(\beta,\alpha)
. The marginal variances are included in beta
and
alpha
. If the factorization method is not "weak"
, this is
NULL
.
Contains the estimated distribution of each random
effect covariance \Sigma_j
; all distributions are Inverse-Wishart.
cov
contains a list of the estimated scale matrices. df
contains a list of the degrees of freedom.
If a Huang-Wand prior is used (see Huang and Wand 2013 or Goplerud
2022b for more details), then the estimated distribution. Otherwise, it is
NULL
. All distributions are Inverse-Gamma. a
contains a list of
the scale parameters. b
contains a list of the shape parameters.
If family="linear"
, this contains a list of the
estimated parameters for \sigma^2
; its distribution is Inverse-Gamma.
a
contains the scale parameter; b
contains the shape
parameter.
If family="negbin"
, this contains the variational
parameters for the log dispersion parameter \ln(r)
. mu
contains
the mean; sigma
contains the variance.
Family of outcome.
Contains the ELBO at the termination of the algorithm.
data.frame
tracking the ELBO per iteration.
Contains the control parameters from vglmer_control
used in estimation.
Variety of internal parameters used in post-estimation functions.
Contains the formula used for estimation; contains the
original formula, fixed effects, and random effects parts separately for
post-estimation functions. See formula.vglmer
for more details.
Goplerud, Max. 2022a. "Fast and Accurate Estimation of Non-Nested Binomial Hierarchical Models Using Variational Inference." Bayesian Analysis. 17(2): 623-650.
Goplerud, Max. 2022b. "Re-Evaluating Machine Learning for MRP Given the Comparable Performance of (Deep) Hierarchical Models." Working paper.
Huang, Alan, and Matthew P. Wand. 2013. "Simple Marginally Noninformative Prior Distributions for Covariance Matrices." Bayesian Analysis. 8(2):439-452.
Ruppert, David, Matt P. Wand, and Raymond J. Carroll. 2003. Semiparametric Regression. Cambridge University Press.
Varadhan, Ravi, and Christophe Roland. 2008. "Simple and Globally Convergent Methods for Accelerating the Convergence of any EM Algorithm." Scandinavian Journal of Statistics. 35(2): 335-353.
Wand, Matt P. and Ormerod, John T. 2008. "On Semiparametric Regression with O'Sullivan Penalized Splines". Australian & New Zealand Journal of Statistics. 50(2): 179-198.
set.seed(234)
sim_data <- data.frame(
x = rnorm(100),
y = rbinom(100, 1, 0.5),
g = sample(letters, 100, replace = TRUE)
)
# Run with defaults
est_vglmer <- vglmer(y ~ x + (x | g), data = sim_data, family = "binomial")
# Simple prediction
predict(est_vglmer, newdata = sim_data)
# Summarize results
summary(est_vglmer)
# Extract parameters
coef(est_vglmer); vcov(est_vglmer)
# Comparability with lme4,
# although ranef is formatted differently.
ranef(est_vglmer); fixef(est_vglmer)
# Run with weaker (i.e. better) approximation
vglmer(y ~ x + (x | g),
data = sim_data,
control = vglmer_control(factorization_method = "weak"),
family = "binomial")
# Use a spline on x with a linear outcome
vglmer(y ~ v_s(x),
data = sim_data,
family = "linear")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.