Comparing ANOVA and Mixed Models: Genstat, ASReml and R

By Dr. Ruth Butler and Dr Vanessa Cave.

Mixed models are now widely used across agricultural, biological and breeding research, but researchers often face a practical question: do different software tools give the same answers, and where do the differences start to matter? In this technical article, Ruth Butler and Vanessa Cave compare Genstat, ASReml-R, lme4 and nlme using two worked examples, looking at model specification, predicted means, standard errors, degrees of freedom and covariance structures.

Why mixed models became central to modern analysisIn recent years, the use of mixed models has exploded, such that mixed-model analyses have become the norm in many fields, largely replacing traditional ANOVA. There are now many available tools for fitting mixed models. So, how do these compare?

Mixed models encompass a broad class of models that include both fixed and random effects. In this blog, we focus on linear mixed models (LMMs), which assume a continuous response and normally distributed random effects and residuals. 

The distinction between fixed and random effects is not always clear-cut: sometimes a variable might be treated as fixed in one analysis, and random in another. However, fixed effects are almost always variables with specifically chosen values (such as treatment factors like ‘fertilizer type’ or variables like ‘dose’); whereas random effects are generally things like ‘blocks’ or ‘plots’ which are randomly selected or arbitrary choices (Piepho et al. 2003 has some discussion on this).

Formally, the simplest linear mixed model has the form: 

y = Xβ + Zu + e

Here, X is the design matrix for the fixed effects β, and Z is the design matrix for the random effects u. The residual vector e is assumed to consist of independent and identically distributed (i.i.d.) normal random variables with mean 0 and variance σ². Likewise, in the simplest mixed model, the random effects u are assumed to be independent and normally distributed with variance σ2
u. This variance is called a ‘variance component’.

More complex models relax these assumptions and extend this basic framework. For example, random effects may be correlated, smoothing splines can be incorporated through additional random effects, and correlation structures can be used to model spatial relationships between adjacent plots or repeated observations on the same experimental unit.

Classical Analysis of Variance (ANOVA) can be viewed as a special case of linear mixed model analysis, albeit with different conceptual interpretations of the model parts and a different fitting algorithm. ANOVA in Genstat can be used provided there is general balance (Payne & Tobias 1992), but other implementations of ANOVA are more limited, fitting only designs with standard balance. Where the data are not balanced, alternative methods are required, and linear mixed models are good alternatives to ANOVA when the data are unbalanced. However, the ‘obvious’ approach of maximum likelihood gives estimates of the variance components that are biased downwards. Restricted (or residual) maximum likelihood (REML) was developed to reduce, and in many cases effectively remove, this bias. For guidance on which method is appropriate for your data, see this VSNi blog.

The foundations of mixed models can be traced back to the work of Ronald Fisher, with significant theoretical developments in the 1960s and 70s (e.g. Patterson & Thompson 1971), with specialised software for fitting mixed models becoming increasingly available during the 80s and 90s (e.g. Gilmour’s TwoD software (1992) for the spatial analysis of field-trial data). The methodology underlying TwoD subsequently formed the basis of ASReml, while similar REML facilities and command structures were incorporated into Genstat in 1990. In the later 90s, REML-based mixed-model procedures were incorporated into other commercial statistical software packages, notably SAS PROC MIXED (Littell et al., 1996). ASReml has been incorporated into R and Python packages (ASReml-R for some time; ASReml-Python in the last year). In parallel with these developments, the R commands lme (from the nlme package) and lmer (from lme4) have been produced: these use very different computational and optimisation approaches to ASReml. Another notable R package is “sommer”, which implements mixed-model fitting using alternative REML algorithms, particularly for quantitative genetics applications.

So, how do these various software packages compare? Are some easier to use? Do they produce similar results? Here we present analysis of two datasets using five packages/commands (Table 1) to help answer these questions. 

How model syntax differs between softwareThe symbols used to define the fixed and random parts of the model are similar, but not identical, between Genstat and R. Both use + and – to add or remove terms from a model. For an interaction term, Genstat uses the conventional dot (A.B) whilst R uses a colon (A:B). Both use * to denote main effect terms plus their interaction. So, A*B expands to A+B+A.B in Genstat and A+B+A:B in R. Both R and Genstat have / as shorthand for nesting: A/B means B nested within each level of A. So, A/B expands to A+A.B in Genstat and A+A:B in R (including in ASReml).

Fitting a mixed model typically involves commands for: 1) specifying the model (i.e., defining the fixed and random effects), 2) fitting the model, and 3) requesting additional output such as variance components, Wald tests, effects, or predicted means. Some of the main commands used for these tasks are shown in Table 1 for Genstat and three R packages.

Table 1: Setup/ main commands plus some further results commands, for each of the software commands/ R packages used. Fixed and Random effect commands: “fm” means fixed model and “rm” random model. “y” is the variable to analyse.

Variance-covariance and correlation models

Additional options to the main commands are available to allow for the inclusion of variance-covariance structures for the residual or other random terms. These structures, such as those used in Case Study 2: Slate Hall (see below), cannot be fitted using ANOVA, and are difficult to fit using lme4 because it provides only limited support for modelling residual variance–covariance structures. However, they can be fitted using Genstat’s REML commands, with a choice of about 18 variance-covariance models. ASReml provides a very comprehensive suite of variance and correlation structures, while the nlme package supports a smaller, but still useful, set of variance-covariance and correlation models.

Why deviances may look differentThe REML log-likelihood and the corresponding deviance are presented in various forms. (Note, the deviance is equal to -2 x the REML log-likelihood.) There are two additive constants that can be included or omitted from the reported deviance. The full deviance is then:

full deviance = deviance + (n-p)*log(2π) – log(det(X′X))

                     = deviance + pi – det

where pi = (n-p)*log(2π) and det = log(det(X′X)). 

Here, X is the design matrix for the fixed effects, n is the total number of usable units in the analysis, and p is the number of parameters fitted in the fixed model. Genstat presents the deviance with both the pi and det terms excluded. The deviance that includes the pi term can be obtained with the Genstat command VAIC with the include option set to ‘pi’:

VAIC [include=pi]

ASReml presents the log-likelihood, labelled ‘loglik’. This does not include the pi or det terms. Thus, the Genstat deviance is equal to the ASReml log-likelihood multiplied by -2. 

lme4 gives the ‘REML criterion at convergence’, which equals the Genstat deviance plus the pi term.

nlme gives the ‘Log-restricted-likelihood:’, which equals -(Genstat deviance + pi)/2.

(Note: The REML model can be fitted using HG commands. The likelihood statistic labelled 

‘-2 × P_β,v(h)’ equals the Genstat deviance + pi.)

Abbreviations used

Below, SEM, SED and df are, respectively, abbreviations for standard error of the mean, standard error of the difference between two means, and degrees of freedom.

Case study 1: Weed Competition

This trial is described under Examples 3.7, 9.5 and 16.1 of Welham et al. (2015). The data are available at https://emitanaka.r-universe.dev/smbdata.

The trial was conducted to assess the effect of irrigation and weed species on the yields of wheat. Treatments were arranged as a ‘split-plot design’, with two irrigation regimes (not irrigated; irrigated) and four weed species (blackgrass, cleavers, chickweed, plus no weeds). The irrigation regimes were applied to large ‘main plots’, with each main plot divided into 4 sub (or split) plots, each of which received one of the weed treatments. There were four replicates of the 8 treatment combinations, laid out in randomized complete blocks (Figure 1). 

The trial has a factorial treatment structure, with Weed × Irrigation fitted as fixed effects and Block, WholePlot, and Subplot representing the hierarchical randomization structure (Block/WholePlot/Subplot). In this dataset, the lowest-level term (Block.WholePlot.Subplot) uniquely identifies each row of the dataset and is therefore equivalent to the observational unit. Not all software recognizes that this term is equivalent to the observational unit, and some packages require it to be omitted from the random model specification for the analysis to fit correctly.

The analyses described below can be carried out in Genstat using the ANOVA and REML menus: 

• For ANOVA, Stats-Analysis of Variance, General…, select Split-plot design from the Design pull-down list. 
• For REML, Stats-Mixed Models (REML)-Linear Mixed Models…, enter the Fixed and Random formulae as immediately above.

Code for this example: The Genstat and R code used for the Weed Competition analysis is available to download here:

Note that below, ‘AOV’ is used as a label to identify the equivalent of a traditional analysis of variance table. However, these tables may summarise statistics other than mean squares and variance ratios, such as Wald tests. Results calculated ‘by hand’, as published in Welham et al. (2015) (‘Book’), are summarised in Table 2. 

Table 2: Case study 1: Summary of analyses results as in Welham et al. (2015).

The test results (p-values), the means, and the SEDs are the same regardless of the command used, Genstat or R. The SEM differs between the ANOVA and REML analyses. The Book and ANOVA give an SEM of 0.3613, whereas all REML methods give an SEM of 0.4063. 

The difference arises because the ANOVA SEM is an effective standard error intended for use in calculating standard errors of treatment comparisons and contrasts, whereas the REML SEM is the standard error of the predicted mean itself.

. (See Section 4.1.3 of the Guide to Genstat, Statistics.) As Welham et al. (2015) p 437, section 16.11 explains: “the difference occurs because the ANOVA SE ignores contributions from the block and wholeplot effects.” With respect to the SEDs, they go on to say they “are the same as the SED from ANOVA because the contributions from structural (random) terms cancel when taking differences”. 

For the R commands, the deviances all match the Genstat deviance once the adjustments for the additive constants, noted above, are made.

The primary differences in the results from the various commands are in the df associated with the SEM and SEDs (Table 3).

Table 3: Summary of degrees of freedom obtained from the various estimation methods, for SEM and SED, for the full Interaction by Species table of means.

Degrees of freedom associated with the SEM and SED vary. For ANOVA and the Book, the df for SEM and the SED for a comparison between means with differing irrigation treatments are given as 6.64, obtained using Satterthwaite’s approximation. ANOVA also has 18: 6.64 is the more general one. The largest d.f (18) are for comparisons between means that have the same main plot (irrigation) treatment, because the comparison does not involve the ‘main plot’ estimated variance. The smaller df is to compare means that have different irrigation treatments, because the main plot variance contributes to the calculation of the SED. For Genstat REML, the denominator degrees of freedom from the Wald tests are used as a default for the SEM (giving 18 df). For R, both predictmeans and emmeans following a lme4 fit give 11.1 df.When followed by predictmeans, nlme gives 18, whilst emmeans 3. 

For the SED, the df in the Book and with ANOVA are as in Table 2. For REML, the Genstat command VLSD gives a range of 3 to 18. For R, after fitting with lme4, predictmeans gives18 d.f and emmeans gives two values (6.64, 18), matching those for ANOVA. After fitting with nlme, predictmeans has 18 and emmeans 3 and 18.

What the Weed Competition example showsResults other than df are the same for all of the methods. There is no universally agreed method for calculating the degrees of freedom associated with the SED and SEM. For the SEM, values of 6.64 or 11.1 seem sensible, but it is not clear how 11.1 was obtained. Methods mentioned in the help for these commands are Kenwood-Roger and Satterthwaite’s method. The differing degrees of freedom reported by Genstat, ASReml, emmeans and predictmeans reflect differences in the methods used to approximate denominator degrees of freedom.

Case Study 2: Slate Hall

This data set is distributed with the ASReml package as dataset ‘shf’ and also with Genstat in the example file Slatehall.gsh (Help-Example Programs- Data Sets…). It appears in the ASREML manual (4.2-4; page 121) and also in the Genstat Guide to REML (page 38 for Release 24.2) and the Guide to statistics (Example 5.4.4a, page 737). The column names vary between the Genstat and ASReml versions of the data. Note also that the Genstat REML Guide version has the yields divided by 100. Here, we present the results using the ASReml column names, plus yields as (g/m²), as in the original Gilmour paper. The data set is also in the R agridat package as kempton.slatehall, but the factors RowBlk and ColBlk are not included.

The trial comprised six replicates for each of 25 varieties of wheat. Each replicate was laid out as a 5×5 array, and the replicates were arranged in two rows by three columns (Figure 2). 

The trial was designed as a balanced lattice square with six replicates, with plots laid out in a 10 × 15 rectangular grid. Unlike Case Study 1, this trial is not orthogonal, but it does possess general balance. It can be analysed using Genstat’s ANOVA, because of the powerful algorithm that Genstat uses.

The Genstat and ASReml documents describe two possible analyses. The first model is a separable AR1(row) x AR1(column) model. That is, correlations are estimated in both directions. lme4 does not have the capability to fit the required correlation model, and extensive web searching and testing of suggested code failed to identify a way to fit the same model using nlme. (With our prompt, code suggested by ChatGPT did not work!). The R package ‘sommer’ can fit the correlation model, but that is not explored here.

The second analysis fits the original design: a balanced lattice square design (BLS). This can be fitted with both lme4 and nlme and results matched those for REML and ASReml.

The spatial and BLS models can be fitted in Genstat using the menus, as described in the REML guide.

Code for this example: The Genstat and R code used for the Slate Hall analysis is available to download here:

In the Genstat code linked above, variable names have been changed from those used in the spreadsheet file to match those used in the shf data frame in the asreml package. Rows and Columns relate to the entire field, while RowBlk and ColBlk number rows and columns within replicates (Rep). In the ASReml code, the “v” as in ar1v specifies a homogeneous variance structure, giving a common variance for all observations.

What the Slate Hall example shows, ss for Case Study 1, the means, SEMs and SEDs match between Genstat REML and ASReml. In addition, the df associated with the SEMs and SEDs were the same for ASReml and Genstat’s REML, with 79 for the BLS model and 80 for the spatial model.

What this means in practice

For the examples considered here, the predicted means, standard errors and test results are the same for each fitting method, with the primary differences occurring in the df associated with the SEMs and SEDs. This seems to be because there is no consensus as to which calculation method to use. However, if, as in Case Study 1 the range of possible degrees of freedom is moderately large, the choice of df can mean that confidence limits or LSDs vary noticeably. The calculations for these statistics use a t-value, which for 3 df is about 3.2 and for 18 df, 2.1. So, confidence limits and LSDs will be about 50% larger with 3 df than with 18: this is a sizeable difference. Any reader using error bars to compare means will find that the results of tests between means will vary substantially depending on the choice of df (One can argue that a focus on “significant” versus “not significant” is not always meaningful, but this is often how results are interpreted in practice.)

Since all of the methods considered here give the same results for the major pieces of information of interest (e.g., means, standard errors, deviances), the choice of method is largely a matter of preference. Genstat’s REML menus are perhaps the most accessible way to fit complex mixed models. However, they do not cover the entire range of models that Genstat can fit, so the use of code may sometimes be required. ASReml-R supports a particularly comprehensive range of variance and correlation structures, many of which are difficult to fit in other software. There are other packages such as SAS and sommer which do allow a wide range of covariance models to be fitted. lme4 is limited in that it can only fit a very small range of covariance structures through random effects (such as random coefficients). nlme can fit a wider range of mainly 1 dimensional correlation models but is still more limited in comparison with Genstat and especially ASReml.

The formulation of the model, and especially the specification of random effects, is perhaps easier in Genstat and ASReml than in lme4 or nlme. The latter use (the somewhat mysterious) (1| rm) to define random effects, which may initially seem unfamiliar to new users. Models with crossed random effects can be very tricky to specify for nlme or lme4. In addition, some correlation structures require spatial coordinates to be supplied as numeric variables rather than factors.

Finally, the speed of fitting models, especially complex models, varies. In our experience, ANOVA is the fastest for generally-balanced designs. ASReml is typically faster than the other software considered here when fitting large and complex mixed models, particularly for the types of models commonly encountered in plant and animal breeding and agricultural research. 

References

Gilmour, A. 1992. TwoD: a program to fit a mixed linear model with two dimensional spatial adjustment for local trend [computer software]. NSW Agriculture Biometric Bulletin (Australia). no. 1.

Patterson, H.D. & Thompson, R. 1971. Recovery of Inter-Block Information when Block Sizes are Unequal. Biometrika 58(3), 545–554.

Payne, R.W. & Tobias, R.D. 1992. General balance, combination of information and the analysis of covariance. Scandinavian journal of statistics, 3–23.

Piepho, H.P., Büchse, A. & Emrich, K. 2003. A Hitchhiker’s Guide to Mixed Models for Randomized Experiments. Journal of Agronomy and Crop Science 189, 310–322.

Welham, S.J., Gezan, S.A., Clark, S.J. & Mead, A. 2015. Statistical methods in biology: design and analysis of experiments and regression. CRC Press, Boca Raton, Florida,Pp 563+xiii.