Our analytics software and consulting services are chosen by seed, plant, aqua and animal breeding companies worldwide to support and inform the development of new varieties, strains, stocks and breeds.

Case Studies

Our software in action

Using mainly the GLM functionality within Genstat, large datasets are quickly and easily analysed with the results factored back into future experimental design.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/elephants_genstat_4e5d18efcb.png)

The conservation of the African elephant provides tremendous opportunities for simultaneously conserving biodiversity and increasing benefits to local communities.

Their role as a flagship species of global significance helps maintain biodiversity of the ecosystems they inhabit. In addition, the cultural and aesthetic values of elephants are also important, not only to African societies, but to the world at large. But this needs to be balanced against the fact that as the African elephant range becomes more and more fragmented and elephants get confined into smaller pockets of suitable habitat, humans and elephants are increasingly coming into contact and in conflict with each other. Tracking and recording the actual numbers of elephants in Africa is difficult; data are collated using a wide variety of techniques from ground or aerial surveys to dung counts and informed and other ‘guesstimation’. Overlay the obvious difficulties in the range and distances elephants can travel (global satellite tracking is helping this process) means putting accurate numbers on elephant populations can be, excuse the pun, a mammoth task. That said, counts are achieved and have shown the elephant population in Kenya has gone through various changes. In 1973, the population was about 167,000 declining to just 19,000 in 19891. The African Elephant Status Report 20072 (AESR 2007), recently published by the African Elephant Specialist Group (AfESG), puts the definite elephant count at 23,353, up 1,317 on the 2002 survey results. Add in probable, possible and speculative counts and the 2006 figure is 31,636; a rise of 2,830 on the 2002 figures. Alongside the elephant population growth, the human population in Africa has escalated way beyond anyone’s predictions, reaching an estimated 924 million3 in 2006. This increasing human population is starting to encroach dramatically into traditionally wildlife rich areas forcing large animals like elephants (and lions, buffalos and leopards incidentally) into smaller and smaller areas. Of course, 6 tonne animals like elephants have a huge requirement for food and water and naturally migrate across the landscape searching for sustenance. This is bringing them into direct conflict with people as new villages, roads, schools, fences, bridges and farms are being built over these natural wildlife corridors. The human-elephant conflict is becoming a serious issue in Kenya as people are now killing, spearing or poisoning elephants that come onto their land to eat their crops. So while elephants are major economic assets in generating tourist revenue, crop raiding by elephants is a cause of increasing levels of conflict between local people and wildlife management strategies across Africa. The economic damage to small-scale farmers can be crippling. Unless a cost-effective method of limiting crop damage is found, the pressure on wildlife managers to implement radical methods such as culling, will inevitably become a reality. As Lucy herself says, “When you meet people who are affected by these invading elephants you can only feel the deepest sympathy that they have lost their entire harvest. This causes genuine stress, damage and sometimes even death to these poor rural villagers.”

#### **What can be done?**

One of the methods that has been tried with variable success is the installation of electric fences to keep wild animals, particularly elephants, out of people’s farms. However, this is very expensive, requires a reliable supply of electricity and cannot be financed by the farmers so they become reliant on donors. Lucy works with a research team called Save the Elephants, which is headed by the elephant expert Dr Iain Douglas-Hamilton OBE; an organisation dedicated to developing a tolerant relationship between elephants and humans. Her research is furthering the aim to find crop-protection methods that can be financed and managed by the farmers themselves. The use of chilli peppers, either as disincentive crop or as a deterrent spray has been pioneered by Osborn (2002) in Zimbabwe with promising, albeit limited, success. Other variations include light-weight fencing hung with cow bells and laced with a potent mix of red chilli and black engine oil, thought to both alarm and irritate the elephants’ huge, sensitive trunks. However, Lucy’s research, based upon Kenya folklore which says that elephants are afraid of bees, would be an innovative and socially appropriate solution.

#### **David and Goliath**

The project explores the use of traditional, wooden beehives as both an elephant deterrent and a social and economic boost to poverty-stricken rural communities through the sustainable harvesting of honey. Her earlier research has shown that elephants react with obvious alarm and retreating behaviour to bee sounds. Further, recent studies show that elephants will actually avoid trees with beehives. By setting hives in acacia trees and recording how elephants react to them, Lucy hopes to be able to demonstrate that beehives can be used as an eco-barrier protecting fields of crops. When she isn’t persuading the British Army to help her fix her overworked Land Rover, Lucy works from a green canvas tent alongside the beautiful Ewaso Ng’iro river, in the Samburu Reserve.Day to day work involves setting up sample plots and selecting trees for beehives to be hung in. Many plots are tested which means driving through acres of park (hence the overworked Landie) and frequent encounters with elephants. Lucy has employed a local assistant Lukas, from the Ndorobo tribe; well known for their beekeeping abilities, and Lukas has 13 hives of his own which he helps to support his family with. There are also 4 long term elephant monitors, 1 MSc student from Nairobi University studying debarking effects of elephants and 1 other Kenyan researching whether or not lions use the same wildlife corridors as elephants. In addition, a team of researchers, all local Samburu are trained to identify all the elephants in the reserve; incredibly they know about 750 elephants by sight alone

#### **Data Analysis**

As well as the bee data, data are also collected from observations into factors such as temperature, wind direction and time of day to establish if the elephant behaviours observed are due to these factors, or if the elephants are indeed responding to the presence of the bees.

Using mainly the GLM functionality within Genstat, large datasets are quickly and easily analysed with the results factored back into future experimental design.

#### **Now, next steps and the future**

In this day and age with all the issues of environment and climate change to be able to find a way of assisting farmers without destroying the natural habitat must be of the utmost importance. Here we have an opportunity to be a part of such a programme; and additionally to prove once again to the non-statistician that data analysis and statistics are a crucial part of everyday living.

[http://elephantsandbees.com/](http://elephantsandbees.com/)

large_Bees_mainimage_desktop.jpg

medium_Bees_mainimage_desktop.jpg

small_Bees_mainimage_desktop.jpg

thumbnail_Bees_mainimage_desktop.jpg

Bees & Elephants

Most breeding programs plan several controlled crosses between outstanding parents to detect favorable alleles in their offspring. The progeny are later evaluated in a field experiment, and this information is used to assess the genetic worth of the parents by fitting parental linear mixed models (LMMs) and obtaining best linear unbiased predictions (BLUPs). These BLUPs, which are the _general combining ability_ (GCA), or 1/2 of the _breeding value_ (BV, with BV = 2 ×× GCA) of each parent, are then used to select the best parents for future crosses or operational deployment.

In many plant breeding programs, a parent is considered in several crosses. In most cases it is easy to assign the sex of a given individual. However, several commercial plant species are monoecious, which means that a given genotype will bear both male and female flowers. In contrast, dioecious species have distinctive male and female plants. Some examples of monoecious species are corn, squash, banana, and many conifers, particularly those of the genus _Pinus_.

In quantitative genetic analyses, monoecious species present a particular challenge, as a given parent can contribute to the estimation of its breeding value (or GCA) as both a male and a female, something that needs to be taken into consideration when a statistical model is fitted.

This is done by _overlaying_ design matrices of the factors associated with male and female parents. We will illustrate this here using an example of loblolly pine (_Pinus taeda_) in which some individuals, depending on the availability of pollen or flowers, were used in several artificial crosses as both male and female in the breeding program.

The data used here originates from a loblolly pine clonal study published by Resende _et al._ (2012). Parents were crossed in a circular mating design, constituting several full-sib families. Individuals from these families were vegetatively propagated (cloned) and established in a series of field trials. A subset of the full dataset, corresponding to diameter at breast height (DBH, inches) measured at 6 years since planting at the Nassau (Florida, USA) site, is used here. We will be using the adjusted mean values for this trait. In this dataset we have a total of 71 families that originated from 43 parents; however, 20 of those parents were used as both males and females in different crosses. There were no reciprocal or self-pollinated crosses planned, but these can occur in other crops and they might need to be identified and modelled properly. In this dataset, each family is formed by between 1 and 16 individuals (with an average of 12.2). We also have pedigree information for the 43 parents.

The first few lines of the phenotypic dataset $\\texttt{PEDPAR\_NASSAU.TXT}$ are presented below:

| tree | mother | father | family | DBH |
| --- | --- | --- | --- | --- |
| 280200 | 20012 | 44112 | 1017 | 9.410 |
| 282200 | 20022 | 44090 | 1020 | 11.953 |
| 284400 | 44126 | 142170 | 1049 | 10.785 |
| 285210 | 50152 | 58040 | 1056 | 9.950 |
| 285600 | 142024 | 44112 | 1006 | 10.325 |
| 286600 | 44012 | 17766 | 1036 | 11.122 |

And the first and last few lines of parental pedigree file $\\texttt{PEDPAR\_NASSAU.TXT}$ are:

|   | genotype | mother | father |
| --- | --- | --- | --- |
| 1 | 14006 | 0 | 0 |
| 2 | 14046 | 0 | 0 |
| 3 | 14060 | 0 | 0 |
| 4 | 14070 | 0 | 0 |
| 5 | 14104 | 0 | 0 |
| 6 | 14106 | 0 | 0 |
| 48 | 142076 | 14114 | 14104 |
| 49 | 142170 | 14106 | 14046 |
| 50 | 142194 | 14216 | 14060 |
| 51 | 202060 | 20030 | 20052 |
| 52 | 202096 | 20034 | 20052 |
| 53 | 222056 | 22034 | 14006 |
| 54 | 222060 | 22020 | 22034 |

In the phenotypic data each row contains the id of the individual tree, together with information about their mother and father, followed by the adjusted phenotypic response DBH. In the pedigree file there is information about the mother and father of each of the parents used in the crosses. For example, note that tree $\\texttt{284400}$ has father $\\texttt{142170}$, and this genotype has as parents genotypes $\\texttt{14106}$ and $\\texttt{14046}$.

We are interested in fitting a parental model, where our objective is to estimate GCA values for each of the parents. (Note that this is different than an animal model, where we estimate a BV for each of the individuals (parents or offspring).) Our parental base model is:

$$

y\_{ijk} = \\mu + f\_i + m\_j + fm\_{ij} + \\epsilon\_{ijk}

$$

where, $y\_{ijk}$ is the response for the kth individual originating from the cross between the ith female with the jth male; $\\mu$ is the overall mean; $f\_i$ is the random effect of the ith female, with $f\_i \\sim N(0, \\sigma^2\_f)$; $m\_j$ is the random effect of the jth male, with $m\_j \\sim N(0, \\sigma^2\_m)$; and $fm\_{ij}$ is the interaction or family effect of the cross between the ith female with the jth male, with $fm\_{ij} \\sim N(0, \\sigma^2\_{fm})$; and $\\epsilon\_{ijk}$ is the random residual with $\\epsilon\_{ijk} \\sim N(0, \\sigma^2)$.

Note that for the parental effects $f\_i$ and $m\_j$ there is a pedigree available, which is used to obtain the numerator relationship matrix $\\mathbf{A}$, that is incorporated in our analyses. Hence, the vectors of parental effects have the assumptions $\\mathbf{f} \\sim N(\\mathbf{0},\\sigma^2\_f \\mathbf{A})$ and $\\mathbf{m} \\sim N(\\mathbf{0},\\sigma^2\_m \\mathbf{A})$.

In the above model, we are assuming that there are distinct variances for females and males, $\\sigma^2\_f$ and $\\sigma^2\_m$, and that these can take any positive value. 

We can fit the above model with ASReml-R using the following code:

```plaintext
ainv &lt;- ainverse(PED)
modP &lt;- asreml(fixed=DBH~1,
             random=~vm(mother,ainv)+vm(father,ainv)+family,
             residual=~idv(units),
             data=PINE)
```

In this code we are generating the numerator relationship matrix $\\mathbf{A}$ using the command $\\texttt{ainverse()}$. The matrix $\\mathbf{A}$ is incorporated into the model using the function $\\texttt{vm()}$. The log-likelihood of the above model fitted to our data is $\\texttt{-744.488}$ and the estimated variance components are:

|   | component | std.error | z.ratio bound | $ch |
| --- | --- | --- | --- | --- |
| vm(mother, ainv) | 0.2774623 | 0.1121623 | 2.473759 P | 0 |
| vm(father, ainv) | 0.3815927 | 0.1348447 | 2.829868 P | 0 |
| family | 0.0000004 | NA | NA B | 0 |
| units!units | 1.8702273 | 0.0930707 | 20.094696 P | 0 |
| units!R | 1.0000000 | NA | NA F | 0 |

Note the difference in variance components between males (fathers, $\\texttt{0.277}$) and females (mothers, $\\texttt{0.382}$). Finding a difference in variance components by sex is common for plant species and it might be correct depending on the origin of the parents, but also other biological reasons might explain these differences. Also note that the family variance is bounded at almost zero.

In our case, we can calculate both the narrow-sense heritability for females and males ($h^2\_f$ and $h^2\_m$) and the dominance ratio ($d^2$), by using the following expressions:

$h^2\_f = 4 \\times \\sigma^2\_f / (\\sigma^2\_f + \\sigma^2\_m + \\sigma^2\_{fm} + \\sigma^2)$

$h^2\_m = 4 \\times \\sigma^2\_m / (\\sigma^2\_f + \\sigma^2\_m + \\sigma^2\_{fm} + \\sigma^2)$

$d^2 = 4 \\times \\sigma^2\_{fm} / (\\sigma^2\_f + \\sigma^2\_m + \\sigma^2\_{fm} + \\sigma^2)$

In some cases, particualrly when the estimations of heritability for both sexes might not be very different, these can be combined into a single estimate, as:

$$

h^2 = 2 \\times (\\sigma^2\_f +\\sigma^2\_m) / (\\sigma^2\_f + \\sigma^2\_m + \\sigma^2\_{fm} + \\sigma^2)

$$

All of the above expressions were obtained in ASReml-R using the command $\\texttt{vpredict()}$ as shown below:

```plaintext
vpredict(modP,h2m~4*V1/(V1+V2+V3+V4))
vpredict(modP,h2f~4*V2/(V1+V2+V3+V4))
vpredict(modP,d2~4*V3/(V1+V2+V3+V4))
vpredict(modP,h2~2*(V1+V2)/(V1+V2+V3+V4))
```

with the results:

|   | Estimate | SE |
| --- | --- | --- |
| h2m | 0.4388000 | 0.1616208 |
| h2f | 0.6034797 | 0.1858860 |
| d2 | 0.0000007 | 0.0000001 |
| h2 | 0.5211399 | 0.1043225 |

As expected, our dominance ratio $d^2$ is very small in this particular trait; hence, it could be ignored. Also, the narrow-sense heritabilities are reasonable for this type of study, with a combined estimate of $\\texttt{0.521}$ but with a moderate approximated standard error ($\\texttt{0.104}$).

Also, we can obtain the BLUPs for the model terms $\\texttt{mother}$ and $\\texttt{father}$, which in this case will correspond to the GCA values. If these are multiplied by two, we will have the BVs. The GCAs are obtained with the code:

```plaintext
summary(modP,coef=TRUE)$coef.random
```

In the output below, we can observe a few of the estimated GCA values for the $\\texttt{mother}$ (top) and $\\texttt{father}$ (bottom).

|   | solution | std.error | z.ratio |
| --- | --- | --- | --- |
| vm(mother, ainv)\_14006 | 0.3600002 | 0.2841548 | 1.2669160 |
| vm(mother, ainv)\_14046 | 0.0000000 | 0.5267472 | 0.0000000 |
| vm(mother, ainv)\_14060 | \-0.2286267 | 0.4864881 | \-0.4699533 |
| vm(mother, ainv)\_14070 | 0.0000000 | 0.5267472 | 0.0000000 |
| vm(mother, ainv)\_14104 | \-0.4620406 | 0.4208766 | \-1.0978056 |
| vm(mother, ainv)\_14106 | 0.0000000 | 0.5267472 | 0.0000000 |

|   | solution | std.error | z.ratio |
| --- | --- | --- | --- |
| vm(father, ainv)\_14006 | 0.0013612 | 0.2827143 | 0.0048148 |
| vm(father, ainv)\_14046 | \-0.5495441 | 0.5510751 | \-0.9972219 |
| vm(father, ainv)\_14060 | 0.0000000 | 0.6177321 | 0.0000000 |
| vm(father, ainv)\_14070 | \-0.3792472 | 0.3257064 | \-1.1643837 |
| vm(father, ainv)\_14104 | \-0.2087355 | 0.5024186 | \-0.4154613 |
| vm(father, ainv)\_14106 | \-0.5495441 | 0.5510751 | \-0.9972219 |
| vm(father, ainv)\_14114 | 0.1476378 | 0.3462045 | 0.4264466 |

Note that for parent $\\texttt{14006}$ the estimates are different as a mother and a father ($\\texttt{0.360}$ against $\\texttt{0.001}$), and their precision will depend on the amount of information a parent has (_i.e._, number of offspring), but also in its relationship with other genotypes on the dataset.

So far, we have not combined the information from a given parent when it is used as a male _and_ as a female. Now, it is of interest to _overlay_ the factors $\\texttt{mother}$ and $\\texttt{father}$. To do this in ASReml-R we require two elements: the use of the statement $\\texttt{and()}$ in the model and the specification of the command $\\texttt{equate.levels}$. This is illustrated in the following code for the same data as above:

```plaintext
modO <- asreml(fixed=DBH~1,
               random=~vm(mother,ainv)+and(vm(father,ainv)),
               equate.levels=c('mother','father'),
               residual=~idv(units),
               data=PINE)
```

First, let’s focus on the $\\texttt{and()}$ function. This is used to request that the design matrices associated with $\\texttt{father}$ are _overlaid_ on the design matrix associated with the immediately previous term (here, $\\texttt{mother}$). Second, the expression $\\texttt{equate.levels}$ requests that the levels for the factor $\\texttt{mother}$ and $\\texttt{father}$ are treated as identical. That is, levels that have the same name in the factors $\\texttt{mother}$ and $\\texttt{father}$ will be treated as the same, and levels that are found in only one of the factors will be added to the new design matrix. Both of the above commands are required in ASReml-R to ensure that consistent BV (or GCA) estimates are obtained.

In order to illustrate the above overlay of matrices, we will use a simple example, in which we have eight individuals originating from crosses of four parents (P1, P2, P3 and P4). These parents were used as female or male, and in some cases we have reciprocal crosses and selfing. The crossing and design matrices for our base model (without overlay) are shown below:

$$

\\begin{matrix}  

indiv & female & male \\\\  

A & P1 & P2 \\\\  

B & P1 & P3 \\\\  

C & P1 & P4 \\\\  

D & P2 & P3 \\\\  

E & P2 & P4 \\\\  

F & P3 & P2 \\\\  

G & P3 & P4 \\\\  

H & P4 & P4 \\\\  

\\end{matrix}  

\\Rightarrow  

\\begin{matrix}  

P1 \\:\\: P2 \\:\\: P3 \\:\\: P4 \\\\  

\\begin{bmatrix}  

1 & 0 & 0 & 0 \\\\  

1 & 0 & 0 & 0 \\\\  

1 & 0 & 0 & 0 \\\\  

0 & 1 & 0 & 0 \\\\  

0 & 1 & 0 & 0 \\\\  

0 & 0 & 1 & 0 \\\\  

0 & 0 & 1 & 0 \\\\  

0 & 0 & 0 & 1 \\\\  

\\end{bmatrix}  

\\end{matrix}  

\\:\\:  

\\begin{matrix}  

P2 \\:\\: P3 \\:\\: P4 \\\\  

\\begin{bmatrix}  

1 & 0 & 0 \\\\  

0 & 1 & 0 \\\\  

0 & 0 & 1 \\\\  

0 & 1 & 0 \\\\  

0 & 0 & 1 \\\\  

1 & 0 & 0 \\\\  

0 & 0 & 1 \\\\  

0 & 0 & 1 \\\\  

\\end{bmatrix}  

\\end{matrix}

$$

Once the LMM model is fitted based on the above matrices (as we did with our previous object $\\texttt{modP}$), we will have a total of four GCAs for females and three for males. If we count the number of ones for each effect, we find that they have between one and three pieces of information to be used in their estimation.

Now, if we proceed to overlay the _**female**_ and _**male**_ incidence matrices, we obtain the following results:

$$

\\begin{matrix}  

indiv & female & male \\\\  

A & P1 & P2 \\\\  

B & P1 & P3 \\\\  

C & P1 & P4 \\\\  

D & P2 & P3 \\\\  

E & P2 & P4 \\\\  

F & P3 & P2 \\\\  

G & P3 & P4 \\\\  

H & P4 & P4 \\\\  

\\end{matrix}  

\\Rightarrow  

\\begin{matrix}  

P1 \\:\\: P2 \\:\\: P3 \\:\\: P4 \\\\  

\\begin{bmatrix}  

1 & 1 & 0 & 0 \\\\  

1 & 0 & 1 & 0 \\\\  

1 & 0 & 0 & 1 \\\\  

0 & 1 & 1 & 0 \\\\  

0 & 1 & 0 & 1 \\\\  

0 & 1 & 1 & 0 \\\\  

0 & 0 & 1 & 1 \\\\  

0 & 0 & 0 & 2 \\\\  

\\end{bmatrix}  

\\end{matrix}

$$

You can observe that there is a single design matrix but only with four columns: one for each parent. In this case, we have more pieces of information to estimate each effect, ranging from three to six. Of special interest is the effect of parent P3, it now uses four pieces of information to estimate its single GCA, whereas before it was being estimated as two effects, each with two pieces of information.

There are other interesting things in the above matrix, such as the value 2 associated with individual H, that originated from the cross P4 $\\times$ P4, which corresponds to a selfing of parent P4. Also, individuals D and F are full-sibs, originating from reciprocal crosses, but in this case there is no distinction between their parental effects. Both, selfing and reciprocal crosses can be considered in the fitted linear mixed model as extensions with the aim of evaluating their importance or significance.

Finally, we can observe our results from fitting the above overlay model to our loblolly pine data. The estimated variance components are:

|   | component | std.error | z.ratio bound | %ch |   |
| --- | --- | --- | --- | --- | --- |
| vm(mother, ainv) | 0.4000453 | 0.1162417 | 3.441494 | P | 0 |
| units!units | 1.8599439 | 0.0917484 | 20.272218 | P | 0 |
| units!R | 1.0000000 | NA | NA | F | 0 |

In this particular case, we are missing the variance associated with the factor $\\texttt{father}$. However, this is still part of the model. As the matrices were overlaid, we have as an output a single parental effect now called $\\texttt{mother}$, which represents both male and female parents. Also, note that we dropped the term associated with $\\texttt{family}$, as in the previous analysis this was bound to zero. The log-likelihood of this model is $\\texttt{-737.328}$, a greater value indicating a better fit than the obtained from the original model ($\\texttt{-744.488}$), and probably a significant improvement.

We can now proceed to obtain our estimate of heritability using the expression:

```plaintext
vpredict(modO,h2~4*V1/(2*V1+V2))
```

which is the traditional formula; however, note that in the denominator we are considering twice the $\\texttt{mother}$ variance as we have a single parental variance component. The results from the above code are presented below:

|   | Estimate | SE |
| --- | --- | --- |
| h2 | 0.601564 | 0.1248922 |

As expected, these results are different to those from our previously fitted model because of the way this model is combining the available information. There is an increase in the heritability when we overlay the factors of mother and father compared to when we do not ($h^2 = 0.521$). This increase is expected as we now have more information to estimate each of the effects. Also, we expect a reduction in the standard errors of the estimated GCA values, which are presented below:

|   | solution | std.error | z.ratio |
| --- | --- | --- | --- |
| vm(mother, ainv)\_14006 | 0.5174180 | 0.2116542 | 2.4446388 |
| vm(mother, ainv)\_14046 | \-0.5055163 | 0.5609615 | \-0.9011606 |
| vm(mother, ainv)\_14060 | \-0.2843599 | 0.5668751 | \-0.5016271 |
| vm(mother, ainv)\_14070 | \-0.3430721 | 0.3366609 | \-1.0190434 |
| vm(mother, ainv)\_14104 | \-0.4849910 | 0.4781210 | \-1.0143688 |
| vm(mother, ainv)\_14106 | \-0.5055163 | 0.5609615 | \-0.9011606 |

In this case we have only $\\texttt{mother}$ effects that identify each of the levels of the parent. In contrast with our previous model, parent $\\texttt{14006}$ has a single GCA value ($\\texttt{0.517}$), and its standard error is smaller than before ($\\texttt{0.21}$ compared to $\\texttt{0.28}$). Some GCA estimates will increase or decrease, and their precision might get better or worse depending on how information is combined.

In this example, we have explored a parental model where the effects of male and female were overlaid. This resulted, according to the log-likelihood value, in a statistically superior model that combines the data better. But more importantly, this model is more biologically accurate given our knowledge of monoecy of the species under study.

### **Author**

Salvador A. Gezan

### **References**

Resende, MFR, Munoz, P, Resende, MDV, Garrick, DJ, Fernando, RL, Davis, JM, Jokela, EJ, Martin, TA, Peter, GF and M. Kirst. (2012). Accuracy of genomic selection methods in a standard data set of loblolly pine (_Pinus taeda_ L.) Genetics 190:1503-1510

### **File to download**

[EQ\_LOBLOLLY](https://vsni1.wpengine.com/wp-content/uploads/2020/09/EQ_LOBLOLLY.zip)

**Notes**: SAG August-2020

Most breeding programs plan several controlled crosses between outstanding parents to detect favorable alleles in their offspring. The progeny are later evaluated in a field experiment.

large_2775807.jpg

medium_2775807.jpg

small_2775807.jpg

thumbnail_2775807.jpg

Overlaying Design Matrices for a Parental Model in ASReml-R

Statistics - serving plant science for generations

Predicting the performance of novel plant genotypes in different environments is essential for developing cultivars with superior economically important traits (such as yield) and resilience to environmental stresses (such as drought).

<table style="background-color:rgb(255, 255, 255);"><tbody><tr><td><p><strong>Case Study</strong></p><figure class="image"><img src="https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/512px_Starr_030419_0163_Pinus_radiata_f440050ede.jpg" alt="radiata pine cones"></figure><p>The <a href="https://www.rpbc.co.nz/">Radiata Pine Breeding Company Ltd.”</a> runs breeding programmes to improve the productivity and quality of wood from radiata pine trees for the Australasian forestry industry. Candidate pine trees (i.e. different genotypes) are tested in a <span style="color:rgb(0,0,255);"><strong>multi-environment trial</strong></span> (MET), a series of experiments conducted across a range of geographic locations in New Zealand and Australia over multiple years. The set of experiments within and across years is designed to provide a range of growing conditions, or “<span style="color:rgb(0,0,255);"><strong>environments</strong></span>”. From these trials, the genetic and environmental effects on the performance of the candidate pine trees can be estimated. This information is then used by radiata pine breeders to select trees for crossing, with the aim of further genetic improvement, and by growers of radiata pine for selecting genotypes predicted to perform well at their site.</p></td></tr></tbody></table>

Plant breeders use multi-environment trial (MET) data to evaluate genotypes across a range of environments. However, the relative performance of the genotypes often varies between environments, a phenomenon known as the genotype by environment (G×E) interaction. The G×E interaction can be exploited to identify genotypes that perform well in all environments (i.e. are suitable for broad use) and those with exceptional performance in specific environments (i.e. are well suited for use in certain growing conditions).

Today, linear mixed models are widely used in the analysis of MET data. The linear mixed model framework accommodates the analysis of genetically and/or experimentally correlated data, heterogeneous variances and unbalanced data sets, enabling the accurate prediction of genotype performance within all environments in the data set. In addition, the ability to formally test statistical hypotheses provides greater insight into the nature of the G×E interaction.

An established method for analysing MET data involves a linear mixed model that adopts multiplicative factor analytic (FA) variance structures for the G×E effects. The FA model provides a parsimonious, yet flexible, method of describing the G×E interaction. For example, it allows for genetic variance heterogeneity between trials and different genetic correlation across trials. It can also be extended to include genetic relationship information (e.g. pedigrees) so that the genetic effects can be partitioned into additive and non-additive components. Furthermore, it typically has higher predictive accuracy than alternative models when there is substantial G×E interaction.

<table style="background-color:rgb(255, 255, 255);"><tbody><tr><td><p><strong>A Brief History of MET Analysis</strong></p><ul><li>Early methods used ANOVA on the two-way table of genotype by environment means. Here, the total variation is partitioned into sources due to genotype, environment and residual variation (a combination of the G×E interaction and within-trial error). Estimates of the average genotype performance across environments can be obtained.<br><span style="color:rgb(255,0,0);"><i><strong>Limitation:</strong></i></span> Doesn’t provide information on the nature of the G×E interaction.</li><li>A greater emphasis on understanding the G×E interaction lead to the development of the AMMI and GGE biplots; these are descriptive tools for visualising the relationships between genotypes and environments.<br><span style="color:rgb(255,0,0);"><i><strong>Limitation:</strong></i></span> Doesn’t provide simple numerical summaries useful for plant selection.</li><li>Today, linear mixed models are widely used for the analysis of MET data. In particular, a linear mixed model with a multiplicative factor analytic (FA) model for the G×E effects.&nbsp; This has been found to perform extremely well in terms of parsimoniously describing the G×E interaction and predictive accuracy.</li></ul></td></tr></tbody></table>

A powerful, efficient and reliable software solution for fitting linear mixed models is [ASReml-R](https://www.vsni.co.uk/software/asreml-r). It is particularly well suited to the analysis of MET data. With over 4500 citations, ASReml-R is a popular choice by MET data analysts due to:

*   flexible syntax that makes fitting complex linear mixed models possible and easy
*   an efficient algorithm for fitting the linear mixed model, which makes it feasible to analyse large and complex data sets.

<table style="background-color:rgb(255, 255, 255);"><tbody><tr><td><p><strong>Case Study Continued…</strong></p><p>Biometricians from the University of Wollongong, Australia, have <a href="https://uploads-ssl.webflow.com/5d4797a2a904cddbbb041f99/5d7ad3db34aefeac46295f95_RPBC%202019%20BV%20Analysis%20Report.pdf">analysed MET data</a> for the Radiata Pine Breeding Company Ltd. Over 360,000 progeny trees were grown in 111 trials. Pedigree information was available for the 2948 parental lines. Their aim was to produce predicted breeding values (i.e. additive genetic effects) of the parental trees. A linear mixed model adopting a fa ctor analytic (FA) model for the G×E effects was fitted using <a href="https://www.vsni.co.uk/software/asreml-r">ASReml-R</a>.</p><p>The breeding values predicted by the statistical model are now publicly available on the <a href="https://www.rpbc.co.nz/">Radiata Pine Breeding Company Ltd.</a> website for use by radiata pine breeders and growers.</p><figure class="image"><img src="https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/RPBC_homepage_410eff9ef7.png" alt="Radiata Pine Breeding Company home page"></figure></td></tr></tbody></table>

For readers new to linear mixed models or ASReml-R, [MMA (Mixed Model Academy)](http://www.mixedmodel.academy/) provides a simple and user-friendly tool to help you

*   formulate an appropriate linear mixed model for your data set
*   learn ASReml-R coding syntax
*   connect with a community of linear mixed model practitioners.

Predicting the performance of plant cultivars

large_Photosynthesis_2880_Lede (1)-min.jpg

medium_Photosynthesis_2880_Lede (1)-min.jpg

small_Photosynthesis_2880_Lede (1)-min.jpg

thumbnail_Photosynthesis_2880_Lede (1)-min.jpg

An important aim when fitting linear mixed models (LMM) is the use of random effect estimates. In some analyses, such as genetic evaluations, the main objective of the analysis is to obtain these estimates. Random effects, also known as BLUPs (best linear unbiased estimations), are obtained after fitting a LMM, where their variance components are estimated by residual maximum likelihood (REML) and these values are used to calculate BLUPs.

Consider the following LMM for a genetic experiment, where we have a randomized complete block design with $b$ blocks and a total of $t$ genotypes evaluated:

$$

y\_{ij} = \\mu + \\beta\_j + a\_i + \\epsilon\_{ij}

$$

where $y{ij}$ is the phenotypic observation of the ith genotype on the jth block; $\\mu$ is the overall mean; $\\beta$ is the fixed effect of the ith block; $a\_i$ is the random effect of the ith genotype, with $a\_i \\sim N(0, \\sigma\_a^2)$ that is associated with some information about pedigree, from which we can obtain a pedigree-based or genomic-based relationship matrix $\\mathbf{A}$; and $\\epsilon\_\\\_{ij}$ is the random residual with $\\epsilon\_i \\sim N(0, \\sigma^2)$.

In the above model, also known as Animal Model, the $a\_i$ random effects correspond to the breeding value (BV) from which selections of future progenitors or individuals will be done. As these effects are based on estimated variance components, then they are often known as EBVs (estimated breeding values). Statistically, the matrix $\\mathbf{A}$ is critical as it is the place were we specify the _correlations or dependencies_ that exist between random effects.

As many breeding and commercial decisions depend on the EBVs it is important to have an assessment of our _confidence_ in these values, and this is done by calculating a statistical measure of precision. Probably the most common statistic used to assess these effects is the calculation of the _correlation between true and predicted random effects_, $r(a,\\hat{a})$, also known as _accuracy_. This can be calculated using the following expression for a given genotype as:

$$

r(a\_i,\\hat{a\_i})=\\sqrt{1-\\frac{PEV(\\hat{a\_i})}{\\hat{\\sigma\_{a}^2}}}

$$

This formula requires the PEV (predictor error variance) of $\\hat{a\_i}$, which corresponds to the square of the standard error (SE) of the random effect estimate, _i.e. $PEV(\\hat{a\_i}) = SE(\\hat{a\_i})^2$_. Most statistical software will report these SE values, including SAS and ASReml. In addition, we have $\\hat{\\sigma\_{a}^2}$, which is the population estimate of the variance associated with the EBVs, and this is our genetic additive variance, which is used to calculate narrow-sense heritability ($\\hat{\\sigma\_{a}^2}$).

According to statistical theory, the values of PEV will range from zero to the ${\\sigma\_{a}^2}$. In cases where we have full information about a given genotype, then the PEV will be close to zero as we have little uncertainty in the true additive value of a genotype, hence,  ${\\sigma\_{a}^2}$ . However, when there is very little information, then the PEV will approximate to ${\\sigma\_{a}^2}$. Note that these extreme values of PEV translate in $r(a,\\hat{a})$ values of 1.0 and 0.0, respectively. Therefore, values closer to 1.0 are an indication of very good quality of a given additive effect estimate.

Sometimes, you will see _reliability_ instead of _accuracy_ reported, where the former is the square of the later. Both statistics are commonly used when reporting the quality of a random effect estimate in genetic studies.

$$

r(a\_i,\\hat{a\_i})^2={1-\\frac{PEV(\\hat{a\_i})}{\\hat{\\sigma\_{a}^2}}}

$$

To illustrate the above calculations, we used genetic data originating from a trial on a population of cultivated two-row spring barley. The original trial measured height (cm) of plants grown in pots during 2011 for a total of 856 lines; however, the available data for this analysis consists of the adjusted mean values of only 478 lines, together with their molecular information. Further details can be found in Oakley _et al_. (2016).

In our example, we obtained the realized relationship matrix from the single nucleotide polymorphism (SNP) marker information available for these 478 lines, and an Animal Model was fitted considering this genomic matrix. This was done in the package R using the library ASReml-R, and the code used to fit this model was:

```plaintext
modelGBLUP &lt;- asreml(fixed=Hadj11~1,
                     random=~vm(lines, ainv),
                     na.action=na.method(y='include'),
                     data=datagFULL)
```

Here, the response variable was $\\texttt{Hadj11}$ and our single random effect is the factor $\\texttt{lines}$ that represents our additive effects, and this is associated with the inverse of the relationship matrix $\\texttt{ainv}$. In addition, in the above code it was requested that observations with missing data (_i.e._, $\\texttt{NA}$) in their response variable will be included in the analyses. This allowed us to obtain EBVs for all individuals, even those that were not phenotyped but that were included in the relationship matrix.

The above analysis, using the following code

```plaintext
vc &lt;- summary(modelGBLUP)$varcomp
```

will output the following estimation of variance components.

|   | component | std.error | z.ratio | bound | %ch |
| --- | --- | --- | --- | --- | --- |
| vm(lines, ainv) | 34.94159 | 7.880138 | 4.434134 | P | 0.4 |
| units!R | 43.10000 | 7.949305 | 5.421857 | P | 0.0 |

These values correspond to a narrow-sense heritability $h^2$ of $0.478 (\\pm 0.097)$ for the adjusted mean of height. Hence, the trait $\\texttt{Hadj11}$ has a strong genetic control. Now, we can proceed to request the BLUP (or EBV) values together with their standard errors, this is done using the code:

```plaintext

BLUP &lt;- summary(modelGBLUP, coef=TRUE)$coef.random
```

Below we can see a subset of this list with the top 10 and bottom 10 genotypes sorted by their SE values.

|   | solution | std.error | z.ratio | PEV | accuracy | reliability |
| --- | --- | --- | --- | --- | --- | --- |
| vm(lines, ainv)\_LENTA | 5.8347 | 3.9482 | 1.4778 | 15.5885 | 0.7442 | 0.5539 |
| vm(lines, ainv)\_PALLAS | 5.1138 | 3.9807 | 1.2846 | 15.8462 | 0.7393 | 0.5465 |
| vm(lines, ainv)\_CHARIOT | 1.1494 | 3.9859 | 0.2884 | 15.8871 | 0.7385 | 0.5453 |
| vm(lines, ainv)\_CLASS | \-7.6429 | 3.9874 | \-1.9168 | 15.8990 | 0.7382 | 0.5450 |
| vm(lines, ainv)\_BONUS | 7.9205 | 3.9897 | 1.9853 | 15.9176 | 0.7379 | 0.5445 |
| vm(lines, ainv)\_TOKEN | \-7.8575 | 3.9930 | \-1.9678 | 15.9440 | 0.7374 | 0.5437 |
| vm(lines, ainv)\_MAJA | 9.9281 | 4.0022 | 2.4806 | 16.0178 | 0.7359 | 0.5416 |
| vm(lines, ainv)\_PRESTIGE | \-8.4205 | 4.0184 | \-2.0955 | 16.1477 | 0.7334 | 0.5379 |
| vm(lines, ainv)\_ASPEN | \-4.1411 | 4.0456 | \-1.0236 | 16.3670 | 0.7291 | 0.5316 |
| vm(lines, ainv)\_CORNICHE | \-1.6122 | 4.0546 | \-0.3976 | 16.4394 | 0.7277 | 0.5295 |

|   | solution | std.error | z.ratio | PEV | accuracy | reliability |
| --- | --- | --- | --- | --- | --- | --- |
| vm(lines, ainv)\_FREJA | 28.0177 | 4.8260 | 5.8056 | 23.2902 | 0.5775 | 0.3335 |
| vm(lines, ainv)\_DASIO | \-8.5969 | 4.8376 | \-1.7771 | 23.4028 | 0.5747 | 0.3302 |
| vm(lines, ainv)\_KESTREL | \-0.4303 | 4.8644 | \-0.0885 | 23.6623 | 0.5682 | 0.3228 |
| vm(lines, ainv)\_GAIRDNER | 5.3170 | 4.9188 | 1.0810 | 24.1944 | 0.5546 | 0.3076 |
| vm(lines, ainv)\_STELLA | 15.1422 | 4.9273 | 3.0731 | 24.2781 | 0.5524 | 0.3052 |
| vm(lines, ainv)\_LOGAN | 5.7601 | 4.9716 | 1.1586 | 24.7170 | 0.5409 | 0.2926 |
| vm(lines, ainv)\_CLARA | 9.4874 | 5.3870 | 1.7612 | 29.0197 | 0.4117 | 0.1695 |
| vm(lines, ainv)\_MARESI | 1.5449 | 5.7448 | 0.2689 | 33.0032 | 0.2355 | 0.0555 |
| vm(lines, ainv)\_FLICK SEJET | \-8.5443 | 5.8523 | \-1.4600 | 34.2490 | 0.1408 | 0.0198 |
| vm(lines, ainv)\_BEKA | 16.5908 | 6.0992 | 2.7202 | 37.2001 | NaN | \-0.0646 |

The EBV estimates for this data are positive and negative, as they represent deviations from the overall height mean. For example, the genotype FREJA has an EBV of $+28.02$ indicating that, if this parent is crossed with a parent of equal genetic worth, its progeny will have an expected deviation from the mean of $+28$ cm. The full range of EBVs is between −$14.54$ cm to $43.30$ cm, a wide interval considering that the mean of the population is $88.03$ cm, and this is a result of the relatively large heritability value found for this trait.

Also, there is an interesting range in SE values from $3.949$ to $6.105$. The width of this range is related to the level of genetic relationships between individuals, that is specified through the matrix $\\mathbf{A}$. Genotypes with many relatives will _share_ or _collect_ information from many phenotypic responses, and therefore they will have lower SE values; in contrast, individuals with few, or no relatives, will tend to have much higher SE values. However, it is also possible that individuals with no phenotypic data will have large SE values, unless they have many relatives to help with the estimation of their EBVs.

Now, we can proceed to calculate the accuracy and reliability using the expressions shown before, with the code

```plaintext
Va <- vc[1,1]
BLUP$PEV <- BLUP$std.error^2
BLUP$accuracy <- sqrt(1 - BLUP$PEV/Va)
BLUP$reliability <- 1 - BLUP$PEV/Va
```

And now we can observe these calculated values for a subset of the top 10 and bottom 10 individuals.

|   | PEV | accuracy | reliability |
| --- | --- | --- | --- |
| vm(lines, ainv)\_LENTA | 15.5885 | 0.7442 | 0.5539 |
| vm(lines, ainv)\_PALLAS | 15.8462 | 0.7393 | 0.5465 |
| vm(lines, ainv)\_CHARIOT | 15.8871 | 0.7385 | 0.5453 |
| vm(lines, ainv)\_CLASS | 15.8990 | 0.7382 | 0.5450 |
| vm(lines, ainv)\_BONUS | 15.9176 | 0.7379 | 0.5445 |
| vm(lines, ainv)\_TOKEN | 15.9440 | 0.7374 | 0.5437 |
| vm(lines, ainv)\_MAJA | 16.0178 | 0.7359 | 0.5416 |
| vm(lines, ainv)\_PRESTIGE | 16.1477 | 0.7334 | 0.5379 |
| vm(lines, ainv)\_ASPEN | 16.3670 | 0.7291 | 0.5316 |
| vm(lines, ainv)\_CORNICHE | 16.4394 | 0.7277 | 0.5295 |

|   | PEV | accuracy | reliability |
| --- | --- | --- | --- |
| vm(lines, ainv)\_FREJA | 23.2902 | 0.5775 | 0.3335 |
| vm(lines, ainv)\_DASIO | 23.4028 | 0.5747 | 0.3302 |
| vm(lines, ainv)\_KESTREL | 23.6623 | 0.5682 | 0.3228 |
| vm(lines, ainv)\_GAIRDNER | 24.1944 | 0.5546 | 0.3076 |
| vm(lines, ainv)\_STELLA | 24.2781 | 0.5524 | 0.3052 |
| vm(lines, ainv)\_LOGAN | 24.7170 | 0.5409 | 0.2926 |
| vm(lines, ainv)\_CLARA | 29.0197 | 0.4117 | 0.1695 |
| vm(lines, ainv)\_MARESI | 33.0032 | 0.2355 | 0.0555 |
| vm(lines, ainv)\_FLICK SEJET | 34.2490 | 0.1408 | 0.0198 |
| vm(lines, ainv)\_BEKA | 37.2001 | NaN | \-0.0646 |

For this trait, the range of correlation between true and predicted breeding values, $r(a,\\hat{a})$, is from $0.134$ to $0.744$, and the lowest accuracy and reliability values correspond to individuals with the largest SE values. Note that the genotype $\\texttt{BEKA}$ has a $\\texttt{NaN}$ for accuracy, which originates from a negative square root, as this BLUP has an SE that is larger than $\\hat{\\sigma\_{a}^2}$.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/reliability_image_2_fbd9816e68.png)

In the above figure we present a boxplot of the accuracy values (left plot) and a scatterplot of the PEV against accuracy (right plot). It is clear from these plots that there is a group of three observations with very poor accuracy that have $r(a,\\hat{a})$ values lower than $0.41$. However, the majority of the observations are found with accuracy values of $0.54$ or higher, where some of them reach values as high as $0.74$.

It is interesting to speculate what the reasons could be for these three lower observations (genotypes $\\texttt{CLARA}$, $\\texttt{MARESI}$ and $\\texttt{FLICK}$ $\\texttt{SEJET}$) together with the $\\texttt{NaN}$ observation (genotype $\\texttt{BEKA}$), to be so different from the rest of the population. The unusual behavior of these four observations is likely to reflect errors in their information, including phenotypic response (maybe a potential outlier), or incorrect genetic relationships (due to pedigree or genotyping errors).

In summary, we have calculated values of accuracy and reliability for the estimation of random effects from an Animal Model. For the data analyzed, these represented a wide range of values, reflecting different levels of precision. For example, a large value of accuracy $r(a, \\hat{a})$, such as $0.90$, will indicate that there is little expected change on the estimation of the EBV if additional information is collected from a given genotype or from its relatives (_e.g._, its offspring). Therefore, $r(a, \\hat{a})$ values provide us with a good measure of our confidence in the genetic evaluation performed.

In our analysis, and in many other genetic analyses, it is critical to consider the incorporation of the matrix $\\mathbf{A}$ of genetic relationships, as it will allow us to improve the estimation of the random effects because we have changed our assumption from independent random effects to a model that considers correlations or dependencies between these random effects.

A routine check of the accuracy and reliability (and indirectly of PEV values), can help to make future improvements on the quality of the information obtained from statistical analyses. For example, we can assess the consequences of incorporating additional relatives or replicates in future evaluations. Alternatively, these checks can also help to assess the benefits of using more complex linear models, such as those that combine multiple traits, measurements or sites.

Finally, the above calculations can be extended to any other random effect on all linear mixed models fitted, and their use is recommended to assess the quality of the estimates beyond what a simple standard error will indicate.

## Author

Salvador A. Gezan

## Reference

Oakey, H; Cullis, B; Thompson, R; Comadran, J; Halpin, C; Waugh, R. 2016. _Genomic selection in multi-environment crop trials_. Genes | Genomes | Genetics 6:1313-1326

**Notes**: SAG Jan-2020

Calculating Accuracy and Reliability of Random Effect Estimates with ASReml-R

large_https://wholegrainscouncil.org/sites/default/files/thumbnails/image/15616367_ml-wheat-field_0.jpg

medium_https://wholegrainscouncil.org/sites/default/files/thumbnails/image/15616367_ml-wheat-field_0.jpg

small_https://wholegrainscouncil.org/sites/default/files/thumbnails/image/15616367_ml-wheat-field_0.jpg

thumbnail_https://wholegrainscouncil.org/sites/default/files/thumbnails/image/15616367_ml-wheat-field_0.jpg

I have been an avid supporter of Genstat since starting my PhD in 1985 (gulp!). As an applied plant biologist, it is perhaps unsurprising that I found procedures that had been developed with colleagues at Rothamsted Research to be highly relevant to the analysis of my own designed experiments. During my thirty-five years in research, teaching and management at higher education institutions, I have continued to rely on Genstat as my core statistical programme both as an investigator and an educator.

My own research is in the area of crop science, and in the last decade or so, justified (to research funders at least) on the basis of the need to ‘sustainably intensify’ our food production \[1\]. For many the phrase is an imperative, i.e. we have to use our land more efficiently and it will have to keep on being able to do this; for others the phrase is an oxymoron, e.g. if land-use efficiency means producing ever higher yields per unit area this will inevitably degrade the very ecosystems and resources on which agricultural production by future generations will depend; a third group would argue that sustainable intensification is unnecessary – if (big if) we were each satisfied with the per capita calorific availability that was the average for humans in 1960, we are already producing enough calories to feed 10 billion people. Hopefully my research has been relevant whichever ‘camp’ you may belong to.

In my early days as a researcher I was primarily interested in the optimising of rates of synthetic inputs (such as fertilizers and crop protection chemicals) to achieve yield and grain quality targets in wheat, to improve economic and environmental sustainability. Experiments were usually complete multi-factorial designs (genetic x management), replicated in randomised blocks usually with some split plotting. Given Genstat’s origins, it is unsurprising that the treatment structure and block structure routines allow a very large degree of flexibility to cope admirably with all types of designed field experiment. Of particular use to me was the ability to split treatment effects into polynomial contrasts within ANOVA \[2, 3, 4\]. ‘Life is not a polynomial’ but this approach provides an initial, statistically sensitive test as to whether there is a response to increasing rates of an input, and some idea as to the complexity of that response. Guided by the ANOVA, the FITCURVE routine provides a straightforward step for exploring more biologically meaningful and parsimonious fits to the responses to inputs \[5,6\]. The standard curves available in Genstat are highly relevant to the applied biologist but further development and modification of responses is supported by FITNONLINEAR \[7,8\].  FITCURVE and FITNONLINEAR have also been useful in fitting responses over time.

Often treatment (such as fungicide) effects on disease or green leaf area on individual days are associated with effects on the errors around the treatment means. This variance heterogeneity is easily identified in the residual plots provided by ANOVA and can be rectified by transformations. Transformed data, however, can be difficult to interpret, whereas fitting responses over time (e.g. logistic or gompertz) can generate new variables such as time scalars which are readily interpretable, and for which the residuals are often normally distributed and lacking variance heterogeneity \[7\].

Over the years the teams of PhD students and researchers I’ve been involved with have continued to be heavily reliant on ANOVA, FITCURVE and FITNONLINEAR although the responses these days are more likely to be due to climate variables (genetic x environment) \[9\]. As datasets and experiments have got bigger and less complete, residual maximum likelihood (REML) has proved more flexible \[10, 11, 12\] than ANOVA, and principal components (PCP) analysis ably condenses large variate x experimental unit matrixes \[13,14\]. The current portfolio of procedures is well matched for advances in applied biology whether involving meta-analyses, ‘big data’, statistics to underpin machine learning, and/or precision agriculture.

As an educator I have taught BSc and MSc statistics and research methods modules with Genstat, and have enabled level 5, 6, & 7 students using the programme to assess class-generated data. The interface is intuitive, and it is easy to demonstrate variations in approaches to highlight the importance of, for example, replication (and the curse of pseudo-replication), blocking structures, covariance, randomisation and the like. Genstat has been the main statistical programme used by the 40 PhD students who have been gracious enough to allow me to be involved in their supervision.

   
1 Reaping the beneﬁts: science and the sustainable intensiﬁcation of global agriculture ISBN: 978-0-85403-784-1;  
2_J. Agric. Sci_ **138**, 317-331, doi: 10.1017\\S0021859602002137;  
3_J. Cereal Sci._ **37**, 295-309, doi: 10.1006/jcrs.2002.0501;  
4_Plant Soil_ **360**, 93-107, doi: 10.1007/s11104-012-1203-x;  
5_J. Sci. Food Agric._ **85**, 727-742, doi: 10.1002/jsfa.2025;  
6_J. Agron. Crop Sci._ **200**, 36-45, doi: 10.1111/jac.12038;  
7_Annals Appl. Biol._ **136**, 77-84, doi: 10.1111/j.1744-7348.2000.tb00011.x;  
8_Field Crops Res._ **95**, 49-63, doi: 10.1016/j.fcr.2005.02.001;  
9Frontiers Plant Science **8**:51, doi: 10.3389/fpls.2017.00051;  
10_J Agric Sci_ **128**, 135-142, doi: 10.1017/S0021859696004054;  
11_Agric. Forest Meteor._ **94**, 159-170, doi: 10.1016/S0168-1923(99)00020-9;  
12_PLoS One_ **11**(5), e0156056, doi: 10.1371/journal.pone.0156056;  
13_J. Sci. Food Agric._ 84, 227-236, doi: 10.1002/jsfa.1657;  
14Euphytica **166**, 249-263, doi: 10.1007/s10681-008-9838-7

Mike Gooding FRSB

(Mike has headed agriculture and been Professor of Crop Science at The University of Reading, Aberystwyth University, and the Royal Agricultural University).

Using Genstat in the ‘sustainable intensification’ of food production

large_spike-wheat-cereals-grain-546332-scaled.jpg

medium_spike-wheat-cereals-grain-546332-scaled.jpg

small_spike-wheat-cereals-grain-546332-scaled.jpg

thumbnail_spike-wheat-cereals-grain-546332-scaled.jpg

Aquaculture is the name given to the farming of salt water and freshwater fish and marine animals. It is a rapidly growing division of farming across the world; reports from the UN’s FAO suggest that it has been growing more rapidly than any other area of animal food production.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/Nofima_AS_Reml_6e610934f8.jpg)

A report published in September 2009 suggested that farm-reared fish accounted for 50% of world fish consumption.

This rise in aquaculture production isn’t much of a surprise when you consider the emphasis on the health-giving properties of fish (particularly Omega-3 from oily fish) and the subsequent encouragement for people to eat more oily fish; and it’s certainly good news for the aquaculture producers across the world. However, growth of this level in any industry, whilst exciting and positive in many ways, also brings with it its own complications. The complexities of all types of farming mean that research and understanding is vital in helping the industry to grow safely as well as rapidly. This is particularly important if the current consumption levels of fish continue: the FAO predict that we will need an extra 27 million tonnes of fish by 2030 to meet the demand.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/shutterstock_52746229_1024x762_e96d0452b2.jpg)

Any farmer will tell you of the importance of vaccination in controlling disease, and this is no different for fish farmers, or specifically Atlantic salmon farmers. Vaccinations against diseases such as furunculosis, vibriosis, cold-water vibriosis, winter ulcer and infectious pancreatic necrosis, which could otherwise decimate a farm, are vital in protecting the fish and farms and controlling these diseases in farmed Atlantic salmon.

However, vaccinations can cause reactions at the injection site such as adhesions and discoloration, reduced appetite and growth, and in severe cases a downgrading of the carcass at market. Not only can this have an economic impact on the salmon farms, but also may create a negative perception of the farms with the general public. As with all aspects of farming the public is concerned not just with the food in question, but also of the welfare of the production animals, which means that breeding companies need to be sure of the actual physical health and safety of their fish.

As with many aspects of farming and agriculture, all aspects of aquaculture are the subject of scientific research to enable better understanding and improved planning for the future, including prevention of disease and improved fish health. The results from challenge tests have indicated substantial genetic variation in disease resistance to a variety of bacterial and viral diseases (including some of those listed earlier), significant but relatively low genetic correlations among the different diseases and a high genetic correlation between the survival of non-vaccinated pre-smolts (smolts are young salmon as they migrate from freshwater to the sea) and the survival of older non-vaccinated full and half sibs under natural break-outs of the diseases. With this in mind Atlantic salmon breeding companies are looking at both the development of efficient vaccinations and the breeding of fish with higher disease resistance.

The results from the challenge test research may suggest that genetic variation in immune systems might lead to genetic variation in innate disease resistance and variation in reactions to the vaccinations. One question researchers in this area need an answer to is whether there is significant genetic correlation between vaccine reactions and innate resistance, as this could mean a change of vaccine dose or a breeding programme to produce fish with the lowest reaction rate to vaccines and the highest natural resistance to diseases.

And it is these issues that a recent project, undertaken by researchers from [Nofima Marin](http://www.nofima.no/marin/en), the [Norwegian School of Veterinary Scienc](http://www.veths.no/)e and [Aqua Gen AS](http://www.aquagen.no/En/Specialfolder/Latest_News/) and funded by the [Research Council of Norway](http://www.forskningsradet.no/servlet/Satellite?c=Page&cid=1177315753906&p=1177315753906&pagename=ForskningsradetEngelsk%2FHovedsidemal) (published in the Elsevier journal Aquaculture 287 (2009) 52-58) has looked at.

Atlantic salmon were taken from farmed strains and from offspring of wild salmon, and reared at different farms. Data was collected on sex, size (at different stages), adhesion scores and melanin scores of fish at harvest size and survival data from challenge tests of their full- and halfsibs. The data collected was analysed using [ASReml](https://www.vsni.co.uk/software/asreml/) to fit a multivariate linear mixed animal model to obtain estimates for the variance and covariance components for the random effects of the studied traits. A single trait sire and dam model, which included the same fixed effects was used to obtain across year-class heritability estimates (using an alpha version of [ASReml 3](https://www.vsni.co.uk/software/asreml/)).

Results from the project indicate that the severity of vaccine injuries in Atlantic salmon could be reduced through selective breeding. Additionally the relatively high correlation that was found between adhesion scores and melanin scores show that there is unlikely to be a need to select against both of these traits. The researchers warned that including an additional breeding trait objective would reduce the genetic gain for other traits, and therefore recommended that other solutions to the vaccine injuries be found, such as improved vaccines or vaccination procedures. The results also indicated a negative genetic correlation between body weight and adhesion or melanin score; however, these estimates are most likely biased towards the negative end of the scale, since vaccine injuries are expected to reduce growth rate for non-genetic reasons.  The study also indicated “no genetic association between resistance to furunuculosis and ISA (infectious salmon anaemia) and vaccine injuries”, which suggest that the vaccine injuries are caused by other aspects of the immune system, not associated with the genetic variation in innate resistance to these diseases.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/shutterstock_726038359_1024x682_5f93bc7e31.jpg)

Research of this kind is extremely important to the aquaculture industry to enable it to continue to develop, expand and improve, not just the industry, but the procedures, processes and systems associated with it. ASReml is widely respected and used in situations such as this.

Our thanks to [Dr Bjarne Gjerde](http://www.nofima.no/marin/en/person/bjarne.gjerde) of [Nofima Marin](http://www.nofima.no/marin/en) for his help with this story.

large_VSNi_aquaculture_image_desktop.jpg

medium_VSNi_aquaculture_image_desktop.jpg

small_VSNi_aquaculture_image_desktop.jpg

thumbnail_VSNi_aquaculture_image_desktop.jpg

Nofima

It’s not news that elephants are endangered; for literally thousands of years ivory has driven hunters to kill these magnificent beasts. Combine this deliberate destruction with increased land use change and habitat transformation by private enterprise, international aid and government expansion across Africa and these magnificent beasts are still under threat even within the national parks and other protected areas set aside for them.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/mike_genstat_28f70ac84d.png)

However, before we get too carried away with horrific idea that elephants may well become extinct; there is another side to this.

Two separate taxa, or races, of elephants exist in Africa today: the savannah elephant, which is well known and has been studied and valued by the international tourist trade for decades, and the diminutive forest elephant, which until recently has remained in obscurity beneath the trees. In the Congo Basin, an area noted for being one of the great wilderness areas with the second largest rainforest on earth – it is home to an outstanding range of biological diversity; forest elephants have until recent decades thrived. These types of wildernesses have always been recognised for their high conservation value, being home to a wide variety of flora and fauna, and as such provide invaluable resources of timber, minerals and other natural resources. The Congo Basin is a typical case where the resources available are in high demand in rich, industrial nations and are also critical to the economies of poor nations. And it is at this point we in the Western developed world need to be very careful about conservation and conservation techniques; the international conservation efforts are aimed at managing and promoting sustainable natural resource exploitation and functional ecosystems. It is incredibly important that we recognise the needs of people in the developing world, and try to ensure they are able to exploit their natural resources in ways that do not fundamentally alter the ecology of the environment, upon which millions of rural people still depend.

The high value of resources in these wilderness areas means that it is only fair that Africans have an opportunity to exploit and benefit from them and are encouraged to do so legally, ethically and with an eye on their own local ecology.

Complications occur as industrial expansion pushes into the wildernesses to obtain these natural resources, an expansion that relies on a network of roads and rail infrastructure. It is this increase in roads that appears to be causing some of the biggest threats to the forest elephants. Illegal elephant killings tend to occur close to roads, since these areas are easily accessible, so the density of forest elephants in the protected areas is determined, not by the size of the area but by how many roads criss-cross through the wilderness.

Elephants are by their nature, shy and avoid human contact as much as possible, so how are they reacting to the increasing number of roads across their home ranges? Will they adopt a “siege” mentality and reduce their home range to avoid roads, or a “skirmish” strategy, where they will continue to range widely, irrespective of the number of roads?

I spoke with Stephen Blake of the Wildlife Conservation Society and Max Planck Institute of Ornithology who has recently headed up a study into this very issue. His team fitted a large sample of forest elephants in the Congo Basin with GPS telemetry collars, so that they could track their roaming behaviour. They collected data from 28 elephants located in conservation areas in the Congo, Central African Republic and Gabon.

Figures 1 & 2 below show data collected from 28 elephants located in the Congo, Central African Republic and Gabon using Genstat.

Figure 1: Relationship between area of MCP home range and the area of road-less polygon; y axis – Home range size (kmsq), x axis roadless polygon size (kmsq).

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/Figure1_1024x966_f96a7a03fc.jpeg)

Figure 2: Maximum linear displacement versus area of road-less polygon; y axis maximum linear distance moved (km), x axis roadless polygon size (kmsq).

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/Figure2_1024x948_1c6fc46651.jpeg)

“The data we collected covered the sex of the elephant, the site (the protected area in which the elephants were collared), the number of roads in the area, the number of days the elephants were collared,” explained Stephen. “I initially analysed the result collected in Genstat, using linear regression for data exploration followed by fitting standard curves to the data (Figures 1 and 2). Genstat was vital in the research programme as I can trust the statistical procedures. I was able to talk directly with developers when I came up against any queries. I have found the Genstat support team to be extremely helpful, I usually get a response to a query within 24 hours and often much sooner. We could have gone much further with the analysis, but in order to understand the influence of roads on forest elephant ranging, the regression options in the Genstat menu showed us clear statistical patterns.”

Results indicated that the forest elephants routinely crossed protected roads, however only one elephant crossed an unprotected road, and of these one did so, as far away from the nearest village as she could, and her speed increased as she crossed the road.

“Our study showed that elephants seem to adopt a siege mentality with the encroaching roads”, said Stephen Blake, “and the issue with this is that over time the elephants will have less and less access to food, in quantity and variety, leading to food competition and aggression; both of which are likely to lead to a reduction in population.”

“The main issue seems to be that elephants are risk averse, hence why they avoid roads where they can. So here we have a problem; central African states need the income earned from the natural resources found in their countries, but not only do we not want to see the elephants themselves disappear, a change like that to any ecosystem is likely to have more wide ranging impacts.”

“There are options for policies to overcome the issues, before it is too late,” says Stephen. “Based on our research, we recommend that no new roads are built on existing forest elephant conservation areas, but to access the natural resources near these sites temporary roads should be built, so that once the resources are exhausted the area can easily return to natural wilderness. The other option is to try to reduce the elephants’ fear of roads, by preventing settlements on the roads, enforcing anti-poaching. These are of course expensive options, but without addressing this we face losing these beautiful creatures.”

large_Monitoring_elephants_mainimage_desktop.jpg

medium_Monitoring_elephants_mainimage_desktop.jpg

small_Monitoring_elephants_mainimage_desktop.jpg

thumbnail_Monitoring_elephants_mainimage_desktop.jpg

MIKE Monitoring Elephants

BirdLife attempted to identify and record all present waders and coastal birds. From there they can identify the areas that need protection and preservation or even rehabilitation.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/birdlife_genstat_1f88f48536.png)

#### **What was their goal?**

Their aim was to raise awareness of the decreasing populations in coastal birds and waders. They hope this will help the protection of the estuarine and tidal ecosystems around Australia. BirdLife attempted to identify and record all present waders and coastal birds. From there they can identify the areas that need protection and preservation or even rehabilitation.

#### **What did they use Genstat for?**

Generalised Linear Models (GLM) were used to determine if statistically significant trends in the maximum counts were occurring. Count numbers were regressed against years to detect trends. Analysis was undertaken using Genstat 13.0 software. Data was assessed for fit to assumptions of GLM (normality, spread etc). For species that were found to greatly violate these assumptions, log transformation was undertaken. Transformations of data were undertaken in Genstat.

large_Shorebirds_mainimage_desktop_flat.jpg

medium_Shorebirds_mainimage_desktop_flat.jpg

small_Shorebirds_mainimage_desktop_flat.jpg

thumbnail_Shorebirds_mainimage_desktop_flat.jpg

Birdlife International

Dealing with Pseudo-Replication in Linear Mixed Models

large_cannock-chase-2_42512571184_o-min (1).jpg

medium_cannock-chase-2_42512571184_o-min (1).jpg

small_cannock-chase-2_42512571184_o-min (1).jpg

thumbnail_cannock-chase-2_42512571184_o-min (1).jpg

Often, after fitting a complex linear mixed model (LMM), we are interested in constraining some of the variance components (VCs). There are many reasons why we might want to do this, such as when we combine data from several sources and we want to have a common parameter, or when we want to simplify a model to evaluate a specific hypothesis relating to these components.

In this snippet we will briefly explore an example where we constrain two VCs to be identical. The example comes from Becker (1984) and corresponds to data from 28 white pine (_Pinus monticola_) families formed by crossing 4 male to 7 female parents. The progeny was evaluated in a randomized complete block design with 4 replicates and planted in large plots. Our response of interest was epicotyl length, the part of the embryonic seedling stem from which the plant’s shoot system develops.

Using ASReml-R, the reference model $\\texttt{(mod0)}$ fitted to this data is:

```plaintext
mod0 <- asreml(fixed=length~rep,
               random=~female + male + female:male,
               data=pine)
```

In the above model we have a fixed effect of replicate ($\\texttt{rep}$), and random effects for the parents and their interaction ($\\texttt{female}$, $\\texttt{male}$ and $\\texttt{female:male}$). This LMM produces the following output in ASReml-R:

|   | component | std.error | z.ratio | bound | %ch |
| --- | --- | --- | --- | --- | --- |
| male | 0.1203614 | 0.1203723 | 0.9999094 | P | 0 |
| female | 0.0209448 | 0.0420690 | 0.4978693 | P | 0 |
| female:male | 0.1368513 | 0.0628149 | 2.1786443 | P | 0 |
| units!R | 0.2004263 | 0.0314941 | 6.3639416 | P | 0 |

In this example, the variance of males ($\\hat{{\\sigma\_m}}^2$ = 0.120) is approximately six times larger than the variance for females ($\\hat{{\\sigma\_f}}^2$ = 0.021). Therefore, the additive variance from males and females are contrastingly different: $\\texttt{V}\_{Am} = 4 \\space \\sigma\_m^2 \\neq 4 \\space \\sigma\_f^2 = \\texttt{V}\_{Af}$.

The above results may seem biologically unreasonable for pines, where the genetic contribution of male and female parents (_i.e_., $\\texttt{V}\_{A}$) is often very similar. Hence, we expect their genetic variances to be almost identical. That is, we think $\\sigma\_f^2 = \\sigma\_m^2$. The inconsistent VC estimates from our reference model could be the result of random sampling of families, or to the limited number of parents considered in this dataset (recall that we have only 4 male and 7 female parents).

In ASReml-R, the constraint $\\sigma\_f^2 = \\sigma\_m^2$ is easy to implement. This is done by creating a matrix, say $\\texttt{M}$, that defines the _grouping_ of the VCs. All VCs in the same group (_i.e_., with the same number) will be constrained to be identical. For our problem, the matrix looks like:

```plaintext
##             V1 V2
## female       1  1
## male         1  1
## female:male  2  1
## units!R      3  1
```

The first column, $\\texttt{V1}$, indicates the _grouping_; hence, there are three VCs two of which ($\\texttt{female}$ and $\\texttt{male}$) will be identical. The second column, $\\texttt{V2}$, is used to specify the multiplicative relationships between components, but this is a topic of another ASReml-R snippet.

The construction of the above matrix $\\texttt{M}$ can be facilitated with the following R code. This starts by requesting the table of VC parameters ($\\texttt{gam}$) from ASReml-R:

```plaintext
modC <- asreml(fixed=length~rep,
               random=~female + male + female:male,
               start.values=TRUE,
               data=pine)
gam <- modC$vparameters.table
```

It then proceeds to generate our matrix $\\texttt{M}$ of constraints using only the column names from the table $\\texttt{gam}$:

```plaintext
M <- as.matrix(data.frame(V1=c(1,1,2,3), V2=c(1,1,1,1)))
dimnames(M)[[1]] <- gam$Component
```

Now, we are ready to fit the model in ASReml-R with the constraint $\\sigma\_f^2 = \\sigma\_m^2$. The matrix $\\texttt{M}$ is incorporated using the option $\\texttt{vcc=M}$:

```plaintext
modC <- asreml(fixed=length~rep,
               random=~female + male + female:male,
               vcc=M,
               data=pine)
```

The output from fitting this constrained model in ASReml-R is:

|   | component | std.error | z.ratio | bound | %ch |
| --- | --- | --- | --- | --- | --- |
| male | 0.0639041 | NA | NA | C | 0.5 |
| female | 0.0639041 | 0.0369027 | 1.731692 | P | 0.5 |
| female:male | 0.1327523 | 0.0602139 | 2.204676 | P | 0.0 |
| units!R | 0.2005090 | 0.0312342 | 6.419535 | P | 0.0 |

In the above table see that $\\hat{{\\sigma\_m}}^2 = \\hat{{\\sigma\_m}}^2$ = 0.064. Also note that the z-ratio for this variance component (1.731) is much higher than the z-ratios for both $\\texttt{male}$ (1.000) or $\\texttt{female}$ (0.498) obtained from our reference model (from the first table), indicating better use of the available information.

We compare our constrained model ($\\texttt{modC}$) against the reference model ($\\texttt{mod0}$) statistically, using a likelihood ratio test (LRT). This is implemented using the function ($\\texttt{lrt.asreml}$):

```plaintext
lrt.asreml(modC, mod0, boundary=FALSE)
```

|   | df | LR-statistic | Pr(Chisq) |
| --- | --- | --- | --- |
| mod0/modC | 1 | 1.03068 | 0.309999 |

The LRT indicates that our null hypothesis of $\\texttt{H}\_0: \\sigma\_f^2 = \\sigma\_m^2$ cannot be rejected in favor of our alternative hypothesis $$\\texttt{H}\_1: \\sigma\_f^2 \\neq \\sigma\_m^2$ (p-value = 0.31). That is, despite the estimated variance for $\\texttt{males}$ being approximatelly six times larger than the variance for $\\texttt{females}$ under $\\texttt{mod0}$, there is no sufficient statistical evidence that they are different. Probably, this result is due to the limited number of parents and crosses considered in this study (_i.e._, low statistical power).

And as a final extension, for the constrained analysis we also can calculate a narrow-sense heritability using the function $\\texttt{vpredict}$, as:

```plaintext
vpredict(modC, h2~4*V2/(V1+V2+V3+V4))
```

|   | Estimate | SE |
| --- | --- | --- |
| h2 | 0.5543991 | 0.2486712 |

The above output provides us with an estimate of $h^2$, 0.554 (SE = 0.249), a reasonably high heritability estimate but with a large amount of uncertainty.

The appeal of this ASReml-R approach for constraining variance components is that, for any complex linear mixed models, constraints can be easily specified in column $\\texttt{V1}$ of matrix $\\texttt{M}$, for most hypothesis of interest, or to evaluate models with specific biological or mathematical constraints.

## Author

Salvador A. Gezan

## References

Becker (1984). _Manual of Quantitative Genetics_, 4th Edition. Academic Enterprises, Pullman, WA.

## Files to download

[White pine.txt](https://vsni1.wpengine.com/wp-content/uploads/2020/05/WHITEPINE.txt)

[Constrain.R](https://vsni1.wpengine.com/wp-content/uploads/2020/05/CONSTRAIN_R-FILE.txt)

**Notes**: SAG May-2020

ASReml-R Snippets: Constraining Variance Components

large_Pinecones-longest-sugar-pine-length-species.jpg

medium_Pinecones-longest-sugar-pine-length-species.jpg

small_Pinecones-longest-sugar-pine-length-species.jpg

thumbnail_Pinecones-longest-sugar-pine-length-species.jpg

KWS is one of the world’s leading suppliers of seeds to the farming industry and it is therefore no surprise that research into plant breeding and seeds is a crucial part of their work.

## KWS value precision and security

For centuries people have been breeding plants, trying to fine-tune individual plants to obtain or change a particular trait. Today plant breeders look at addressing concerns such as drought or environment tolerance, pesticide, herbicide, fungi or bacteria tolerance, improving nutritional value or increasing yields.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/KWS_genstat_M_Sanova_AS_Reml_997e9254db.jpg)

This, by definition, ranges from “classical” propagation techniques to the modern methods of molecular breeding and genetical modification. However, one thing is very apparent: the scientists engaged in any form of plant breeding need to be sure their data analysis is based on sound, solid and secure statistics.

KWS is one of the world’s leading suppliers of seeds to the farming industry and it is therefore no surprise that research into plant breeding and seeds is a crucial part of their work. The research and development teams at KWS focus on yield, seed quality, resistance to disease and pests, enhancing the nutrient quality and improved processing capability of the plants.

Recently the researchers at KWS decided to change their data analysis software portfolio by introducing Genstat. During the last 2 years main emphasis was put on establishing Genstat for statistical analysis of phenomic data from their field.

Genstat was chosen, because it was possible to analyse data from trials with large Alpha Designs with the REML techniques. Plus, Genstat’s REML algorithm dealt with missing values in a more superior way than the previous software. REML itself can be used to analyse models with several types of error variation (multi-level models) and to fit models to correlated data like repeated measurements. Genstat’s powerful command language streamlined the data analyses, reducing the calculation time by at least 70%.

Additionally Genstat’s efficient data handling meant that far less data pre-processing was needed. Genstat allows the full integration into automated analysis pipelines, thus decreasing calculation times by 30%. This allows the flexible addition of new statistical methods into the pipelines.

KWS researchers use a wide variety of different factors in their trials, using Block Designs, Alpha Designs, Lattice, Split Plot, Multi factorial Experiments. Genstat allows the researchers to easily design and analyse these and other types of experiments. Genstat also effortlessly handles large and complex data sets, which will be essential for the analysis of genomic data and their association with phenomic informations. Combine these factors with Genstat’s ability to easily export and import data from other databases and the plant breeder has very powerful analysis tools.

As well as the extensive technical capabilities of Genstat, KWS chose Genstat because of its history within agricultural science and research; it was originally designed by statisticians working on agricultural research at Rothamsted Research; where Fisher, Yates and Nelder (to mention but 3) developed statistical techniques that are central to modern statistics. Genstat’s connection with Rothamsted continues today. Developments in each new version of Genstat reflect the needs of the agricultural scientist, including ANOVA, design of experiments a host of multivariate analysis techniques and linear mixed model analysis. A huge variety and complexity of data analysis is possible in Genstat, and an understanding of the requirements of experiments in this area is reflected in the terminology and thinking behind Genstat. But Genstat’s history not only shows its suitability to any form of agricultural research but also shows the stability of its performance; it has been tried and tested for over 30 years.

Coupled with the power and background of Genstat, VSNi were able to provide a bespoke consultancy service to KWS to assist with their specific needs; this is possible because our developers have been working with agricultural researchers and statistics for years. The VSNi statisticians were able to talk to the researchers at KWS, in their language to find more efficient ways to run their analyses.

In a subject area that demands precision and security Genstat provides the data analysis solution that removes the unknowns surrounding the research. You can trust the statistics within Genstat, because it is developed by people who know and understand the issues in agricultural research today.

large_KWS_wheatfield_final.jpg

medium_KWS_wheatfield_final.jpg

small_KWS_wheatfield_final.jpg

thumbnail_KWS_wheatfield_final.jpg

Everything we do reflects our vision of ‘Safer food for the nation’. We aim to ensure that food produced or sold in the UK is safe to eat, consumers have the information they need to make informed choices about where and what they eat and that regulation and enforcement is risk-based and focused on improving public health.

![alt text](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/foodstandards_genstat_88ed61ea86.jpg)

Genstat has been deployed in many important areas of FSA activity, including the following:-

1\. Exposure Assessment

There are many types of chemical whose presence in food can pose some level of risk for consumers. These include environmental contaminants (e.g. lead or cadmium), as well as certain pesticides and food additives. The FSA needs to estimate the likely level of human exposure to such chemicals, either from individual foods, or from the diet as a whole. Our exposure assessment models combine occurrence data (on the level of a given contaminant in foods) with extensive human intake data (from individual food diaries). Genstat provides excellent tools for the necessary data import, indexing, manipulation and analysis and produces fast, reliable results. It has the flexibility to perform more complex modelling when required: e.g.  probabilistic assessments of risk that allow for variability and uncertainty in some of the key input variables.

2\. Risk Assessment

The FSA uses Hierarchical Generalized Linear Models to model the probability of elevated levels of naturally occurring biotoxins in shellfish. These estimates are used to prioritise biotoxin sampling activity around the shores of the UK, indicating where, and at what times of year, the level of sampling should be raised, and also where and when it can be reduced. Genstat has excellent procedures for performing models of this kind and for extracting model coefficients in a form that can be used to make a range of policy-relevant inferences.

large_rhubarb_final.jpg

medium_rhubarb_final.jpg

small_rhubarb_final.jpg

thumbnail_rhubarb_final.jpg

Food Standards Agency

Statistical inference is the process of drawing some conclusions about a population based on the sample data at hand. This is what we do when we fit a linear mixed model to our data, where some population parameters (_e.g._, treatment means) are estimated and used to draw conclusions about our study. Under statistical inference we want to expand our findings to the complete population from where our data originated.

The three most common types of statistical inference are: 1) point estimation, 2) interval estimation, and 3) hypothesis testing. Point estimation is where we estimate an unknown parameter that is presented as a single numeric value. In contrast, under interval estimation, we represent the estimated unknown parameter using a range of numeric values that is highly likely to contain the true population value. This range is known as the confidence interval. Finally, under hypothesis testing, we assess the evidence provided by the sample data in favor, or against, a given claim about the population. This is usually done using a statistical test.

All of the above types of statistical inference are commonly used for linear mixed models (LMM), and these can be implemented on fixed or random effects, functions of these effects (_e.g._, predicted means), and variance components or their functions (e.g., heritability).

In this tutorial we will present an example where statistical inference is done for a variance component. In this example, the aim of our LMM is to estimate a correlation variance component ($\\rho$), and we are interested in: 1) knowing what is the value of this correlation (point estimate), 2) reporting a 95% confidence interval (interval estimation), and 3) evaluating if this parameter is significantly different from zero (hypothesis testing).

The dataset used originates from a study that evaluated IQ scores for twins (McClave and Sincich, 2009). In this research, IQ scores of 32 pairs of identical twins were measured, where one of the twins ($\\texttt{Twin A}$) was reared by their natural parent, and the other ($\\texttt{Twin A}$) was reared by a relative or another person. Studies with twins are ideal because they share identical genetics, and therefore, their differences will be due exclusively to non-genetic factors, like in this case the type of rearing. For this study, the main objective is to evaluate if there are statistically significant differences between the average IQ scores between these two rearing levels.

The dataset for this example can be found at the end of this document in the file $\\texttt{TWINS.zip}. The first few lines are:

| PairID | TwinR | IQ |
| --- | --- | --- |
| 112 | A | 113 |
| 112 | B | 109 |
| 114 | A | 94 |
| 114 | B | 100 |
| 126 | A | 99 |
| 126 | B | 86 |

Here, we have a factor that identifies the pair of twins $\\texttt{PairID}$, and another that specifies the rearing type $\\texttt{TwinR}$($\\texttt{A}$ = reared by neutral person, $\\texttt{B}$ = reared by relative or other person). The raw means (and standard deviations) for groups $\\texttt{A}$ and $\\texttt{B}$ are 93.22 (15.17) and 96.13 (13.66), respectively. These values indicate small differences in IQ between these groups. In addition, we can also plot the IQ scores for each pair of twins, as shown in the below. We note a reasonable relationship between the IQ scores from the same twin pair, with a Pearson’s correlation of 0.816.

![Plot of IQ scores for pairs of twins](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/case_study_likelihood_1_5fd1cdef73.png)

Now, let’s fit our LMM. First, we will define our model in a generic way:  
$$  
y\_{ij} = \\mu + TwinR\_i + \\epsilon\_{ij}  
$$

where $y\_{ij}$ is the IQ measurement for the ijth twin; $\\texttt{TwinR}$ is the fixed effect of rearing type ($\\texttt{A}$ or $\\texttt{B}$), and $\\epsilon\_{ij}$ is the residual error, with $\\epsilon\_{ij} \\sim N(0, \\sigma^2)$.

For analysis, the above model needs to indicate that the errors, $\\epsilon\_{ij}$, are not independent from each other. In this case we need to specify that observations belonging to the same twin pair are correlated. This will be done in ASReml-R by a using the variance component structure $\\texttt{corv}$. For a single twin pair this structure has the following form:  
$$  
\\begin{bmatrix} \\sigma^2 & \\rho \\sigma^2 \\ \\rho \\sigma^2 & \\sigma^2 \\end{bmatrix} = \\sigma^2 \\begin{bmatrix} 1 & \\rho \\ \\rho & 1 \\end{bmatrix}  
$$

This variance structure is described by 2 variance components, one associated with the error variance, $\\sigma^2$, and the other, $\\rho$, representing the correlation between residual errors of the same twin pair. It is important to note that these parameters have natural constraints: $\\sigma^2$ can only be positive, and $\\rho$ has to be between -1 and +1.

In ASReml-R we can fit this model and produce basic output using the following code:

```plaintext
modT&lt;-asreml(fixed=IQ~TwinR,
             residual=~id(PairID):corv(TwinR),
             data=TWINS)
predict(modT,classify='TwinR')$pvals
wald.asreml(modT,denDF='numeric')$Wald
summary(modT)$varcomp
```

The above model has the fixed effect of $\\texttt{TwinR}$, but under $\\texttt{residual}$ we have defined a matrix $\\texttt{R}$ of variance-covariance between all residuals. Here, we have an identity matrix of dimension 32 × 32 ($\\texttt{id}$) multiplied (using the Kronoecker product, $\\otimes$) by a 2 × 2 matrix ($\\texttt{corv}$) defined as above. Therefore, this forms a $\\texttt{R}$ matrix of dimension 64 × 64 corresponding to the number of observations in the dataset, as shown below.  
$$  
\\mathbf{R}= \\sigma^2 \\begin{bmatrix} 1 & \\rho & 0 & 0 & \\cdots & 0 & 0 \\\\ \\rho & 1 & 0 & 0 & \\cdots & 0 & 0 \\\\ 0 & 0 & 1 & \\rho & \\cdots & 0 & 0 \\\\ 0 & 0 & \\rho & 1 & \\cdots & 0 & 0 \\\\ \\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\ 0 & 0 & 0 & 0 & \\cdots & 1 & \\rho \\\\ 0 & 0 & 0 & 0 & \\cdots & \\rho & 1 \\\\ \\end{bmatrix}  
$$

After fitting this model, we obtain a log-likelihood value of -183.650, and the following means with their SEMs (standard error of the mean) for each of the levels of $\\texttt{TwinR}$:

| TwinR | predicted.value | std.error | status |
| --- | --- | --- | --- |
| A | 93.21875 | 2.568317 | Estimable |
| B | 96.12500 | 2.568317 | Estimable |

The p-value of 0.074 ($F\_{1,31}$ = 3.42) from ANOVA indicates that there is no significant difference (at a 5% level) between the mean IQ of the two rearing levels. And finally, we can also report the estimated variance components:

|   | component | std.error | z.ratio | bound | %ch |
| --- | --- | --- | --- | --- | --- |
| PairID:TwinR!R | 1.0000000 | NA | NA | F | 0.0 |
| PairID:TwinR!TwinR!cor | 0.8125831 | 0.0610009 | 13.320844 | U | 0.0 |
| PairID:TwinR!TwinR!var | 211.0800432 | 48.7987205 | 4.325524 | P | 0.1 |

In relation to our population parameter of interest, $\\rho$, we have a point estimate of 0.813 with a standard error of 0.061. These are both REML (restricted maximum likelihood) variance component estimates, and they will have an asymptotic Normal distribution. For interval estimation, we can calculate an approximated 95% confidence interval (CI) for this population parameter using:  
$$

\\hat{\\rho} \\pm z\_{0.975} \\times SE(\\hat{\\rho}) = 0.813 \\pm 1.96 \\times 0.061 = \[0.693; 0.932\]

$$

In the above calculation we have used a z-value (assuming a Normal distribution), and our 95% CI for $\\hat{\\rho}$ is \[0.693; 0.932\]. This CI is symmetric and centered on the point estimate of 0.813, and it does not contain the value zero, indicating that zero might not be a likely value for this population parameter. In order to formally evaluate this, we define the following null and alternative hypotheses:

$$

H\_0: \\rho = 0

$$

$$

H\_1: \\rho \\ne 0

$$

This is a two-sided hypothesis, and we can use $\\hat{\\rho}$ and its SE obtained from the table of variance components presented earlier, and perform a Wald Test, of the form:

$$

(\\hat{\\rho} – 0)/{SE}(\\hat{\\rho}) = 0.813 / 0.061 = 13.321

$$

The above is a z-test and results in a p-value of &lt; 0.0001, implying that the null hypothesis should be rejected. This _Wald Test_ is the most common way to evaluate parameters that are asymptotically normal, something that occurs with large number of observations. However, in this case, the correlation parameter, $\\rho$, is been estimated with a small number of points: only 32 pairs, resulting in a liberal and approximated CI.

Most specialized literature recommends the use of likelihood-based statistical approaches, specifically the _likelihood ratio test_ (LRT) to statistically compare contrasting models. In order to evaluate the significance of a parameter, the LRT compares two models: one with the variance component estimated from the data, and another where the variance component is fixed to the value specified under the null hypothesis. In our example, we will evaluate the _full model_ (_i.e._, the one with the $\\rho$ component estimated), against a _restricted model_ (_i.e._, the one that has $\\rho = 0$). Our full model is the one presented earlier ($\\texttt{modT}$), and we will need to represent our restricted model based on the following variance structure:

$$

\\begin{bmatrix} \\sigma^2 & 0 \\\\ 0 & \\sigma^2 \\end{bmatrix} = \\sigma^2 \\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix}

$$

In this structure the correlation is assumed tobe exactly zero; hence, we are assuming that all residuals are independent of one another. This model can be fitted with the following ASReml-R code that replaces $\\texttt{cov}$ by $\\texttt{idv}$:

```plaintext
modT0&lt;-asreml(fixed=IQ~TwinR,
             residual=~id(PairID):idv(TwinR),
             data=TWINS)
```

The LRT evaluates our hypothesis by comparing the log-likelihood values and using the formula

$$

d = 2 \\times (logL\_{restr} – logL\_{full})

$$

where $logL\_{restr}$ and $logL\_{full}$ are the log-likelihood values for the restricted and full models, respectively. The $d$ statistic follows an approximated $\\chi^2$ distribution with degrees of freedom equal to the difference in the number of estimated variance components between both models (in this example 1 degree of freedom).

In addition, we can also calculate some _goodness-of-fit statistics_ to compare these models. The most common ones used for LMMs are the Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC):

$$

AIC = -2 \\times logL + 2 \\times t

$$

$$

BIC = -2 \\times logL + 2 \\times t \\times log(v)  
$$

where $t$ is the number of variance components in the model, and $v$ is the residual degrees of freedom. Both AIC and BIC are reasonable goodness-of-fit statistics for comparing nested and non-nested models. However, the LRT can only be used to statistically compare nested models.

The log-likelihood, AIC and BIC values for our two models ($\\texttt{modT}$ and $\\texttt{modT0}$) are presented in the following table.

| Model | logL | AIC | BIC |
| --- | --- | --- | --- |
| Restricted | \-200.385 | 402.77 | 404.897 |
| Full | \-183.65 | 371.301 | 375.555 |

Based on the above figures, we find that the value $d$ is 33.470. For a 5% significance level, the critical value for the $\\chi^2$ distribution with 1 degree of freedom is 3.841. Hence, we find that our full and restricted models are statistically different. This leads us to reject the null hypothesis that the correlation is zero. We can also compare the models using the AIC and BIC values, where the smaller the value the better the model. In this case, both criteria point to the full model as the better model.

We can implement the LRT within ASReml-R by using the following command:

```plaintext
lrt.asreml(modT0, modT, boundary=FALSE),
```

That has the option $\\texttt{boundary=FALSE}$ indicating that this is a two-sided test, and the output from this command is:

|   | df | LR-statistic | Pr(Chisq) |
| --- | --- | --- | --- |
| modT/modT0 | 1 | 33.46966 | 0 |

As you can see, the $d$ value of 33.47 is associated with a very small p-value of &lt; 0.0001. Hence, in summary, we have strong evidence to conclude that our population $\\rho$ parameter is very likely to be non-zero.

Recall that above we calculated a 95% CI for our correlation parameter based on the reported information from the variance components table, resulting in the interval \[0.693; 0.932\]. However, as with our Wald Test, this CI is also based on an asymptotic Normal distribution, which is a strong assumption given the low number of observations. An added difficulty with the approximated CIs is that in some cases they can have limits outside of the natural range of the true parameter, which in our case is \[-1; 1\].

Alternatively, we can calculate a _likelihood-based_ confidence interval. In order to do this, we first need to fit a model with the parameter of interest fixed to a range of values. Then, we observe the log-likelihood values for all of these models and identify, based on the LRT, those that will not be significantly different from our full model.

In our example, we fit our LMM hundreds of times where the correlation parameter is fixed to values ranging from 0.01 to 0.99 in steps of 0.001. To start the process using ASReml-R, we first fit the base model as:

```plaintext
mod.LR&lt;-asreml(fixed=IQ~TwinR,
             residual=~id(PairID):corv(TwinR),
             start.values=TRUE,
             data=TWINS)
```

In the above code, we are not actually fitting the model. Rather, by using the option $\\texttt{start.values=TRUE}$ we are requesting the setup table of the variance components. This table can be seen using:

```plaintext
vc.table &lt;- mod.LR$vparameters.table
```

| Component | Value | Constraint |
| --- | --- | --- |
| PairID:TwinR!R | 1.0000 | F |
| PairID:TwinR!TwinR!cor | 0.1000 | U |
| PairID:TwinR!TwinR!var | 105.5401 | P |

Here, we have the same variance components reported from the $\\texttt{summary(modT\_\\\\$varcomp)}$ command. However, we also have their starting values and if they have a constraint. These constraints are $\\texttt{F}$ (fixed) for a base variance component, $\\texttt{U}$ (unconstrained) for the correlation component and $\\texttt{P}$ (positive) for the variance components.

We can change, for example, our $\\rho$ parameter’s starting value to 0.903 and its constraint to $\\texttt{F}$ using

```plaintext
vc.table[2,2] &lt;- 0.903
vc.table[2,3] &lt;- 'F'
```

he modified $\\texttt{vc.table}$ is incorporated within a model fit using:

```plaintext
mod.LR&lt;-asreml(fixed=IQ~TwinR,
             residual=~id(PairID):corv(TwinR),
             G.param=vc.table,
             R.param=vc.table,
             data=TWINS)
```

This is the same as our original full model, but we have added the lines $\\texttt{G.param=vc.table}$ and $\\texttt{R.param=vc.table}$. After we fit the above model, we obtain the following variance component estimates:

|   | component | std.error | z.ratio | bound | %ch |
| --- | --- | --- | --- | --- | --- |
| PairID:TwinR!R | 1.000 | NA | NA | F | 0 |
| PairID:TwinR!TwinR!cor | 0.903 | NA | NA | F | 0 |
| PairID:TwinR!TwinR!var | 304.443 | 54.67787 | 5.567938 | P | 0 |

Notice from the above table that the correlation parameter $\\rho$ has been fixed at the value of 0.903. Also, this model results in a log-likelihood of -185.550, which is not as favorable as the one reported under the _unconstrained_ model $\\texttt{modT}$ of -183.650. (Note: We want to maximize the log-likelihood and therefore a higher value is better.)

Finally, we let’s fit the above model over the complete range of $\\rho$ values (from 0.01 to 0.99). For our Twin dataset, the plot of the correlation parameter against its log-likelihood value is presented below.

![Plot of the correlation parameter against its log-likelihood value ](https://web-global-media-storage-production.s3.eu-west-2.amazonaws.com/case_study_likelihood_2_c4004c332f.png)

In the above plot the maximum log-likelihood value is associated with a correlation estimate of 0.816 (as given in $\\texttt{modT}$). The red line represents the threshold line, below which we have significant differences compared to the full (unconstrained) model. This threshold line is calculated by:

```plaintext
loglik.thr &lt;- modT$loglik - chi.crit/2
```

where $\\texttt{modTloglik}$ is the log-likelihood value of the full model, and $\\texttt{chi.crit}$ corresponds to the critical $\\chi^2$ value for 1 df at the 5% significance level (= 3.841). Hence, any model that has a log-likelihood value smaller than -185.571 will be significantly different from our full model.

he dotted black vertical lines correspond to the places where the curve is intersected by our red threshold line. In this case, the fixed correlation values are 0.652 and 0.903, corresponding to our limits. Therefore, we can conclude that the 95% likelihood-based confidence interval for $\\rho$ is \[0.652; 0.903\]. This approach provides us with a non-symmetric CI, that also restricts the parameter of interest to its natural range.

In summary, we have performed different types of statistical inference in the context of linear mixed models, particularly associated with their variance component estimates. The use of likelihood-based inferences in hypothesis testing and confidence intervals is always recommended, as it provides more accurate inferences, and it has been reported that LRT has higher power than other statistical tests. One disadvantage of this inferential approach is the need to fit several, even hundreds, of models, in contrast with the use of asymptotic distribution assumptions that only requires a single model fit.

In our full model ($\\texttt{modT}$) we have used the $\\texttt{corv}$ variance structure, but we could have used more complex and flexible variance structures, such as $\\texttt{corh}$. For our data, this structure will have the following form:

$$

\\begin{bmatrix}  

\\sigma\_A^2 & \\rho \\sigma\_A \\sigma\_B \\\\  

\\rho \\sigma\_A \\sigma\_B & \\sigma\_B^2  

\\end{bmatrix}

$$

This variance structure has extended $\\texttt{corv}$ from a single residual variance to specific variances to each of the levels of rearing. This might be a reasonable assumption, and it could be evaluated and tested using LRT by comparing these models. We leave this to the reader to explore, but this is an extension of many modelling options of linear mixed models under ASReml-R.

## Author

Salvador A. Gezan

## Reference

McClave, JT; Sincich, T. 2009. _Statistics_, 12th Edition. Prentice Hall. Upper Saddle River, NJ.

## Files to download

[TWINS.zip](https://vsni1.wpengine.com/wp-content/uploads/2020/03/TWINS.zip)

**Notes:** SAG Feb-2020

Likelihood-based statistical inference for variance components in linear mixed models

Software

Sector

Learning

Case studies

About us

Contact

A world leader in the advancement and application of algorithmic and analytical content for the smart/precision biotech sector

Visit our other sites

h3analytics.com

biosci.global

vsnc.com.cn

Follow us

Copyright © 2000-2021 VSN International Ltd. | Privacy Policy | End-User License Agreement