vglmer {vglmer}

Variational Inference for Hierarchical Generalized Linear Models

Description

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.

Usage

vglmer(formula, data, family, control = vglmer_control())

Arguments

formula

lme4 style-formula for random effects. Typically, (1 + z | g) indicates a random effect for each level of variable "g" with a differing slope for the effect of variable "z" and an intercept (1); see "Details" for further discussion and how to incorporate splines.

data

data.frame containing the outcome and predictors.

family

Options are "binomial", "linear", or "negbin" (experimental). If "binomial", outcome must be either binary (\{0,1\}) or cbind(success, failure) as per standard glm(er) syntax. Non-integer values are permitted for binomial if force_whole is set to FALSE in vglmer_control.

control

Adjust internal options for estimation. Must use an object created by vglmer_control.

Details

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.

Value

This returns an object of class vglmer. The available methods (e.g. coef) can be found using methods(class="vglmer").

beta

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.

alpha

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.

joint

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.

sigma

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.

hw

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.

sigmasq

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.

ln_r

If family="negbin", this contains the variational parameters for the log dispersion parameter \ln(r). mu contains the mean; sigma contains the variance.

family

Family of outcome.

ELBO

Contains the ELBO at the termination of the algorithm.

ELBO_trajectory

data.frame tracking the ELBO per iteration.

control

Contains the control parameters from vglmer_control used in estimation.

internal_parameters

Variety of internal parameters used in post-estimation functions.

formula

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.

References

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.

Examples

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")