glmm | R Documentation |

This function fits generalized linear mixed models (GLMMs) by approximating the likelihood with ordinary Monte Carlo, then maximizing the approximated likelihood.

glmm(fixed, random, varcomps.names, data, family.glmm, m, varcomps.equal, weights=NULL, doPQL = TRUE,debug=FALSE, p1=1/3,p2=1/3, p3=1/3, rmax=1000,iterlim=1000, par.init, zeta=5, cluster=NULL)

`fixed` |
an object of class " |

`random` |
an object of class " |

`varcomps.names` |
The names of the distinct variance components in order of |

`data` |
an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which |

`family.glmm` |
The name of the family. Must be class |

`m` |
The desired Monte Carlo sample size. See a note under "Details." |

`varcomps.equal` |
An optional vector with elements 1 through the number of distinct variance components. Denotes variance components are to be set equal by assigning them the same integer. The length of varcomps.equal must be equal to the length of the list of random effects formulas. If omitted, varcomps.equal assumes no variance component should be set equal. |

`weights` |
An optional vector with length equal to the length of the response vector. This argument makes the specified observations more or less informative for the model. See a note under "Details." |

`doPQL` |
logical. If |

`debug` |
logical. If |

`p1` |
A probability for mixing the random effects generated from three distributions. |

`p2` |
A probability for mixing the random effects generated from three distributions. |

`p3` |
A probability for mixing the random effects generated from three distributions. |

`rmax` |
The maximum allowed trust region radius. This may be set very large. If set small, the algorithm traces a steepest ascent path. This is an argument for |

`iterlim` |
A positive integer specifying the maximum number of trust iterations to be performed before the trust program is terminated. This is an argument for |

`par.init` |
An optional argument. A single vector that specifies the initial values of the fixed effects and variance components. The parameters should be inputted in the order that |

`zeta` |
A scalar that specifies the degrees of freedom for the t-distribution from which random effects are generated. |

`cluster` |
An optional argument. A cluster created by the user to enable computations to be done in parallel. See "Details" for more information on the default value. |

Let *beta* be a vector of fixed effects and let *u* be a vector of random effects. Let *X* and *Z* be design matrices for the fixed and random effects, respectively. The random effects are assumed to be normally distributed with mean 0 and variance matrix *D*, where *D* is diagonal with entries from the unknown vector *nu*. Letting *g* be the link function, *g(mu)=X beta+ ZU*. If the response type is Bernoulli or Binomial, then the logit function is the link; if the response type is Poisson, then the natural logarithm is the link function.

Models for glmm are specified symbolically. A typical fixed effects model has the form `response ~ terms`

where `response`

is the (numeric) response vector and `terms`

is a series of terms which specifies a linear predictor for response. A terms specification of the form `first + second`

indicates all the terms in `first`

together with all the terms in `second`

with duplicates removed.

A specification of the form `first:second`

indicates the set of terms obtained by taking the interactions of all terms in `first`

with all terms in `second`

. The specification `first*second`

indicates the cross of `first`

and `second`

. This is the same as `first + second + first:second`

.

The terms in the formula will be re-ordered so that main effects come first, followed by the interactions, all second-order, all third-order and so on: to avoid this, pass a `terms`

object as the formula.

If you choose `binomial.glmm`

as the `family.glmm`

, then your response should be a two-column matrix: the first column reports the number of successes and the second reports the number of failures.

The random effects for glmm are also specified symbolically. The random effects model specification is typically a list. Each element of the `random`

list has the form `response ~ 0 + term`

. The 0 centers the random effects at 0. If you want your random effects to have a nonzero mean, then include that term in the fixed effects. Each variance component must have its own formula in the list.

To set some variance components equal to one another, use the `varcomps.equal`

argument. The argument `varcomps.equal`

should be a vector whose length is equal to the length of the random effects list. The vector should contain positive integers, and the first element of the `varcomps.equal`

should be 1. To set variance components equal to one another, assign the same integer to the corresponding elements of `varcomps.equal`

. For example, to set the first and second variance components equal to each other, the first two elements of `varcomps.equal`

should be 1. If `varcomps.equal`

is omitted, then the variance components are assumed to be distinct.

Each distinct variance component should have a name. The length of `varcomps.names`

should be equal to the number of distinct variance components. If `varcomps.equal`

is omitted, then the length of `varcomps.names`

should be equal to the length of `random`

.

The package uses a relevance-weighted log density weighting scheme, similar to that described in Hu and Zidek (1997). This means that the `weights`

argument functions in the following manner: a model built for a data set with three observations (obs. 1, obs. 2, obs. 3) and weighting scheme 1, 1, 2 gives the same results as the model built for a data set with four observations (obs. 1, obs. 2, obs. 3, obs. 3 - the last two observations are identical) and weighting scheme 1, 1, 1, 1.

Monte Carlo likelihood approximation relies on an importance sampling distribution. Though infinitely many importance sampling distributions should yield the correct MCMLEs eventually, the importance sampling distribution used in this package was chosen to reduce the computation cost. When `doPQL`

is `TRUE`

, the importance sampling distribution relies on PQL estimates (as calculated in this package). When `doPQL`

is `FALSE`

, the random effect estimates in the distribution are taken to be 0, the fixed effect estimates are taken to be 0, and the variance component estimates are taken to be 1.

This package's importance sampling distribution is a mixture of three distributions: a t centered at 0 with scale matrix determined by the PQL estimates of the variance components and with `zeta`

degrees of freedom, a normal distribution centered at the PQL estimates of the random effects and with a variance matrix containing the PQL estimates of the variance components, and a normal distribution centered at the PQL estimates of the random effects and with a variance matrix based on the Hessian of the penalized log likelihood. The first component is included to guarantee the gradient of the MCLA has a central limit theorem. The second component is included to mirror our best guess of the distribution of the random effects. The third component is included so that the numerator and the denominator are similar when calculating the MCLA value.

The Monte Carlo sample size `m`

should be chosen as large as possible. You may want to run the model a couple times to begin to understand the variability inherent to Monte Carlo. There are no hard and fast rules for choosing `m`

, and more research is needed on this area. For a general idea, I believe the `BoothHobert`

model produces stable enough estimates at *m=10^3* and the `salamander`

model produces stable enough estimates at *m=10^5*, as long as `doPQL`

is `TRUE`

.

To decrease the computation time involved with larger Monte Carlo sample sizes, a parallel computing component has been added to this package. Generally, the larger the Monte Carlo sample size `m`

, the greater the benefit of utilizing additional cores. By default, `glmm`

will create a cluster that uses a single core. This forces all computations to be done sequentially rather than simultaneously.

To see the summary of the model, use summary().

`glmm`

returns an object of class `glmm`

is a list containing at least the following components:

`beta ` |
A vector of the Monte Carlo maximum likelihood estimates (MCMLEs) for the fixed effects. |

`nu ` |
A vector of the Monte Carlo maximum likelihood estimates for the variance components. |

`loglike.value` |
The Monte Carlo log likelihood evaluated at the MCMLEs |

`loglike.gradient` |
The Monte Carlo log likelihood gradient vector at the MCMLEs |

`loglike.hessian` |
The Monte Carlo log likelihood Hessian matrix at the MCMLEs |

`mod.mcml` |
A list containing the weighted (if applicable) fixed effect design matrix, the list of weighted (if applicable) random effect design matrices, the weighted (if applicable) response, and the number of trials (for the Binomial family). |

`call` |
The call. |

`fixedcall` |
The fixed effects call. |

`randcall` |
The random effects call. |

`x` |
The unweighted design matrix for the fixed effects. |

`y` |
The unweighted response vector. |

`z` |
The unweighted design matrix for the random effects. |

`weights` |
The weights vector. |

`family.glmm` |
The name of the family. Must be class |

`varcomps.names` |
The vector of names for the distinct variance components. |

`varcomps.equal` |
The vector denoting equal variance components. |

`umat` |
A matrix with |

`pvec` |
A vector containing |

`beta.pql` |
PQL estimate of |

`nu.pql` |
PQL estimate of |

`u.pql` |
PQL predictions of the random effects. |

`zeta` |
The number of degrees of freedom used in the t component of the importance sampling distribution. |

`cluster` |
A cluster created by the user for use during the approximation of the log-likelihood value, gradient and hessian. |

`debug` |
If |

The function `summary`

(i.e., `summary.glmm`

) can
be used to obtain or print a summary of the results. The generic accessor function `coef`

(i.e., `coef.glmm`

) can be used to extract the coefficients.

Christina Knudson

Geyer, C. J. (1994) On the convergence of Monte Carlo maximum likelihood calculations. *Journal of the Royal Statistical Society, Series B*, **61**, 261–274.
doi: 10.1111/j.2517-6161.1994.tb01976.x.

Geyer, C. J. and Thompson, E. (1992) Constrained Monte Carlo maximum likelihood for dependent data. *Journal of the Royal Statistical Society, Series B*, **54**, 657–699.
doi: 10.1111/j.2517-6161.1992.tb01443.x.

Knudson, C. (2016). Monte Carlo likelihood approximation for generalized linear mixed models. PhD thesis, University of Minnesota. http://hdl.handle.net/11299/178948

Knudson, C., Benson, S., Geyer, J., Jones, G. (2021) Likelihood-based inference for generalized linear mixed models: Inference with the R package glmm. *Stat*, 10:e339. doi: 10.1002/sta4.339.

Sung, Y. J. and Geyer, C. J. (2007) Monte Carlo likelihood inference for missing data models. *Annals of Statistics*, **35**, 990–1011.
doi: 10.1214/009053606000001389.

Hu, F. and Zidek, J. V. (2001) The relevance weighted likelihood with
applications. In: Ahmed, S. E. and Reid N. (eds).
*Empirical Bayes and Likelihood Inference*.
Lecture Notes in Statistics, **148**, 211–235. Springer, New York.
doi: 10.1007/978-1-4613-0141-7_13.

#First, using the basic Booth and Hobert dataset #to fit a glmm with a logistic link, one variance component, #one fixed effect, and an intercept of 0. The Monte Carlo #sample size is 100 to save time. library(glmm) data(BoothHobert) set.seed(1234) mod.mcml1 <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE) mod.mcml1$beta mod.mcml1$nu summary(mod.mcml1) coef(mod.mcml1) # This next example has crossed random effects but we assume # all random effects have the same variance. The Monte Carlo #sample size is 100 to save time. #data(salamander) #set.seed(1234) #onenu <- glmm(Mate~Cross, random=list(~0+Female,~0+Male), #varcomps.names=c("only"), varcomps.equal=c(1,1), data=salamander, #family.glmm=bernoulli.glmm, m=100, debug=TRUE, doPQL=TRUE) #summary(onenu)

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.