R: Mixed-effect modeling functions.

gremlin {gremlin}

R Documentation

Mixed-effect modeling functions.

Description

Fit and setup functions for linear mixed-effect model (Gaussian responses).

Usage

gremlin(
  formula,
  random = NULL,
  rcov = ~units,
  data = NULL,
  ginverse = NULL,
  Gstart = NULL,
  Rstart = NULL,
  Bp = NULL,
  Gcon = NULL,
  Rcon = NULL,
  maxit = 20,
  algit = NULL,
  vit = 1,
  v = 1,
  control = gremlinControl(),
  ...
)

gremlinR(
  formula,
  random = NULL,
  rcov = ~units,
  data = NULL,
  ginverse = NULL,
  Gstart = NULL,
  Rstart = NULL,
  Bp = NULL,
  Gcon = NULL,
  Rcon = NULL,
  maxit = 20,
  algit = NULL,
  vit = 1,
  v = 1,
  control = gremlinControl(),
  ...
)

## S3 method for class 'gremlin'
getCall(x, ...)

## S3 method for class 'gremlin'
update(object, ...)

gremlinSetup(
  formula,
  random = NULL,
  rcov = ~units,
  data = NULL,
  ginverse = NULL,
  Gstart = NULL,
  Rstart = NULL,
  Bp = NULL,
  Gcon = NULL,
  Rcon = NULL,
  maxit = 20,
  algit = NULL,
  vit = 1,
  v = 1,
  control = gremlinControl(),
  ...
)

## S3 method for class 'gremlin'
is(x)

## S3 method for class 'grMod'
is(x)

mkModMats(
  formula,
  random = NULL,
  rcov = ~units,
  data = NULL,
  subset = NULL,
  ginverse = NULL,
  na.action = na.pass,
  offset = NULL,
  contrasts = NULL,
  Xsparse = TRUE,
  ...
)

Arguments

`formula`	A `formula` for the response variable(s) and fixed effects.
`random`	A `formula` for the random effects.
`rcov`	A `formula` for the residual covariance structure.
`data`	A `data.frame` in which to look for the terms in `formula`, `random`, and `rcov`.
`ginverse`	A `list` of (preferably sparse) inverse matrices that are proportional to the covariance structure of the random effects. The name of each element in the list should match a column in `data` that is associated with a random term. All levels of the random term should appear as `rownames` for the matrices.
`Gstart`	A `list` of matrices with starting (co)variance values for the G-structure or random effects terms.
`Rstart`	A `list` of matrices with starting (co)variance values for the R-structure or residual terms.
`Bp`	A prior specification for fixed effects.
`Gcon`, `Rcon`	A `list` of matrices with constraint codes for the G-structure/random effects or R-structure/residual effects terms, respectively. Must be a `character` of either `"F"` for fixed, `"P"` for positive, or `"U"` for unbounded.
`maxit`	An `integer` specifying the maximum number of likelihood iterations.
`algit`	A `character` vector of length 1 or more or an expression to be evaluated that specifies the algorithm to use for proposing (co)variances in the next likelihood iteration.
`vit`	An `integer` value specifying the verbosity of screen output on each iteration. A value of zero gives no iteration specific output and larger values increase the amount of information printed on the screen.
`v`	An `integer` value specifying the verbosity of screen output regarding the model fitting process. A value of zero gives no details and larger values increase the amount of information printed on the screen.
`control`	A `list` returned by `gremlinControl` containing specific named values from that function. See `gremlinControl`.
`...`	Additional arguments to be passed to control the model fitting.
`x`, `object`	An object of class `gremlin`.
`subset`	An expression for the subset of `data` to use.
`na.action`	What to do with NAs.
`offset`	Should an offset be specified.
`contrasts`	Specify the type of contrasts for the fixed effects.
`Xsparse`	Should sparse matrices be used for the fixed effects design matrix.

Value

A list containing an object of class grMod and, if a model was fit (gremlin or gremlinR) then an object containing details of the REML iterations (object itMat). An object of class grMod contains:

call

The model call.

modMats

A list of the model matrices used to construct the mixed model equations.

y

The response vector.

ny

The number of responses.

ncy

The number of columns of the original response.

X

The fixed effects design matrix.

nb

The number of columns in X.

Zr

The residual design matrix.

Zg

A list of the design matrices for each random term.

nG

The number of parameters in the G structure.

listGeninv

A list of generalized inverse matrices.

logDetG

The log-determinants of the generalized inverse matrices - necessary to calculate the log-likelihood.

rfxIncContrib2loglik

A numeric value containing the sum of the log determinants of the random effects that do not change between log-likelihood iterations (i.e., the part of the log determinants of (co)variance matrices to be estimated that have been factored out).

ndgeninv

A logical indicating which terms in the random formula have generalized inverses associated with them (non-diagonal matrices in the Kronecker product.

dimsZg, nminffx, rfxlvls, nminfrfx

Numeric vectors or scalars describing the numbers of random effects or some function of random and fixed effects.

conv, bounds

(Co)variance component constraints and boundaries of the allowable parameter space for each component.

thetav

A vector of the (co)variance parameters to be estimated by REML withe the attribute “skel” giving the skeleton to recreate a list of matrices from this vector.

thetaG, thetaR

Vectors indexing the random and residual (co)variances, respectively, in a list of (co)variance matrices (i.e., theta).

nu

A list of transformed (co)variance matrices to be fit by REML. If a residual variance has been factored out of the mixed model equations, nu contains the ‘lambda’ parameterization with expresses the (co)variance components as ratios of variance parameters with the residual variance. The ‘nu’ scale (co)variances are the ones actually fit by REML.

sigma2e

The estimate of the factored out residual variance from the mixed model equations (i.e., the ‘lambda’ scale) \sigma^{2}_{e}.

p

An integer for the total number of (co)variances to be estimated.

lambda

A logical indicating whether the ‘lambda’ scale parameterization has been used.

uni

A logical to indicate if the model is univariate or not.

W, tWW, RHS, Bpinv

Sparse matrices of class Matrix that form the mixed model equations and do not change between iterations of REML. These are the column binded ‘X’ and ‘Z’ design matrices for fixed and random effects, the cross-product of W, the Right-Hand Side of the mixed model equations, and the inverse of the fixed effect prior matrix (zeroes on the diagonal if no priors have been specified). Note, these may be NULL if lambda=FALSE, because the NULL objects are not used or do change between REML iterations.

sLc

A Matrix containing the symbolic Cholesky factorization of the coefficient matrix of the Mixed Model Equations.

sln

A one column matrix of solutions in the mixed model equations.

Cinv_ii

A one column matrix of variances for the solutions to the mixed model equations. These are obtained from the diagonal of the inverse Coefficient matrix in the mixed model equations. If lambda is TRUE then these are on the lambda scale and must be multiplied by sigma2e to be converted to the original data scale.

r

A one column matrix of residual deviations, response minus the values expected based on the solutions, corresponding to the order in modMats$y.

AI

A matrix of values containing the Average Information matrix, or second partial derivatives of the likelihood with respect to the transformed (co)variance components (‘nu’). The inverse of this matrix gives the sampling (co)variances of these transformed (co)variance components.

dLdnu

A single column matrix of first derivatives of the transformed (co)variance parameters (‘nu’) with respect to the log-Likelihood.

maxit

See the parameter described above.

algit

A character vector of REML algorithms to use in each iteration.

vit

See the parameter described above.

v

See the parameter described above.

cctol

A numeric vector of convergence criteria thresholds. See gremlinControl for more details.

ezero, einf

numeric values for the effective numbers to use as “zero” and maximum negative or positive numbers. Values less than ezero are treated as zero and fixed to this value. Values less than -1*einf or greater than einf are restricted to these values. See gremlinControl for more details.

step

A numeric value indicating the step-reduction to use. See gremlinControl for more details.

itMat

A matrix of details about each iteration. Rows indicate each REML iteration (rownames reflect the REML algorithm used) and columns contain:

nu, theta: (Co)variance parameters.
sigma2e: See ‘sigma2e’ described above.
tyPy, logDetC: Estimates for two these two components of the log of the REML likelihoods. These are obtained from Cholesky factorization of the coefficient matrix of the mixed model equations.
loglik: The REML log-likelihood.
itTime: Time elapsed for each REML iteration.

Functions

mkModMats: Generates model matrices.

Author(s)

matthewwolak@gmail.com

Examples

  grSire <- gremlin(WWG11 ~ sex, random = ~ sire, data = Mrode11)
  # Now drop sire random effects and use the `anova` method to compare models
  grLM <- update(grSire, random = ~ 1)  #<-- use `~1` to drop all random effects
    anova(grSire, grLM)

  # Modular functions
  ## get model matrices for a mixed model
  mM11 <- mkModMats(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11)

  ## setup model, but do not evaluate the log-likelihood
  grSetup <- gremlinSetup(WWG11 ~ sex - 1, random = ~ sire, data = Mrode11)
  ## maximize the restricted maximum likelihood
  grOut <- remlIt(grSetup)
  summary(grOut)

## Not run: 
  # Following the example from Mrode 2005, chapter 11.
  library(nadiv)  #<-- to construct inverse of the numerator relatedness matrix
  Ainv <- makeAinv(Mrode11[, 1:3])$Ainv

  gr11lmm <- gremlin(WWG11 ~ sex - 1,
random = ~ calf,
data = Mrode11,
ginverse = list(calf = Ainv),
Gstart = matrix(0.2), Rstart = matrix(0.4),  #<-- specify starting values
maxit = 15,    #<-- maximum iterations
     v = 2, vit = 1,  #<-- moderate screen output (`v`) every iteration (`vit`)
     algit = "AI")  #<-- only use Average Information algorithm iterations
  summary(gr11lmm)

  # Compare the model to a Linear Model with no random effects
  ## Use `update()` to update the model
  gr11lm <- update(gr11lmm, random = ~ 1)  #<-- `~1`=drop all random effects
  summary(gr11lm)

  # Do analysis of variance between the two models
  ## See AIC or evaluate likelihood ratio against a Chi-squared distribution
  anova(gr11lm, gr11lmm)

## End(Not run)

[Package gremlin version 1.0.1 Index]