asreml {asreml} | R Documentation

## Fit the linear mixed model.

### Description

asreml estimates variance components under a general linear mixed model by residual maximum likelihood (REML).

### Usage

    asreml(
      fixed = y ~ 1,
      random = ~NULL,
      sparse = ~NULL,
      residual = ~NULL,
      G.param = list(),
      R.param = list(),
      data = sys.parent(),
      na.action = na.method(),
      subset,
      weights,
      predict = predict.asreml(),
      vcm = vcm.lm(),
      vcc = matrix(NA),
      family = asr_gaussian(),
      asmv = NULL,
      mbf = list(),
      group = list(),
      equate.levels = character(0),
      start.values = FALSE,
      knot.points = list(),
      pwr.points = list(),
      wald = list(),
      prune = list(),
      combine = list(),
      uid = list(),
      mef = list(),
      last = list(),
      model.frame = TRUE,
      ...
    )

### Arguments

`fixed` | A formula object specifying the fixed terms in the model, with the response on the left of a ~ operator, and the terms, separated by + operators, on the right. If `data` is given, all names used in all formulae should appear in the data frame. A model with the intercept as the only fixed effect can be specified as ~1; there must be at least one fixed effect specified. If the response (y) evaluates to a matrix then a factor `trait` with levels `dimnames(y)[[2]]` is added to the model frame, and must be explicitly included in the model formulae.
`random` | A formula object specifying the random effects in the model. This argument has the same general characteristics as `fixed`, but there can be no left side to the ~ operator. Variance structures imposed on random terms are specified using special model functions.
`sparse` | A formula object, specifying the fixed effects for which the full variance-covariance matrix is not required. This argument has the same general characteristics as `fixed`, but there can be no left side to the ~ expression. Wald statistics are not available for sparse fixed terms in order to reduce the computing load.
`residual` | A formula object specifying the residual model; any term specified on the left of the ~ expression is ignored. The default is ~`units`, where the reserved word `units` is defined as `seq(1,nrow(data))` and is automatically included in the model frame. Variance models for the residual component of the model can be specified using special model functions. For single-section univariate models, the residual variance model determines the computational mode: If the residual variance model specifies a correlation structure (includes `id()`), then the model is fitted on the gamma scale, otherwise the model is fitted on the sigma scale. The default is `id(units)` if not explicitly specified.
`G.param` | Either, a list object derived from the `random` formula, holding initial parameter estimates and boundary constraints for each term, or a character string naming a comma delimited file with a header line and three columns for the variance component name, initial value and constraint code, respectively. This file can be created using the `start.values` argument; the internal list object is then generated from the contents of this file. On termination, `G.param` is updated with the final `random` component estimates.
`R.param` | Either, a list object derived from the `random` formula, holding initial parameter estimates and boundary constraints, or a character string naming a comma delimited file with a header line and three columns for the variance component name, initial value and constraint code, respectively. This file can be created using the `start.values` argument; the internal list object is then generated from the contents of this file. On termination, `R.param` is updated with the final `residual` component estimates.
`data` | A data frame in which to interpret the variables named in `fixed`, `random`, `sparse` and `residual`. If the `data` argument is missing, the default is `sys.parent()`. The data frame is converted internally to a `data.table` object and returned as such by `model.frame`.
`na.action` | A call to `na.method()` specifying the action to be taken when missing values are encountered in the response (`y`) or explanatory variables (`x`). The function definition for `na.method` is: `function(y=c("include","omit","fail"),` `x=c("fail","include","omit"))` The default action is to include (and estimate) missing values in the response, and raise an error if there are missing values in the explanatory variables.
`subset` | A logical vector identifying which subset of the rows of `data` should be used in the fit. All observations are included by default.
`weights` | A character string or name identifying the column of `data` to use as weights in the fit.
`predict` | A list object specifying the classifying factors and related options when forming predictions from the model. This list would normally be the value returned by a call to the method `predict` for `asreml` objects.
`vcm` | A matrix defining relationships among variance parameters. The matrix has a row for each original variance parameter and a column for each new parameter. The default is the identity matrix, that is, no action. See `vcm.lm` for further information and an example.
`vcc` | Equality constraints between variance parameters; a two-column numeric matrix with a `dimnames` attribute. The first column defines the grouping structure of equated components, that is, components within an equality group are given the same numeric index, and the second column contains the scaling coefficients. The `dimnames()[[1]]` attribute must match the component names in the `asreml` parameter vector; see `start.values`. The parameters are scaled relative to the first parameter in its group, so the scaling of the first parameter in each group is one. For example, the following `vcc` matrix
1 | 1
2 | 1
2 | 2
3 | 1
is equivalent to the `vcm` matrix
1 | 0 | 0
0 | 1 | 0
0 | 2 | 0
0 | 0 | 1
`family` | A list of functions and expressions for defining the link and variance functions. Optionally a list of such structures for a multivariate analysis involving non-normal variates. Currently this is restricted to a bivariate model where the first variate (excluding the multinomial distribution) is non-normal. Supported families are gaussian, inverse Gaussian, binomial, negative binomial, poisson, Gamma and multinomial. Family objects are generated from the `asreml` family functions which prefix the usual function names with `"asr_"`; for example `asr_gaussian()`, `asr_binomial()` etc. In addition to the link argument, these functions take an additional `dispersion` argument and a `total` argument where relevent; for example: `asr_binomial(dispersion=1.0, total=counts)`. The default for `asr_gaussian()` is `dispersion=NA`, which implies that `asreml` will estimate the dispersion parameter, otherwise the scale is fixed at the nominated value.
`asmv` | A character string or name specifying the column in the data that identifies the traits in a multivariate analysis. If not `NULL`, `asmv` implies that the data for a multivariate analysis is set up as though it were for a univariate analysis with the response in a single variate.
`mbf` | A named list specifying sets of covariates to be included with one or more `mbf()` model functions. Each component of the list must in turn contain components named `key` and `cov`, where `cov` is a character string naming the data frame holding the covariates, and `key` is a character vector of length 2 naming the columns in `data` and `cov`, respectively, used to match corresponding records in the two data frames. The default is an empty list.
`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. The component names can then appear in asreml model formulae using the `grp()` special function. The default is an empty list.
`equate.levels` | A character vector of factor names whose levels are to be equated. If factor `A` has levels a,b,c,d and factor `B` has levels a,b,c,e, the effect of `equate.levels(A, B)` is that both `A` and `B` have 5 levels, with `as.numeric(A)` = 1,2,3,4 and `as.numeric(B)` = 1,2,3,5. This may be necessary if using the `and` model function to overlay columns of the model's design matrix in forming a compound term. The default is a zero length character vector.
`start.values` | If `TRUE`, `asreml` exits prior to the fitting process and returns a list of length 3: the `G.param` and `R.param` lists, and a data frame containing variance parameter names, initial values and boundary constraints. Initial values or constraints can then be set in the list or data frame objects. If a character string, then a file of that name is created and the data frame object containing initial parameter values is written out in comma separated form. This file can be edited externally and subsequently specified in the `G.param` or `R.param` arguments.
`knot.points` | A named list where each component is a vector of user supplied knot points for a particular spline term; the component name is the object of the `spl()` model function.
`pwr.points` | A named list with each component containing a vector of distances to be used in a one-dimensional power model. The component names must correspond to the `object` arguments of the power function model terms.
`wald` | A named list with four components: `denDF`, `ssType`, `Ftest` and `kenadj`. denDF a character string from the set `["none"`, `"numeric"`, `"algebraic"`, `"default"]` specifying the calculation of approximate denominator degrees of freedom. The default `"none"` is to suppress the computations. Algebraic computations are not feasible in large analyses, use `"default"` to autommatically choose numeric or algebraic computations depending on problem size. The denominator degrees of freedom are calculated according to Kenward and Roger, 1997 for terms in the `fixed` model formula. ssType can be `"incremental"` for incremental sum of squares (the default) or `"conditional"` for F-tests that respect both structural and intrinsic marginality. Ftest a one sided formula of the form \sim `test_term | background_terms` specifying a conditional Wald test of the contribution of `test_term` conditional on those fixed terms listed in `background_terms`, and the those in the `random` and `sparse` model formulae. kenadj a character string from the set `["none"`, `"expected"`, `"observed]"`, specifying any adjustment to the variance matrix for fixed effects.
`prune` | A named list with each component generated from a call to `Subset()`. `prune`, in conjunction with `Subset` and the model function `sbs()`, forms a new factor from an existing one by selecting a subset of its levels. The function `Subset` is defined as: `function(f, x)` where `f` is the name of an existing factor and `x` is a character or numeric vector of levels to select. The name of the list component is the new factor that may appear in the model formulae as the argument to the `sbs()` model function. For example, `prune=list(A=Subset(Site, c(2,3)))` creates a new factor `A` by selecting the second and third levels of `Site`, and would be included in the model as `sbs(A)`. While the actions of `prune` can be duplicated outside `asreml`, `sbs()` is necessary if the `asreml` method `predict()` is to be used.
`combine` | A named list with each component generated from a call to `Levels()`. `combine`, in conjunction with `Levels` and the model function `gpf()` forms a new factor from an existing one by merging a subset of its levels. The function `Levels` is defined as: `function(f, x)` where `f` is the name of an existing factor and `x` is a vector of length `length(levels(f))` defining the levels of `f` to merge. The name of the list component is the new factor that may appear in the model formulae as the argument to the `gpf()` model function. For example, if `Site` has levels `"1"`, `"2"` and `"3"`, `combine=list(A=Levels(Site, c("1","2","1")))` creates a new factor `A` with levels `"1"` and `"2"` by merging levels `"1"` and `"3"` of `Site`, and would be included in the model as `gpf(A)`. While the actions of `combine` can be duplicated outside `asreml`, `gpf()` is necessary if the `asreml` method `predict()` is to be used.
`uid` | A named list with each component generated from a call to `Units()`. `uid`, in conjunction with `Units` and the model function `uni()` forms a new factor by selecting a subset of records for an existing one. The function `Units` is defined as: `function(f, n=0)` where `f` is the name of an existing factor and `n` is a character or numeric scalar that determines which records are selected. The default, `n=0`, forms a factor with a level for each record where `f` is non-zero (strictly, `f` != 0). Otherwise, a factor with a level for each record in `data` where `f` has the value `n` is formed. For example, `uid = list(A = Units(group, 1))` creates a new factor `A` with levels from `row.names(data)` where `group` = 1, and would be included in the model as `uni(A)`. While the actions of `uid` can be duplicated outside `asreml`, `uni()` is necessary if the `asreml` method `predict()` is to be used.
`mef` | A named list linking a relationship matrix (or its inverse) as specified in the `vm()` special function with the original matrix of subject x regressor (typically molecular marker) scores. If not an empty list (the default), `mef` flags the computation of the regressor (marker) effects from the subject effects. For example, `mef=list(MM = snp.mat)` links the relationship matrix `MM` to the original marker scores in `snp.mat`. The `mef` list would typically be set from a call to the `asreml` `meff()` method.
`last` | A named list restricting the order equations are solved in the sparse partition for the nominated model terms. Each component of the list is named by a model term and contains a scalar n specifying that the first n levels of the term be solved after all others in the sparse set. It is intended for use when there are multiple fixed terms in the sparse equations so that `asreml` will be consistent in which efects are identified as singular. A maximum of three factor/level pairs can be specified.
`model.frame` | If `TRUE` (the default) the model frame (a `data.table` object with additional attributes derived from the model specification) is included in the returned object. The model frame is required by the `asreml` `summary`, `plot`, `resid` and `fitted` methods. In large analyses, the model frame is likely to be a large object. If `model.frame` is a character string, the model frame is saved in a file as an `RDS` object by a call to `saveRDS()`, and named by the supplied string with the extension `.RDS`. If the model frame is not included in the returned `asreml` object, this `RDS` file is searched for by the methods noted above.
`...` | Additional arguments to `asreml` from `asreml.options()`: `maxit`, `workspace`, `pworkspace`, `fixgammas`, `trace`, `aom`.

### Details

Models for `asreml` are specified symbolically in the formula objects `fixed, random, sparse` and `residual`. A typical model has the form `response` \sim `terms`, `fixed` only, or \sim `terms` for `random`, `sparse` and `residual`, where `response` is the (usually numeric) response vector and `terms` is a linear predictor for `response`. An exception is raised if the response is a factor and `family` is not multinomial.

The formulae objects are parsed in the context of the data frame, all internal data structures are constructed in R or compiled code, and the model is fitted by calls to the underlying `Fortran` `REML` routines (Gilmour et al., 1995). Variance models for random model terms are specified using special functions in the `random` and `residual` formulae. If not specified, the variance models default to (scaled) identity structures. A table of special model functions is included below; see the reference guide or appropriate vignette for further details and examples of their use. Some of these model functions require the formula arguments to be partially evaluated before the final model frame is computed; it is recommended that all names used in the formulae be resolvable in a data frame named by the `data` argument.

If the response is a matrix, a multivariate linear model is fitted to the columns unless `family` = `asr_multinomial()` is declared.

The terms in the `fixed` formula are re-ordered by default so that main effects preceed interactions in increasing order. The option `keep.order` (see `asreml.options`) can be used to modify this behaviour.

A formula has an implied intercept term. To remove the intercept use `y` ~ `-1 + ...`. This is only effective in the `fixed` formula; in all other formula arguments any reference to the intercept is ignored. Note that currently there must be at least one fixed effect in the model.

In addition to the formal arguments, various options can be set with `asreml.options`; these are stored in an environment for the duration of the `R` session.

`asreml` uses either a `"gamma"` (ratio) or `"sigma"` (component) scale parameterization for estimation depending on the residual model specification. The current default for single section analyses is the gamma parameterization if the error model specifies a correlation structure. In this case, all scale parameters are estimated as a ratio with respect to the residual variance, with correlation parameters unchanged. If the residual model specifies a variance structure then variance parameters are estimated on the sigma scale. For models with more than one residual section, `asreml` always estimates variance parameters on the sigma scale.

### Value

An object of class `asreml` containing the results of the fitted model. Instances of generic methods such as `plot()`, `predict()` and `summary()` return various derived results of the fit; `resid()`, `coef()` and `fitted()` extract some of its components. See `asreml.object` for the components of the returned list.

### Special functions

Special model functions are used in `asreml` formulae objects to create new or modify existing model terms, or more often to specify the variance model associated with one or more terms. These functions can be broadly categorised as constructor-type functions, or (default) identity, time-series, general-structure, 1d metric, 2d metric, known relationship, general variance structures that span more than one term (str), or user-defined structures.

The special model functions that are available in `asreml` model formulae are introduced below; see the user guide, relevent vignette or the man pages for selected functions for more details or illustration.

The symbols used in the following tables are defined as:

obj | a factor in `data`.
n | `length(levels(obj))`.
p | number of parameters estimated by the base function.
v | number of parameters estimated by the homogeneous function form.
h | number of parameters estimated by the heterogeneous function form.

Constructor type functions

Call | Description
con(obj) | Apply sum-to-zero constraints to factor `obj`.
C(obj, contr) | Define contrasts among the levels of `obj` from the coefficients in `contr`.
lin(obj) | Fits factor `obj` as a variate.
pow(obj, p, offset) | Create the term `(offset+obj)`^p.
pol(x, t) | Orthogonal polynomials to degree t; -t omits the intercept polynomial.
leg(x, t) | Legendre polynomials to degree t; -t omits the intercept polynomial.
spl(x, k) | The random component of a cubic spline; optionally k knot points.
dev(x) | Fit variate `x` as a factor; typically used for spline deviations.
ma(obj) | Forms a moving average (1) design matrix from factor `obj`
at(obj, vec) | Form conditioning covariables for the levels in `obj` given in `vec`.
dsum(~term | obj, ~model) | Direct sum of `term` for the levels of `obj` with variance structure `model`.
| Used in `residual` to define multiple sections.
and(obj, k) | Add k times the design matrix for `obj` to the previous columns.
grp(name) | Include the term defined by `name` in the `group` argument in the model.
mbf(name) | Include the covariates defined by `name` in the `mbf` argument as a factor.
sbs(name) | Include the term defined in the `prune` argument in the model.
gpf(name) | Include the term defined in the `combine` argument in the model.
uni(name) | Include the term defined in the `uid` argument in the model.

Default identity

Call | Description | p | v | h
id(obj) | identity | 0 | 1 | n

Time series type models

Call | Description | p | v | h
ar1(obj) | autoregressive order 1 | 1 | 2 | 1+n
ar2(obj) | autoregressive order 2 | 2 | 3 | 2+n
ar3(obj) | autoregressive order 3 | 3 | 4 | 3+n
sar(obj) | symmetric autoregressive | 1 | 2 | 1+n
sar2(obj) | symmetric autoregressive order 2 | 2 | 3 | 2+n
ma1(obj) | moving average order 1 | 1 | 2 | 1+n
ma2(obj) | moving average order 2 | 2 | 3 | 2+n
arma(obj) | autoregressive-moving average | 2 | 3 | 2+n

General structure models

Call | Description | p | v | h
cor(obj) | simple correlation | 1 | 2 | 1+n
corb(obj, b) | banded correlation; b bands | b | b+1 | b+n
corg(obj) | general correlation | n(n-1)/2 | 1+n(n-1)/2 | n(n+1)/2
diag(obj) | heterogeneous variance | n |  |
us(obj) | unstructured variance | n(n+1)/2 |  |
sfa(obj, k) | factor analytic; k factors | kn+n |  |
fa(obj, k) | sparse factor analytic | kn+n |  |
facv(obj, k) | factor analytic, covariance form | kn+n |  |
rr(obj, k) | reduced rank variant of fa | kn+n |  |
chol(obj, k) | Cholesky order k | (k+1)(n-k/2) |  |
cholc(obj, k) | Cholesky | (k+1)(n-k/2) |  |
ante(obj, k) | antedependence order k | (k+1)(n-k/2) |  |
mthr(obj) | multinomial models only |  |

Metric based models in 1D or 2D

Call | Description | p | v | h
exp(x) | exponential 1D | 1 | 2 | 1+n
iexp(x, y) | isotropic exponential 2D | 1 | 2 | 1+n
aexp(x, y) | anisotropic exponential 2D | 2 | 3 | 2+n
gau(x) | gaussian 1D | 1 | 2 | 1+n
igau(x, y) | isotropic gaussian 2D | 1 | 2 | 1+n
agau(x, y) | anisotropic gaussian 2D | 2 | 3 | 2+n
ieuc(x, y) | isotropic euclidean 2D | 1 | 2 | 1+n
isp(x, y) | isotropic spherical 2D | 1 | 2 | 1+n
cir(x, y) | isotropic circular 2D | 1 | 2 | 1+n
mtrn(x, y, ...) | Matern class 2D | * | * | *

^*See the user guide for an extended discussion of the Matern class.

Known relationship structures

Call | Description
vm(obj, source, singG) | Create a term based on `obj` with known variance structure in `source`.
ide(obj) | Identity term based on `obj` with levels as for `vm(obj)`.

General variance structures

Call | Description
str(~terms, ~model) | Apply the direct product variance structure in `~model`
| to the set of terms in `~terms`.

User defined structures

Call | Description
own(obj, fun, init, type) | Call `fun` with the parameter estimates in `init`
| to compute the variance matrix and its derivatives.

### References

Gilmour AR, Thompson R, Cullis BR (1995). “AI, An Efficient Algorithm for REML Estimation in Linear Mixed Models.” Biometrics, 51, 1440–1450. Kenward MG, Roger JH (1997). “The Precision of Fixed Effects Estimates from Restricted Maximum Likelihood.” Biometrics, 53, 983–997.

### See Also

`asreml.options` `asreml.object` `asr_families`

### Examples

    ## Not run:
    data(oats)
    oats.asr <- asreml(yield ~ Variety*Nitrogen, random = ~ Blocks/Wplots, data=oats)

    ## End(Not run)

[Package asreml version 4.1.0.176 Index]