Linear mixed models: obtaining estimates of random effects with ASReml-R

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

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

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 blocks and a total of genotypes evaluated:

where is the phenotypic observation of the ith genotype on the jth block; is the overall mean; is the fixed effect of the ith block; is the random effect of the ith genotype, with that is associated with some information about pedigree, from which we can obtain a pedigree-based or genomic-based relationship matrix ; and is the random residual with .

In the above model, also known as Animal Model, the 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  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, , also known as accuracy. This can be calculated using the following expression for a given genotype as:

This formula requires the PEV (predictor error variance) of , which corresponds to the square of the standard error (SE) of the random effect estimate, i.e. . Most statistical software will report these SE values, including SAS and ASReml. In addition, we have , 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 ().

According to statistical theory, the values of PEV will range from zero to the . 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,   . However, when there is very little information, then the PEV will approximate to . Note that these extreme values of PEV translate in 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.

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:

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

Here, the response variable was and our single random effect is the factor that represents our additive effects, and this is associated with the inverse of the relationship matrix . In addition, in the above code it was requested that observations with missing data (i.e., ) 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

vc <- summary(modelGBLUP)$varcomp

will output the following estimation of variance components.

 componentstd.errorz.ratiobound%ch
vm(lines, ainv)34.941597.8801384.434134P0.4
units!R43.100007.9493055.421857P0.0

These values correspond to a narrow-sense heritability of for the adjusted mean of height. Hence, the trait 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:

BLUP <- 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.

 solutionstd.errorz.ratioPEVaccuracyreliability
vm(lines, ainv)_LENTA5.83473.94821.477815.58850.74420.5539
vm(lines, ainv)_PALLAS5.11383.98071.284615.84620.73930.5465
vm(lines, ainv)_CHARIOT1.14943.98590.288415.88710.73850.5453
vm(lines, ainv)_CLASS-7.64293.9874-1.916815.89900.73820.5450
vm(lines, ainv)_BONUS7.92053.98971.985315.91760.73790.5445
vm(lines, ainv)_TOKEN-7.85753.9930-1.967815.94400.73740.5437
vm(lines, ainv)_MAJA9.92814.00222.480616.01780.73590.5416
vm(lines, ainv)_PRESTIGE-8.42054.0184-2.095516.14770.73340.5379
vm(lines, ainv)_ASPEN-4.14114.0456-1.023616.36700.72910.5316
vm(lines, ainv)_CORNICHE-1.61224.0546-0.397616.43940.72770.5295
 solutionstd.errorz.ratioPEVaccuracyreliability
vm(lines, ainv)_FREJA28.01774.82605.805623.29020.57750.3335
vm(lines, ainv)_DASIO-8.59694.8376-1.777123.40280.57470.3302
vm(lines, ainv)_KESTREL-0.43034.8644-0.088523.66230.56820.3228
vm(lines, ainv)_GAIRDNER5.31704.91881.081024.19440.55460.3076
vm(lines, ainv)_STELLA15.14224.92733.073124.27810.55240.3052
vm(lines, ainv)_LOGAN5.76014.97161.158624.71700.54090.2926
vm(lines, ainv)_CLARA9.48745.38701.761229.01970.41170.1695
vm(lines, ainv)_MARESI1.54495.74480.268933.00320.23550.0555
vm(lines, ainv)_FLICK SEJET-8.54435.8523-1.460034.24900.14080.0198
vm(lines, ainv)_BEKA16.59086.09922.720237.2001NaN-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 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 cm. The full range of EBVs is between − cm to cm, a wide interval considering that the mean of the population is 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 to . The width of this range is related to the level of genetic relationships between individuals, that is specified through the matrix . 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

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.

 PEVaccuracyreliability
vm(lines, ainv)_LENTA15.58850.74420.5539
vm(lines, ainv)_PALLAS15.84620.73930.5465
vm(lines, ainv)_CHARIOT15.88710.73850.5453
vm(lines, ainv)_CLASS15.89900.73820.5450
vm(lines, ainv)_BONUS15.91760.73790.5445
vm(lines, ainv)_TOKEN15.94400.73740.5437
vm(lines, ainv)_MAJA16.01780.73590.5416
vm(lines, ainv)_PRESTIGE16.14770.73340.5379
vm(lines, ainv)_ASPEN16.36700.72910.5316
vm(lines, ainv)_CORNICHE16.43940.72770.5295
 PEVaccuracyreliability
vm(lines, ainv)_FREJA23.29020.57750.3335
vm(lines, ainv)_DASIO23.40280.57470.3302
vm(lines, ainv)_KESTREL23.66230.56820.3228
vm(lines, ainv)_GAIRDNER24.19440.55460.3076
vm(lines, ainv)_STELLA24.27810.55240.3052
vm(lines, ainv)_LOGAN24.71700.54090.2926
vm(lines, ainv)_CLARA29.01970.41170.1695
vm(lines, ainv)_MARESI33.00320.23550.0555
vm(lines, ainv)_FLICK SEJET34.24900.14080.0198
vm(lines, ainv)_BEKA37.2001NaN-0.0646

For this trait, the range of correlation between true and predicted breeding values, , is from to , and the lowest accuracy and reliability values correspond to individuals with the largest SE values. Note that the genotype has a for accuracy, which originates from a negative square root, as this BLUP has an SE that is larger than .

alt text

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 values lower than . However, the majority of the observations are found with accuracy values of or higher, where some of them reach values as high as .

It is interesting to speculate what the reasons could be for these three lower observations (genotypes , and ) together with the observation (genotype ), 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 , such as , 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, 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 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

Related Reads