Agricultural Statistics & Data Science · Design and analysis of experiments · also Variance Component Estimation, Mixed Model Analysis, REML and BLUP
Models with more than one source of random variation: ANOVA, ML and REML estimation of variance components, BLUE and BLUP, Henderson's mixed-model equations, and spatial models that correct field trials for fertility trends and neighbour correlation.
In the exams: AICE-JRF Ag Stat Unit 4 'random, fixed and mixed effects' (present); AICE-JRF Ag Stat Unit 3 genetic variance (present); PhD course in most programmes
Field, animal and laboratory data rarely meet the assumptions of the classical analysis of variance exactly: blocks are random samples of the field, varieties are tested in unbalanced sets, measurements are repeated over time, plots next to each other are correlated and sites differ in variability. Linear mixed models handle all of this within one framework. This subject covers mixed models with fixed and random effects, the ANOVA (Henderson) methods of estimating variance components, maximum likelihood and REML, Henderson's mixed-model equations with BLUE and BLUP, denominator degrees of freedom for unbalanced data (Satterthwaite and Kenward-Roger), covariance structures for repeated measures, spatial analysis of field trials with nearest-neighbour and AR1 x AR1 models, spline-based trend adjustment such as SpATS, model selection by likelihood-ratio tests and information criteria, and the main software.
It builds on the linear models subject (general linear model, generalised least squares, fixed, random and nested effects and expected mean squares) and on the incomplete block designs subject (recovery of inter-block information, which is a mixed-model analysis with random blocks). Its methods are used directly in the multi-environment trials subject, the statistical genetics subjects (heritability and breeding values) and wherever repeated measurements occur.
In examinations, ICAR AICE-JRF/SRF Agricultural Statistics names random, fixed and mixed effects and nested and crossed classifications in Unit 4 and genetic variance and its partitioning in Unit 3; variance component estimation, REML and BLUP are standard in PhD courses and in ASRB scientist interviews. Questions test expected mean squares and ANOVA estimators, the difference between ML and REML, the form of the mixed-model equations and shrinkage, degrees-of-freedom approximations, the number of parameters in covariance structures, the AR1 x AR1 correlation, and boundary-corrected likelihood-ratio tests.
The linear models subject introduced fixed and random effects and the one-way random model. A linear mixed model writes all of this in one equation: y = X beta + Z u + e, where beta is a vector of fixed effects with design matrix X, u is a vector of random effects with design matrix Z, and e is the residual. The random parts are assumed normal with mean zero, u ~ N(0, G) and e ~ N(0, R), independent of each other. The variance of y is then V = Z G Z' + R. Everything that makes field and animal data awkward fits into G and R. Blocks of an alpha design, sires in a progeny test, locations sampled from a region and plants within plots go into u with variance components in G. Correlation between neighbouring plots, between repeated measurements on the same animal, or unequal variances between sites go into R. The fixed effects are the treatments or contrasts the experiment was designed to compare. Ordinary ANOVA is the special case G = 0 (or random effects handled through expected mean squares) and R = sigma squared I.
An effect is treated as random when its levels can be regarded as a sample from a population and inference is wanted about that population, or when there are many levels each with little information, so that borrowing strength across levels helps. Blocks are usually random in incomplete block designs, which recovers inter-block information as described in the incomplete block designs subject. Genotypes are fixed when the aim is to compare named varieties for release, and random when the aim is to estimate genetic variance or heritability or to select the best lines from a large breeding population, where best linear unbiased prediction (BLUP) shrinks noisy estimates towards the mean. The decision changes what is estimated and how: fixed effects give best linear unbiased estimates (BLUEs) of their levels; random effects give variance components and BLUPs. It also changes the tests, because the appropriate error for a fixed effect depends on which other effects are random.
An RCBD with random blocks is y = mu + treatment + block + e, with blocks in Z and treatments in X; with complete blocks the treatment BLUEs equal the ordinary treatment means, and the mixed model only changes the interpretation of the block variance. A split plot is y = mu + A + B + AB + rep + rep.A + e, where the random rep.A term (main plots) carries what the classical analysis calls error (a); the mixed model reproduces the two error terms automatically and keeps working when plots are missing. An alpha design is y = mu + genotype + rep + rep.block + e with random rep.block, which recovers inter-block information. A progeny test is y = mu + site + sire + e with random sires, whose variance gives a quarter of the additive genetic variance in a half-sib design. Seeing designs this way makes clear that mixed models are not an exotic extra but the natural statement of most agricultural experiments, and that the classical ANOVA tables of the design subjects are special cases that hold when the data are balanced.
Counts, proportions and scores are often analysed with generalised linear mixed models (GLMMs), which combine the link functions of generalised linear models with random effects. Germination counts per pot with random blocks, disease scores with random sites, and the number of piglets born with random sows are typical. GLMMs are estimated by approximate likelihood methods (Laplace approximation, adaptive quadrature or penalised quasi-likelihood), and their variance components are on the link scale. They are introduced in the categorical data subject and are mentioned here so that the analyst knows that a normal mixed model is not the only option when residuals are clearly non-normal.
The oldest way to estimate variance components is to compute an ANOVA, write the expected value of each mean square as a linear combination of the variance components, and solve the equations obtained by setting the observed mean squares equal to their expectations. In a balanced one-way random model with a groups of n, E(MSA) = sigma squared + n sigma_a squared and E(MSE) = sigma squared, so sigma squared-hat = MSE and sigma_a squared-hat = (MSA - MSE)/n. For balanced data these estimators are unbiased, have minimum variance among quadratic unbiased estimators under normality, and coincide with REML estimates whenever they are non-negative.
With unbalanced data the expectations of sums of squares depend on the pattern of subclass numbers. Henderson (1953) set out three methods. Method 1 uses sums of squares computed as in balanced analyses but with unequal numbers; it is easy but valid only for random models (no fixed effects other than the mean). Method 2 first adjusts the data for estimated fixed effects and then applies Method 1; it fails when fixed and random factors interact. Method 3, the fitting-constants method, uses reductions in sums of squares from fitting sub-models, and is valid for mixed models, though not unique when several orders of fitting are possible. All ANOVA-type estimators share two problems: they can be negative, because they are differences of mean squares, and with unbalanced data they lack the optimality they have under balance. A negative estimate usually signals a small component and is often truncated to zero, but truncation biases the estimate. These problems led to likelihood-based methods.
In a balanced two-stage nested random model with a = 10 sires, n = 8 progeny each, suppose MS(sires) = 36 and MS(within) = 12. Then sigma_s squared-hat = (36 - 12)/8 = 3 and sigma_w squared-hat = 12, so the intraclass correlation is 3/15 = 0.2. In a half-sib design the sire component estimates one quarter of the additive genetic variance, so heritability is estimated as 4 x 3/15 = 0.8, a figure that would be suspiciously high for most traits and would prompt a check of the data and design. If MS(sires) had been 10, below MS(within), the ANOVA estimate of sigma_s squared would be -0.25, which would usually be reported as zero, with a note.
A mixed model can be read in two ways. Conditionally, given the random effects, the observations are independent with mean X beta + Z u and variance R. Marginally, averaging over the random effects, the observations have mean X beta and variance V = ZGZ' + R, so observations sharing a random effect are correlated. In a randomised block with random blocks, the marginal correlation between two plots in the same block is sigma_b squared / (sigma_b squared + sigma squared). Fixed effects are estimated from the marginal model by generalised least squares; random effects are predicted from the conditional model. Keeping both views in mind clarifies why random blocks induce correlation and why BLUPs depend on V.
Under normality the log-likelihood of y in the mixed model depends on beta and the variance parameters through V. Maximum likelihood (ML), developed for mixed models by Hartley and Rao (1967), maximises it over all parameters jointly. ML estimates of variance components are non-negative (when the parameter space is constrained) and work for any data pattern, but they are biased downward because they ignore the degrees of freedom used to estimate the fixed effects. In the simplest case of a sample of n observations with mean unknown, the ML estimate of the variance is the sum of squares divided by n, not n - 1.
Patterson and Thompson (1971) proposed restricted (or residual) maximum likelihood (REML), which maximises the likelihood of a set of error contrasts, linear combinations of the data whose distribution does not depend on the fixed effects. REML accounts for the loss of degrees of freedom: in the simple case it divides the sum of squares by n - 1, and in balanced ANOVA designs it reproduces the ANOVA estimators when they are non-negative. REML is now the standard method for estimating variance components in agricultural and animal research and is the default in most mixed-model software. Both ML and REML are computed iteratively, by Fisher scoring, expectation-maximisation, derivative-free methods or the average-information algorithm used in ASReml, which converges quickly for large problems. Convergence failures and estimates on the boundary (a component estimated as zero) are common with small data sets and over-complex models and should be read as signals that the model is asking more than the data can tell.
Because the REML likelihood depends on the fixed-effects design, REML likelihoods of models with different fixed effects cannot be compared. Models that differ in fixed effects are compared with ML likelihoods (or by Wald or F tests), while models with the same fixed effects but different random or covariance structures are compared with REML likelihoods. This rule is the most frequently broken in applied mixed-model work.
For a single random sample y_1, ..., y_n with unknown mean and variance, the ML estimate of the variance is sum of (y_i - y-bar) squared / n and the REML estimate is the same sum divided by n - 1. With n = 5 and a sum of squares of 40, ML gives 8 and REML gives 10. The difference shrinks as n grows, but in designs with many fixed effects relative to the number of observations (such as variety trials with many entries and few replicates) the downward bias of ML can be large, which is why REML is preferred. In a balanced one-way random model, REML gives sigma squared-hat = MSE and sigma_a squared-hat = max(0, (MSA - MSE)/n), the truncated ANOVA estimator, whereas ML divides the between-group sum of squares by a rather than a - 1.
An estimate of exactly zero for a variance component, a very large standard error for a component, a correlation estimated at plus or minus one, or a warning that the information matrix is singular all indicate that the model is over-parameterised for the data. The remedies are to simplify the random or covariance structure, to check the data for coding errors (for example a block labelled inconsistently across replicates), and to collect more levels of the random factor in future: variance components estimated from four or five levels are inherently imprecise.
REML can be pictured as a two-step process: first remove the fixed effects by working with residuals from an ordinary least squares fit (more exactly, with any full set of n - p linearly independent error contrasts), then apply maximum likelihood to those residuals. Because the residuals carry no information about beta, the variance estimates are not distorted by the estimation of the fixed effects. The REML likelihood is the same whichever set of error contrasts is chosen, which is why the method is well defined.
Given G and R, the generalised least squares estimator of beta and the predictor of u can be obtained together from Henderson's mixed-model equations (MME): [X'R^(-1)X, X'R^(-1)Z; Z'R^(-1)X, Z'R^(-1)Z + G^(-1)] [beta-hat; u-hat] = [X'R^(-1)y; Z'R^(-1)y]. When R = sigma squared I and G = sigma_u squared I, the equations reduce to the ordinary normal equations with lambda = sigma squared / sigma_u squared added to the diagonal of Z'Z. The solution beta-hat is the BLUE of beta, identical to GLS with V = ZGZ' + R, and u-hat is the BLUP of u, which can also be written u-hat = G Z' V^(-1)(y - X beta-hat). The MME avoid inverting V, which is n x n, and only require G^(-1), which made large animal-breeding evaluations feasible.
BLUPs are shrunk towards zero relative to the corresponding fixed-effect estimates. In a balanced one-way model, the BLUP of a group effect is [n sigma_a squared / (n sigma_a squared + sigma squared)] times the deviation of the group mean from the overall mean. With sigma_a squared = 8, sigma squared = 10 and n = 5, the shrinkage factor is 40/50 = 0.8, so a group whose mean is 5 units above the overall mean has a BLUP of 4. Groups with less data are shrunk more. Shrinkage is what makes BLUP the best predictor for selection: ranking lines by BLUPs maximises the expected merit of those selected, because lines that look good only through luck with few observations are pulled back. In animal breeding the same principle, with the relationship matrix A in G = A sigma_a squared, gives estimated breeding values that use information on all relatives. When the purpose is to compare a few fixed varieties precisely, however, BLUEs are reported, since shrinkage would bias the comparison.
Variance components are the building blocks of genetic parameters: heritability on a plot basis is sigma_g squared / (sigma_g squared + sigma_e squared); on a line-mean basis the error variance is divided by the number of replicates, which is why mean-basis heritabilities are higher. Repeatability, intraclass correlations and the proportion of variance due to blocks, sites or years are all ratios of components and are only as good as the estimates behind them, which with few levels of a random factor (four sites, three years) are imprecise. Standard errors of components, from the inverse of the information matrix, should be reported.
Suppose two groups have 5 and 2 observations, with group means 4 units and 4 units above the overall estimate, sigma squared = 10 and sigma_a squared = 8. The shrinkage factors are 5 x 8 / (5 x 8 + 10) = 0.8 and 2 x 8 / (2 x 8 + 10) = 0.615, so the BLUPs are 3.2 and 2.46: the group with less data is shrunk more, even though both observed deviations are equal. In the mixed-model equations this appears as the addition of lambda = 10/8 = 1.25 to the diagonal elements of Z'Z, which are the group sizes 5 and 2, so the denominators become 6.25 and 3.25 and the solutions are 5 x 4/6.25 = 3.2 and 2 x 4/3.25 = 2.46, the same values.
In breeding programmes, selection on BLUPs rather than raw means increases the expected genetic gain, particularly when lines differ in the amount of data (number of sites, years or relatives) available. The prediction error variance of each BLUP, obtained from the inverse of the coefficient matrix of the MME, gives the reliability of each estimated breeding value, which breeders use to decide how much further testing a candidate needs.
In the animal model, every animal in the pedigree has a random additive genetic effect, with G = A sigma_a squared, where A is the numerator relationship matrix computed from the pedigree. Henderson's rules for computing A^(-1) directly from the pedigree, without inverting A, made national genetic evaluations of dairy cattle possible from the 1970s onward. The same structure with a genomic relationship matrix computed from markers gives genomic BLUP, which is developed in the statistical genomics subject.
The sampling variance of an estimated variance component depends heavily on the number of levels of the random factor. With five sites, a site variance has only four degrees of freedom of information and its standard error is commonly as large as half the estimate or more; heritabilities from a single small trial are similarly imprecise. Pooling information over several trials, or using prior information in Bayesian analyses, gives more reliable parameters for planning breeding programmes.
In balanced designs, each fixed effect is tested against a single mean square with known degrees of freedom: the main-plot factor of a split plot against error (a), the sub-plot factor against error (b). With unbalanced data, missing plots or complex covariance structures, the variance of a fixed-effect contrast is a combination of several variance components, and the test statistic is only approximately t or F distributed with uncertain denominator degrees of freedom. Using the residual degrees of freedom overstates precision, sometimes badly.
Satterthwaite's approximation gives the degrees of freedom of a linear combination of mean squares: if L = sum of c_i MS_i, with MS_i on nu_i degrees of freedom, then nu is approximately L squared / sum of [(c_i MS_i) squared / nu_i]. In the split-plot comparison of two main-plot levels at the same sub-plot level, the variance involves (b - 1) Eb + Ea; with b = 3, Eb = 4 on 16 degrees of freedom and Ea = 12 on 6, nu = 20 squared / (8 squared/16 + 12 squared/6) = 400/28 = 14.3. In mixed-model software the same idea is applied to the estimated variance of any contrast. The Kenward-Roger method (1997) goes further: it adjusts the estimated variance-covariance matrix of the fixed effects for the uncertainty in the estimated variance parameters, which otherwise makes standard errors too small, and then uses a Satterthwaite-type degrees-of-freedom approximation for the adjusted F statistic. It performs well in small, unbalanced experiments and is the recommended default in many guides for agricultural trials. For balanced designs both methods reproduce the classical degrees of freedom.
For split-plot and other designs with several error strata, reports should state the degrees-of-freedom method, because it changes p-values and confidence intervals noticeably when the main-plot error has few degrees of freedom. When software offers residual degrees of freedom by default (as some generalised mixed-model programs do), the analyst should change the option for small experiments. Large-sample methods (Wald chi-square tests) are adequate only when all variance components are estimated from many levels.
Suppose a contrast's estimated variance is (1/4)(MS1 + MS2) with MS1 = 20 on 4 degrees of freedom and MS2 = 6 on 24. The Satterthwaite degrees of freedom are (20 + 6) squared / (20 squared/4 + 6 squared/24) = 676 / (100 + 1.5) = 6.7, much nearer the smaller of the two degrees of freedom because the mean square with fewer degrees of freedom dominates the sum. Using 24 or 28 degrees of freedom instead would give a t critical value too small and an interval that is too narrow.
Likelihood-based mixed-model analysis uses all available observations and is valid when data are missing at random, that is, when the probability of being missing depends only on observed quantities. A plot lost to a flood or an animal culled for a reason related to its earlier recorded weights satisfies this; an animal removed because of an unrecorded illness affecting the response may not. The classical missing-plot techniques of the basic designs subject are unnecessary in a mixed-model analysis, which handles imbalance directly.
Measurements taken repeatedly on the same units, such as weekly weights of animals, plant heights through a season or yields from successive harvests, are correlated within units, and the correlation usually declines as measurements get further apart in time. The classical split-plot-in-time analysis, described in the factorial and split-plot subject, assumes compound symmetry: all pairs of times equally correlated. Its F tests for time and time x treatment are exact only under sphericity (the Huynh-Feldt condition), and corrections such as Greenhouse-Geisser or Huynh-Feldt epsilon adjust the degrees of freedom when sphericity fails. The mixed-model approach instead models the covariance directly in R.
Compound symmetry (CS) has a common variance and a common covariance: 2 parameters. First-order autoregressive, AR(1), has variance sigma squared and correlation rho^|i - j| between times i and j: 2 parameters, suited to equally spaced times where correlation decays. Heterogeneous AR(1) lets each time have its own variance: t + 1 parameters for t times. Toeplitz allows a separate correlation for each lag: t parameters. Antedependence structures suit unequally spaced or developmental data. Unstructured (UN) estimates every variance and covariance: t(t + 1)/2 parameters, 10 for 4 times and 55 for 10 times. The choice balances fit against parsimony: unstructured is always valid but can be unstable with few subjects and many times; simple structures may be wrong. Structures with the same fixed effects are compared by REML likelihoods and information criteria (chapter 9). For unequally spaced times, spatial-type structures (power or exponential functions of the time difference) generalise AR(1). Random coefficient (random slope) models are an alternative that gives each unit its own growth line.
Consider body weights of 40 animals recorded at 6 monthly ages under four diets. Compound symmetry uses 2 covariance parameters, AR(1) 2, heterogeneous AR(1) 7, Toeplitz 6 and unstructured 21. If the REML log-likelihoods give AIC values of 1520 (CS), 1488 (AR(1)), 1480 (heterogeneous AR(1)) and 1484 (unstructured), heterogeneous AR(1) is preferred: variances change with age and correlation decays with time between weighings, and the unstructured model gains little for its extra parameters. The diet x age interaction would then be tested with the chosen structure and Kenward-Roger degrees of freedom.
An alternative to modelling the residual covariance is to give each experimental unit its own intercept and slope over time, with the intercepts and slopes as correlated random effects. This random coefficient (random regression) model implies a covariance among repeated measures that changes with time, and it describes individual growth or lactation curves directly. In animal breeding, random regression models on lactation stage are used to evaluate test-day milk records. The choice between random coefficients and residual structures, or a combination, is made by fit and by the scientific question.
When measurement times are irregular, for example plant heights at 14, 21, 35 and 63 days, structures defined on the time difference rather than on the order of measurement are needed. The spatial power structure, with correlation rho raised to the power of the time difference, reduces to AR(1) when times are equally spaced; exponential and Gaussian correlation functions are alternatives. Ignoring unequal spacing and fitting an ordinary AR(1) treats the 14-day and 28-day gaps as if they were equal.
Blocks capture fertility differences only if fertility changes at the scale and in the direction of the blocks. In practice field trends are smooth and irregular, and plots that are close together are more alike than plots in the same block but far apart. Spatial analysis models this directly, using the positions of plots in rows and columns, and typically reduces the error of variety comparisons substantially beyond what blocking achieves, especially in large breeding trials.
The earliest spatial method, nearest-neighbour analysis (Papadakis, 1937), adjusts each plot using the residuals of its neighbours: plot yields are first expressed as deviations from their variety means, the average deviation of the neighbouring plots is used as a covariate, and the analysis of covariance removes the local fertility component. Later work by Wilkinson, Besag and others refined the idea into first-difference and smoothing models. Nearest-neighbour methods are simple and effective but less flexible than modern mixed-model approaches.
Gilmour, Cullis and Verbyla (1997) proposed modelling plot errors as a separable first-order autoregressive process in rows and columns, AR1 x AR1, combined with an independent measurement error (the nugget) and with global trends and extraneous effects added where needed. Under AR1 x AR1 the correlation between two plots i rows apart and j columns apart is rho_r^i x rho_c^j; with rho_r = 0.6 and rho_c = 0.5, plots two rows and one column apart have correlation 0.36 x 0.5 = 0.18. The modelling process is iterative: fit AR1 x AR1, inspect the sample variogram of residuals (a three-dimensional plot of semivariance against row and column displacement), and add linear row or column trends, random row and column effects for operational patterns (such as harvest direction or irrigation lines) and other terms until the variogram looks stationary. Varieties are then compared with their BLUEs or BLUPs from the spatial model. This approach is standard in Australian and many other national variety-testing systems.
Field trials carry patterns created by operations as well as by soil: harvesting in serpentine order can create alternating column effects, irrigation runs create row gradients, and a sowing machine with several seed units can create periodic column effects. These are modelled as fixed or random row and column effects, linear trends or cyclic terms, alongside the AR1 x AR1 process for natural variation. Recording the field layout and operations is therefore part of the analysis, not only of the field book.
The sample variogram plots half the average squared difference between residuals against their separation in rows and columns. For independent residuals it is flat at the residual variance; for spatially correlated residuals it rises from a low value at small separations (the nugget) to a plateau (the sill). A variogram that keeps rising indicates a global trend not yet modelled; ridges along rows or columns indicate row or column effects; oscillations suggest periodic operational effects. Modelling stops when the variogram of the fitted model's residuals looks flat or matches the fitted correlation structure.
Published comparisons of spatial and non-spatial analyses of large variety trials commonly report clear reductions in the average standard error of variety differences, with the size of the gain depending on the field and the design; in trials with strong trends the gain can be large, in uniform fields small. Because the gain is trial-specific, programmes often fit both and choose by likelihood or information criteria, keeping design terms in all models.
An alternative to autoregressive correlation is to model the fertility surface as a smooth function of row and column positions. Two-dimensional penalised splines (P-splines) represent the surface as a sum of many local basis functions with a penalty on roughness; written as a mixed model, the penalty becomes a variance component estimated by REML, so the amount of smoothing is chosen by the data. The SpATS approach (Spatial Analysis of field Trials with Splines, Rodriguez-Alvarez and colleagues, 2018) uses an ANOVA-type decomposition of the two-dimensional spline (main smooth trends in rows and in columns plus their smooth interaction) and reports effective dimensions that show how complex the fitted surface is.
Spline models need less manual model building than AR1 x AR1 and are robust as a routine method for large numbers of trials, which is why many breeding programmes use them in automated pipelines; AR1 x AR1 with extraneous terms can capture local correlation that splines smooth over. Both are better than ignoring spatial variation, and both should be combined with the randomisation-based structure of the design (replicates and blocks). Diagnostics are the same: residual plots in field layout, variograms and the comparison of standard errors of variety differences with those of the non-spatial analysis.
In a P-spline model, the effective dimension of each smooth component summarises how many degrees of freedom the fitted surface uses: a value near zero means the component is almost flat, a large value means a complex surface. Reporting effective dimensions for the row trend, column trend and their interaction shows the user what kind of spatial pattern the trial had, and comparing across trials helps a programme understand its fields.
The two approaches are not exclusive. A model can include a smooth spline surface for large-scale trend and an AR1 x AR1 or independent residual for local variation, and some software fits both. The trade-off is complexity: more terms give more convergence problems and less transparent results, so the guiding principle is to add components only when diagnostics show they are needed.
Two nested models can be compared by the likelihood-ratio statistic, twice the difference in log-likelihoods (equivalently the difference in -2 log L), referred to a chi-square distribution with degrees of freedom equal to the difference in the number of parameters. For random effects and covariance parameters, REML likelihoods are used, with the same fixed effects in both models. Testing whether a single variance component is zero is a special case: the null value lies on the boundary of the parameter space, and the null distribution of the statistic is a 50:50 mixture of a point mass at zero and chi-square with 1 degree of freedom, so the correct p-value is half the naive chi-square(1) p-value. A statistic of 3.2, with naive p about 0.074, has a mixture p of about 0.037.
For non-nested structures, such as AR(1) against Toeplitz or AR1 x AR1 against a spline model, Akaike's information criterion AIC = -2 log L + 2k and the Bayesian criterion BIC = -2 log L + k log(n) are used, where k counts the variance parameters (with REML, models must share fixed effects; some definitions of n for REML BIC use n - p). Smaller is better. Information criteria are guides: two structures with similar AIC can give the same variety rankings, and the practical question is whether the choice changes conclusions. Overly complex covariance models may converge poorly or overfit, which a lower AIC alone does not reveal.
Suppose a model with random rows has -2 REML log-likelihood 812.4 and the model without them 815.6, both with the same fixed effects. The statistic is 3.2 on one variance parameter. The naive chi-square(1) p-value is about 0.074, which would suggest dropping the row effect at the 5% level; the correct 50:50 mixture p-value is about 0.037, which suggests keeping it. Because random terms that reflect the design (replicates, blocks) are usually kept regardless of significance, such tests are used mainly for optional terms like extraneous row and column effects.
Fixed effects are best tested with F tests using Kenward-Roger or Satterthwaite degrees of freedom in the final model, rather than with likelihood-ratio tests, which tend to be anti-conservative in small samples. If a likelihood-ratio test of fixed effects is used, both models must be fitted by ML, not REML, because REML likelihoods of models with different fixed effects are not comparable.
A useful working rule is to start from the design-based model (replicates, blocks, main plots) plus a plausible residual structure, then add or drop optional terms one at a time, checking diagnostics after each step. Changing many terms at once makes it hard to see which change mattered and often causes convergence failures.
In R, lme4 (function lmer) fits linear mixed models with crossed and nested random effects efficiently by REML or ML, and its companion glmer fits generalised linear mixed models; lme4 does not model residual correlation or heterogeneous residual variances directly. The lmerTest package adds Satterthwaite and Kenward-Roger degrees of freedom to lmer models (the latter through pbkrtest). The older nlme package (function lme) handles fewer crossed random effects but allows residual correlation structures (such as corAR1 or corSpatial) and variance functions (such as varIdent for different variances by group), which suits repeated measures and heterogeneous variances. Other packages include glmmTMB, sommer for genetic and spatial models, and SpATS for spline-based spatial analysis.
ASReml, developed by Gilmour, Cullis, Thompson and colleagues, is a commercial program (also available as an R package, asreml-R) built around the average-information REML algorithm. It handles very large animal and plant breeding models, complex G and R structures such as factor-analytic and AR1 x AR1 models, and pedigree or genomic relationship matrices, and is widely used in breeding programmes and variety testing. SAS PROC MIXED, with its RANDOM and REPEATED statements, fits a wide range of covariance structures and offers DDFM=SATTERTHWAITE and DDFM=KENWARDROGER options for denominator degrees of freedom; PROC GLIMMIX extends it to non-normal data. Different programs can give slightly different results for the same model because of parameterisation, constraints on variance parameters (allowing or forbidding negative estimates) and convergence criteria; reports should state the software, the estimation method and the degrees-of-freedom method.
Whatever the program, a mixed-model analysis should be scripted, with the data, model statements, software version and options recorded, so that the analysis can be repeated and audited. Variety-testing and breeding programmes that analyse thousands of trials a year depend on such pipelines, with automated checks for convergence, boundary estimates and outliers.
For crossed random effects with simple residuals, lme4 is fast and widely understood; for residual correlation and heterogeneous variances on modest data sets, nlme or glmmTMB; for large breeding data, factor-analytic and spatial models, ASReml or sommer; for organisations standardised on SAS, PROC MIXED and GLIMMIX. The statistical model should be decided first and the software chosen to fit it, not the reverse.
A complete report states: the model (fixed and random terms, residual structure), the estimation method (REML or ML), the software and version, how denominator degrees of freedom were computed, the estimated variance components with standard errors, how the model was chosen (likelihood-ratio tests, information criteria, diagnostics), and whether BLUEs or BLUPs are reported for treatments or genotypes. Readers need all of these to judge and reproduce the results.
Students need mixed models for nearly every modern thesis: unbalanced field trials, alpha and lattice designs with random blocks, repeated measurements on animals or plants, and heritability estimation. In examinations they must write expected mean squares and ANOVA estimators, explain REML, write the mixed-model equations, compute a BLUP shrinkage and know how many parameters a covariance structure has.
A five-year biometrician analyses dozens of trials a season with mixed models: random blocks and spatial terms for variety trials, REML variance components and BLUPs for breeding selections, repeated-measures structures for animal growth and lactation data. Daily decisions include which covariance structure to keep, whether a variance component is worth retaining, which degrees-of-freedom method to report and how to diagnose convergence problems.
A senior head of biometrics sets the standard analysis pipelines of a breeding or testing programme: spatial models for all trials, factor-analytic models across sites, BLUP-based selection, software and computing infrastructure, and reporting rules for variance components and standard errors. They decide when the programme adopts new methods such as spline-based spatial models or genomic relationship matrices.
A PhD with long experience develops mixed-model methodology (covariance models, spatial and genomic models, small-sample inference), teaches it, and advises programmes on analysis strategy. As a reviewer, they judge whether variance components, BLUPs and tests in a paper rest on an appropriate model.
Derivations of the REML likelihood and the mixed-model equations are outlined, not proved. Generalised linear mixed models, Bayesian mixed models and genomic prediction models are left to other subjects. Software details (function names, options) are given as at the time of writing and change between versions; check current documentation. The boundary correction for likelihood-ratio tests is given for a single variance component only.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions. Every numerical answer was computed with a script; the values needed are stated in the question. Software names and options are as at the time of writing.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.