LMMsolve: Solve Linear Mixed Models
In LMMsolver: Linear Mixed Model Solver

LMMsolve

R Documentation

Solve Linear Mixed Models

Description

Solve Linear Mixed Models using REML.

Usage

LMMsolve(
  fixed,
  random = NULL,
  spline = NULL,
  group = NULL,
  ginverse = NULL,
  weights = NULL,
  data,
  residual = NULL,
  family = gaussian(),
  offset = 0,
  tolerance = 1e-06,
  trace = FALSE,
  maxit = 250,
  theta = NULL,
  grpTheta = NULL
)

Arguments

`fixed`	A formula for the fixed part of the model. Should be of the form "response ~ pred"
`random`	A formula for the random part of the model. Should be of the form "~ pred".
`spline`	A formula for the spline part of the model. Should be of the form "~ spl1D()", ~ spl2D()" or "~spl3D()". Generalized Additive Models (GAMs) can also be used, for example "~ spl1D() + spl2D()"
`group`	A named list where each component is a numeric vector specifying contiguous fields in data that are to be considered as a single term.
`ginverse`	A named list with each component a symmetric matrix, the precision matrix of a corresponding random term in the model. The row and column order of the precision matrices should match the order of the levels of the corresponding factor in the data.
`weights`	A character string identifying the column of data to use as relative weights in the fit. Default value NULL, weights are all equal to one.
`data`	A data.frame containing the modeling data.
`residual`	A formula for the residual part of the model. Should be of the form "~ pred".
`family`	An object of class family specifying the distribution and link function.
`offset`	An optional numerical vector containing an a priori known component to be included in the linear predictor during fitting.
`tolerance`	A numerical value. The convergence tolerance for the modified Henderson algorithm to estimate the variance components.
`trace`	Should the progress of the algorithm be printed? Default `trace = FALSE`.
`maxit`	A numerical value. The maximum number of iterations for the algorithm. Default `maxit = 250`.
`theta`	initial values for penalty or precision parameters. Default `NULL`, all precision parameters set equal to 1.
`grpTheta`	a vector to give components the same penalty. Default `NULL`, all components have a separate penalty.

Details

A Linear Mixed Model (LMM) has the form

y = X \beta + Z u + e, u ~ N(0,G), e ~ N(0,R)

where y is a vector of observations, \beta is a vector with the fixed effects, u is a vector with the random effects, and e a vector of random residuals. X and Z are design matrices.

LMMsolve can fit models where the matrices G^{-1} and R^{-1} are a linear combination of precision matrices Q_{G,i} and Q_{R,i}:

G^{-1} = \sum_{i} \psi_i Q_{G,i} \;, R^{-1} = \sum_{i} \phi_i Q_{R,i}

where the precision parameters \psi_i and \phi_i are estimated using REML. For most standard mixed models 1/{\psi_i} are the variance components and 1/{\phi_i} the residual variances. We use a formulation in terms of precision parameters to allow for non-standard mixed models using tensor product splines.

Value

An object of class LMMsolve representing the fitted model. See LMMsolveObject for a full description of the components in this object.

Examples

## Fit models on john.alpha data from agridat package.
data(john.alpha, package = "agridat")

## Fit simple model with only fixed effects.
LMM1 <- LMMsolve(fixed = yield ~ rep + gen,
                data = john.alpha)

## Fit the same model with genotype as random effect.
LMM1_rand <- LMMsolve(fixed = yield ~ rep,
                     random = ~gen,
                     data = john.alpha)

## Fit the model with a 1-dimensional spline at the plot level.
LMM1_spline <- LMMsolve(fixed = yield ~ rep + gen,
                       spline = ~spl1D(x = plot, nseg = 20),
                       data = john.alpha)

## Fit models on multipop data included in the package.
data(multipop)

## The residual variances for the two populations can be different.
## Allow for heterogeneous residual variances using the residual argument.
LMM2 <- LMMsolve(fixed = pheno ~ cross,
                residual = ~cross,
                data = multipop)

## QTL-probabilities are defined by the columns pA, pB, pC.
## They can be included in the random part of the model by specifying the
## group argument and using grp() in the random part.

# Define groups by specifying columns in data corresponding to groups in a list.
# Name used in grp() should match names specified in list.
lGrp <- list(QTL = 3:5)
LMM2_group <- LMMsolve(fixed = pheno ~ cross,
                      group = lGrp,
                      random = ~grp(QTL),
                      residual = ~cross,
                      data = multipop)

LMMsolver documentation built on Sept. 11, 2024, 7:31 p.m.