Introducing the linear mixed model

Introduction

The two main uses of the linear mixed model are

to estimate treatment effects and test for differences, and

to predict genetic merit of individuals.

In the first case, treatment effects are estimated as fixed effects and reflect what happened in the experiment. Random effects then reflect sources of variation fitted to get the most efficient hypothesis tests.

In the second case, the random genetic effects are of most interest, having been adjusted for other fixed or random sources of variation.

If y denotes the n × 1 vector of observations, the linear mixed model can be written as

y = X τ + Z u + e

where τ is the p × 1 vector of fixed effects, X is an n × p design matrix of full column rank (after deleting any redundant singular columns) which associates observations with the appropriate combination of fixed effects, u is the q × 1 vector of random effects, Z is the n × q design matrix which associates observations with the appropriate combination of random effects, and e is the n × 1 vector of residual errors.

The data may be collected from a controlled experiment, or by sampling an existing population.

Introducing the linear mixed model

The linear mixed model is also called a linear mixed effects model. It is assumed that
u ~ N(0, G (σ_g)
is independent of (uncorrelated with)
e ~ N(0, R(σ_r)
where the matrices G and R are functions of parameters σ_g and σ_r, respectively.

The variance matrix for y is then of the form
var ( y) = Z G (σ_g)Z' + R (σ_r)
which we will refer to as the sigma parameterization of the G and R variance structures, and the variance structure parameters in σ_g and σ_r will be referred to as sigmas although they are not necessarily all variances. The variance models given by G and R are referred to as G structures and R structures respectively.

We illustrate these concepts using the simplest linear mixed model.

A one-way classification

Consider a one-way classification comprising a single random effect u, and a residual error term e. The two random components of this model, namely u and e, are each assumed to be independent and identically distributed (IID) and to follow a normal distribution such that u~ N(0, σ²_uI_q) and e~ N(0, σ²_eI_n). Hence the variance of u has the form
var (y) = σ²_uZZ' + σ²_eI_n

This model has two variance structure parameters or sigmas: the variance component σ²_u associated with u, and the variance component σ²_eassociated with e. Mapping this equation back to the general expression, we have σ_g=σ²_u, G(σ_g) = σ²_uI_q, σ_r =σ²_e and R(σ_r) = σ²_eI_n .

G structure for the random model terms}

Typically, τ and u are composed of several model terms, that is, τ can be partitioned as τ=[τ'₁: τ'_t]' and u can be partitioned as u = [u'₁: u'_b]' with X and Z partitioned conformably as X = [X₁: X_t] and Z = [Z₁: Z_b].

For u thus partitioned, we impose a direct sum structure on the matrix G, written G = ⊕_i=1:b G_i giving

G₁	0	:	0	0
0	G₂	:	0	0
:	:	:	:	:
0	0	:	G_b-1	0
0	0	:	0	G_b

where ⊕ is the direct sum operator, each G_i is of size q_i and q = ∑_iq_i.

The default assumption is that each random model term generates one component of this direct sum; var(u_i) = G_i for i = 1 : b. In this case, the random effects from any two distinct model terms are uncorrelated.

Traditional Variance components models

Extending the one-way example to a linear mixed model with more than one (b>1) random terms (typically known as a variance components mixed model), the random effects u_i in u, and the residual errors e, are assumed uncorrelated and to each be normally distributed with mean zero and variance given by var(u_i)= σ²_{u_i} I_{q_i} and var(e) = σ²_eI_n where I_{q_i} and I_n are identity matrices of dimension q_i and n, respectively. In this case
var(y)= ∑_i=1:b σ²_{u_i} Z_iZ'_i + σ²_eI_n

Partitioning the residual error term

As for the fixed and random model terms, it is often appropriate to partition the vector of residual errors e according to some conditioning factor. We use the term section to describe this partitioning. A common reason for using sections in e is because subsets of the residuals have different variances. For example, in the analysis of multi-environment trials (METs) it is natural to expect that each trial will have a different (possibly spatial) error structure. In this case, for s sections we have e = [e'₁: e'_s]' assuming that the data vector is ordered by section, and where e_j represents the vector of errors for the j^th section.

For e partitioned as e = [e₁,e₂, :,e_s]' we allow the matrix R to have a similar direct sum structure, with R = ⊕_j=1:s R_j giving

R₁	0	:	0	0
0	R₂	:	0	0
:	:	:	:	:
0	0	:	R_s-1	0
0	0	:	0	R_s

for s sections and the data ordered by section. It may be necessary to re-order (sort) the data units to achieve this structure.

While the default variance structure for a section is IID (σ²_{e_j} I_{n_j} ), a more complex structure is often appropriate. ASReml offers a wide range of variance structures.

Gamma parameterization for the linear mixed model

In the sigma parameterization of model (2.3) the scale parameters in both G(σ_g and R(σ_r are variances. However, in a univariate analysis with one data section, the one residual variance parameter may be factored out of G(σ_g) and R(σ_r) so that
var (y)= σ²_eΣ_i=1:b γ_{u_i} Z_iZ'_i + I_n
for γ_{u_i} being the ratio of the variance component for the random term u_i relative to error variance, that is, γ_{u_i} = σ²_{u_i} /σ²_e. We will refer to this as the gamma parameterization, and the individual variance structure parameters in γ_g and γ_r will be referred to as gammas. ASReml defaults to the gamma parameterizations when there is just one residual variance as it it easier to set initial values in terms of these ratios. This parameterization was used in some of the early development of REML for the traditional mixed model.

This leads to an alternative specification of the random effects:
u ~ N(0, σ²_eG (σ_g)) and e ~ N(0, σ²_eR(σ_r))
where γ_g and γ_r are the parameters for the scaled (by σ²_e) G and R matrices, respectively; R(γ_r) is actually a correlation (or identity) matrix. The variance matrix for y is then
var(y)=σ²_eZ G(γ_g)Z' + R(γ_r).

Parameter types

Each sigma in σ_g and σ_r and each gamma in γ_g and γ_r has a parameter type, for example, variance, variance ratio, autocorrelation parameter, factor loading, distance. Furthermore, the parameters in σ_g, σ_r, γ_g and γ_r can have different types. For example, the spatial analysis of a simple column trial would involve a variance or variance ratio and a spatial autocorrelation parameter.

Variance structures for the random model terms

The random model terms u_i in u define the random effects and associated design matrices, Z_i∈Z , but a variance structure is required for each term before the model can be fitted. If nothing is specified, σ²_iI_{q_i} is assumed. Other structures are specified by using functions applied to the model term, or to components of the model term in an interaction, to define G_i. This produces a consolidated model term that simultaneously defines both the design matrix (Z_i) and variance model (G_i).

A random model term may comprise either a single factor or be an interaction of several component factors to give a compound model term. Consider a compound model term represented by A.B, where the component factors A and B have m and n levels respectively and the ˙ . operator forms a model term with levels corresponding to the combinations of all levels of A with all levels of B. The effects ab_ij for A.B are generated with the levels of B nested in the levels of A, i.e. the levels of B cycling fastest: (ab)= (ab₁₁, ab₁₂, : ab_1n, ab₂₁, ab₂₂, : ab_2n, : ab_m1, ab_m2, : ab_mn)'.

If we propose variance structures A = [A_ij] for A and B = [B_kl] for B, the variance structure for A.B is given by the direct product A⊗B. This means that the covariance between two effects ab_ik and ab_jl in (ab) is constructed as the product of the covariance between a_i and a_j ( A_ij ) and the covariance between b_k and b_l (B_kl ) giving cov(ab_ik , ab_jl )= A_ij B_kl .

Note that for direct product structures, only one of the component 'variance structure' matrices should contain estimable variances. The other can be an Identity, a correlation matrix or some other known matrix with no estimable variance parameter; you can't estimate the two elements of a product individually if you only know the product.

Direct product structures

Mathematically, the direct product of two matrices A^{(m × p)} and B^{(n × q)} is A⊗B giving

	A_₁₁B	:	A_1p B
:	:	:
	A_{_m1}B	:	A_{_mp}B

Variance structures formed as direct product construction are known as separable variance structures and we call the assumption that a separable variance structure is plausible the assumption of separability.

A small direct product structure

If A has 3 levels and B has 2 levels, then the term A.B would have the 6 levels: (ab)= (ab₁₁, ab₁₂, ab₂₁, ab₂₂, ab₃₁, ab₃₂)'. Writing Σ_A and Σ_B as

A₁₁	A₁₂	A₁₃
A₂₁	A₂₂	A₂₃	and	B₁₁	B₁₂
A₃₁	A₃₂	A₃₃		B₂₁	B₂₂

then cov(ab₂₁ , ab₃₂ ) = A₂₃ B₁₂

Direct products in R structures

Consider a vector of common errors associated with an experiment. The usual least squares assumption (and the default in ASReml) is that these are independently and identically distributed (IID). However, if e was from a field experiment laid out in a rectangular array of r rows by c columns, we could arrange the residuals as a matrix and might consider that they were autocorrelated within rows and columns. Writing the residuals as a vector in field order, that is, by sorting the residuals rows within columns (plots within blocks) the variance of the residuals might then be
σ²_eΣ_r(ρ_r)⊗ Σ_c(ρ_c)
where Σ_r(ρ_r) and Σ_c(ρ_c) are correlation matrices for the row model (order r, autocorrelation parameter ρ_r) and column model (order c, autocorrelation parameter ρ_c), respectively. More specifically, a two-dimensional separable autoregressive spatial structure (AR1 ⊗ AR1) is sometimes assumed for the common errors in a field trial analysis (see Gogel (1997) and Cullis {em et al}. (1998) for examples). In this case Σ_r and Σ_c are like

ρ_r

ρ_c

ρ_r²

ρ_r

and

ρ_c²

ρ_c

ρ_r^r-1

ρ_r^r-2

ρ_r^{r - 3}

ρ_c^c-1

ρ_c^c-2

ρ_c^{c - 3}

Alternatively, the residuals might relate to a multivariate analysis with n_t traits and n units and be ordered traits within units. In this case an appropriate variance structure might be I_n ⊗ Σ where Σ^{(n_t× n_t} is a general or unstructured variance matrix.

Direct products in G structures

Likewise, the random terms in u in the model may have a direct product variance structure. For example, for a field trial with s sites, g varieties and the effects ordered varieties within sites, the model term site.variety, assumed random, may have the variance structure Σ ⊗ I_g where Σ is the variance matrix for sites. This would imply that the varieties are independent random effects within each site, have different variances at each site, and are correlated across sites. Whenever a random term is formed as the interaction of two factors you should consider whether the IID assumption is sufficient or if a direct product structure might be more appropriate.

Common consolidated model terms

This section lists the common variance structures used in ASReml. Some common variables used in mixed models are:
range/row/column: index field plots in field order
site/location/expt: associate plots in the same place
animal/sire/dam/genotype/entry genetically related individuals

ar1(row).ar1v(column) The ar1() variance function defines an autoregressive correlation matrix with one correlation parameter. Appending v to the name scales the matrix by a second variance parameter. This model term, typically used for the residual in the spatial analysis of a field trial, has a row autocorrelation, a column autocorrelation and a variance parameter. When analysing data from several sites together, the model term becomes sat(site).ar1(row).ar1v(column) which is expanded to define a separate model term for each site.

units.us(Trait) units is an internally defined factor indexing the rows of the data file. Trait is an internally defined factor indexing the variables involved in a multivariate analysis. The us() variance function has a variance parameter for each trait/variable and a covariance for each pairwise combination; ten parameters for four traits.

nrm(animal) The nrm() variance function associates a known (previously defined) Numerator Relationship Matrix (NRM) with the animal factor. When used alone, a variance parameter, the genetic variance, is estimated. The grm() variance function associates a known (previously defined) Genomic Relationship Matrix (GRM) with its argument.

us(Trait).nrm(animal) estimates a genetic variance for each trait, and a genetic covariance for each pair of traits.

coruh(site).genotype The coru() variance function defines a uniform correlation matrix. The h estimates a separate genetic variance for each site, to scale the correlation matrix. In this model term, genotypes are treated as independent.

xfa1(site).grm(genotype) The xfa1() variance function defines a factor analytic variance structure (ΓΛ'+Ψ) with 1 factor (Γ has 1 column, called loadings, Ψ is diagonal with values called specific variances). This is an intermediate model having more parameters than coruh() but less than us(). It is also a sparser model to fit than coruh() or us() when there are many sites (more than 4). The model term also includes genomic relationships among genotypes in the second dimension of the site.genotype table.

rrd1(site).grm(genotype) This is expanded to rr1(site).grm(genotype) + diag(site).grm(genotype). The rr1() variance function is the xfa1() variance function without specific variances; it estimates the covariances. The diag() variance function has all covariances zero. So the sum is equivalent to the regular xfa1() variance function. However, the mixed model equations have a sparser form when interacted with the dense grm(genotype) component and so runs faster, especially for many sites.

Return to index