mhglm: Fitting Moment Hierarchical Generalized Linear Models
In mbest: Moment-Based Estimation for Hierarchical Models

Description Usage Arguments Details Value Note References See Also Examples

mhglm is used to fit a moment hierarchical generalized linear model of one level. mhglm_ml is used to fit a moment hierarchical generalized linear model of arbitrary number of levels (including one level).

mhglm(formula, family = gaussian, data, weights, subset,
      na.action, start = NULL, etastart, mustart, offset,
      control = list(), model = TRUE, method = "mhglm.fit",
      x = FALSE, z = FALSE, y = TRUE, group = TRUE,
      contrasts = NULL)

mhglm.fit(x, z, y, group, weights = rep(1, nobs),
          start = NULL, etastart = NULL, mustart = NULL,
          offset = rep(0, nobs), family = gaussian(), 
          control = list(), intercept = TRUE, dispersion = NULL)

mhglm_ml(formula, family = gaussian, data, weights, subset,
         na.action, start = NULL, etastart, mustart, offset,
         control = list(), model = TRUE, method = "mhglm_ml.fit",
         x = FALSE, z = FALSE, y = TRUE, group = TRUE,
         contrasts = NULL)

mhglm_ml.fit(x, z, y, group, weights = rep(1, nobs),
             start = NULL, etastart = NULL, mustart = NULL,
             offset = rep(0, nobs), family = gaussian(), 
             control = list(), intercept = TRUE)

`formula, family, data, weights, subset, na.action, start, etastart, mustart, offset, model, contrasts, intercept`	These arguments are analogous to the similarly-named arguments for the `glm` and `glm.fit` functions.
`control`	a list of parameters for controlling the fitting process. For `mhglm.fit` this is passed to `mhglm.control`.
`method`	the method to be used in fitting the model. The default method `"mhglm.fit"` uses moment-based estimates; the alternative `"model.frame"` returns the model frame and does no fitting.
`x, z, y, group`	For `mhglm` and `mhglm_ml`: logical values indicating whether the response vector, model matrices, and grouping factor used in the fitting process should be returned as components of the returned value. For `mhglm.fit`: `x` is a fixed effect design matrix of dimension `n * p`, `z` is a random effect design matrix of dimension `n * q`, `y` is a vector of observations of length `n`, and `group` is a grouping factor vector of length `n`. For `mhglm_ml.fit`: `x` is a fixed effect design matrix of dimension `n * p`, `z` is a list of L elements, with L the depth of nested hierarchies, each element of `z` is a random effect design matrix of dimension `n * q_i`, with q_i the feature dimension on tree depth i, `y` is a vector of observations of length `n`, and `group` is a list of L elements (same L as `z`), each element of `group` is a grouping factor vector of length `n`.
`dispersion`	If `NULL`, will estimate from data; otherwise use this argument as dispersion parameter.

These functions are analogues of glm and glm.fit, meant to be used for fitting hierarchical generalized linear models. A typical predictor has the form response ~ terms + (reterms | group) where response is the (numeric) response vector, terms is a series of terms which specifies a linear predictor for response, reterms is a series of terms with random coefficients (effects), and group is a grouping factor; observations with the same grouping factor share the same random effects.

mhglm and mhglm.fit only allow one random effect term, along with a single level of hierarchy. mghlm_ml and mhglm_ml.fit allow multiple random effect terms so long as levels of random effects are hierarchically nested. If the random effect design matrices are the same for each level, a predictor has the form response ~ terms + (reterms | g1/.../gQ). If the random effects design matrices differ from level to level, colons are used to delineate the nesting structure; for example, response ~ fe + (re1 | g1) + (re2 | g2:g1) + (re3 | g3:g2:g1).

mhglm allows || in the formula response ~ terms + (reterms || group) to indicate that random effects are independent, that is the random effects covariance matrix has non-zero value only on its diagonal. mhglm_ml currently does not support ||, to indicate indpenedent random effects, set control=list(diagcov = TRUE).

mhglm returns an object of class inheriting from "mhglm". mhglm_ml returns an object of class inheriting from "mhglm_ml".

The function summary can be used to obtain or print a summary of the results.

The generic accessor functions fixef, ranef, VarCorr, sigma, fitted.values and residuals can be used to extract various useful features of the value returned by mhglm or mhglm_ml.

If the moment-based random effect covariance is not positive-semidefinite, then a warning will be issued, and a projection of the estimate to the positive-semidefinite cone will be used instead.

P. O. Perry (2017) "Fast moment-based estimation for hierarchical models."

N. Zhang, K. Schmaus, and P. O. Perry (2018) "Fitting deeply-nested hierarchical models to a large book review dataset using moment-based estimators."

terms.mhglm, model.matrix.mhglm, and predict.mhglm for mhglm methods, and the generic functions fitted.values, residuals, summary, vcov, and weights.

Generic functions fixef, ranef, VarCorr, and sigma for features related to mixed effect models.

glmer (package lme4) for fitting generalized linear mixed models with likelihood-based estimates.

library(lme4)
## The following examples are adapted from lme4:
(fm1 <- mhglm(Reaction ~ Days + (Days | Subject), gaussian, sleepstudy))

(fm2 <- mhglm(Reaction ~ Days + (Days || Subject), gaussian, sleepstudy))

(gm <- mhglm(cbind(incidence, size - incidence) ~ period + (1 | herd),
             data = cbpp, family = binomial))

## The following examples are for multilevel models
g_data <- mhglm_sim(n = 30, m_per_level = c(10, 5, 2), sd_intercept = c(1, 1, 1),
                    sd_slope = c(1, 1, 1), family = "gaussian", seed = 12345)

(m1 <- mhglm_ml(y ~ 1 + x + (1 + x | g1/g2/g3), gaussian, g_data))
# or equivalent
(m2 <- mhglm_ml(y ~ 1 + x + (1 + x | g1) + (1 + x | g2:g1) + (1 + x | g3:g2:g1),
                gaussian, g_data))