In animal breeding, it is essential to know genetic parameters such as heritability, with the aim of being able to predict genetic values (GV) and efficiently direct selection programs. A mixed model refers to those cases where the researcher considers fixed and random factors in a statistical model. Models widely used in the area of animal genetic improvement are the reproductive model and the animal model, which consider the reproductive or animal factor as random and a group of non-genetic effects as fixed. These mixed models allow us to obtain both heritability values (h²) for a trait, as well as genetic predictions such as the expected progeny difference (EPDs) or the predicted transmission ability (PTA) for each animal. An example of birth weight (BW) in cattle was used to calculate the VG, h2 and e2 using a mixed model, with a fixed and a random factor. The ANOVA, ML and REML methods were used to calculate h², e² and the VG first using all the information and subsequently assuming the last lost data, under a reproductive model and an animal model. The results found using the 3 methods were the same for REML and ANOVA in balanced data and different for the 3 methods in unbalanced data, where in the unbalanced case the ANOVA estimated a negative variance component, therefore, it can be concluded that estimate genetic values and parameters using ANOVA, ML and REML, but with the risk of estimating negative variance components using ANOVA or null (or overestimated) heritabilities with likelihood-based methods when the data structure or model is not the same correct.

Keywords: heritability, ANOVA, REML, ML, mixed models

In animal breeding, it is essential to know genetic parameters such as heritability in order to be able to predict genetic values (GV) and efficiently conduct selection programs. Genetic parameters are ratios between estimated population variances, known as variance components, which are calculated using linear models containing fixed and random factors, generally known as mixed models.¹ For the correct estimation of parameters and genetic values, it is necessary to have a broad knowledge of estimation using mixed models. Therefore, this article reviews the estimation of variance components and genetic values using ANOVA, ML and REML under a reproductive model and an animal model, explaining the virtues and limitations of each method in balanced and unbalanced data.

Theoretical framework mixed models

A mixed model refers to those cases where the researcher considers both fixed and random factors in a statistical model.² A model widely used in the area of animal breeding is the reproductive model or Sire Model, which considers the reproductive factor as random and a group of non-genetic effects as fixed.² The reproductive model allows obtaining both heritability values (h²) for a trait, as well as genetic predictions such as the expected difference of progeny (DEPs) or the predicted transmission ability (PTA) for each breeder.³ In matrix algebra the reproductive model takes the following form:

$y=Xb+Zs+e$

Where 𝑦 is a vector for the data, 𝑋 is an incidence matrix relating the data to the fixed effects, 𝑏 is a vector of unknown parameters for the fixed effects, 𝑍 is an incidence matrix relating the data to the random effects, 𝑠 is a vector of unknown predictions for each player, and 𝑒 is a vector of residuals.

The covariance structure of the above model is:

$VAR\left[\frac{s}{e}\right]=\left[\begin{array}{cc}I{\sigma }^2{}_s& 0\\ 0& I{\sigma }^2{}_s\end{array}\right]$

Where 𝐼 is an identity matrix, 𝜎_𝑠² is the variance between breeders and 𝜎_𝑒² is the residual variance.

The Henderson normal equations, necessary to find the genetic values of the breeders, for the above model are given by³:

$\left[\begin{array}{cc}X\text{'}X& X\text{'}Z\\ Z\text{'}X& Z\text{'}Z+I\alpha \end{array}\right]\left[\begin{array}{c}{b}_i\\ s_i\end{array}\right]=\left[\begin{array}{c}{X}^{\prime }y\\ Z\text{'}y\end{array}\right]$

Where 𝛼 is a ratio of the residual variance to the variance between breeders:

$\alpha =\frac{{\sigma }_e^2}{{\sigma }_s^2}$

According to Román and Aranguren,⁴ it is possible to substitute 𝐼𝛼 by 𝐴⁻¹𝛼 in the normal Henderson equations, with the objective of improving predictions using all the parentage information between males, therefore, the new equations are:

$\left[\begin{array}{cc}X\text{'}X& X\text{'}Z\\ Z\text{'}X& Z\text{'}Z+A^{-1}\alpha \end{array}\right]\left[\begin{array}{c}{b}_i\\ s_i\end{array}\right]=\left[\begin{array}{c}X\text{'}y\\ Z\text{'}y\end{array}\right]$

Where 𝐴⁻¹ is the inverse of the kinship matrix.

And the covariance structure taking into account the introduction of 𝐴 is⁵:

$VAR\left[\begin{array}{c}s\\ e\end{array}\right]-\left[\begin{array}{cc}A{\sigma }_s^2& 0\\ 0& I{\sigma }_e^2\end{array}\right]$

Another model widely used in genetic evaluation is the animal model, which uses all the parentage information in the pedigree, and unlike the reproductive model, allows obtaining genetic predictions of all the animals in the herd, whether or not data is present or not:

$y=Xb+Za+e$

Where 𝑎 is a vector of genetic predictions for each animal the covariance structure of the above model is as follows:

$VAR\left[\begin{array}{c}a\\ e\end{array}\right]=\left[\begin{array}{cc}A{\sigma }_a^2& 0\\ 0& I{\sigma }_e^2\end{array}\right]=\left[\begin{array}{cc}G& 0\\ 0& R\end{array}\right]$

Where 𝐺 is a variance and covariance matrix for the random effects and 𝑅 is a matrix of residuals.

The Henderson normal equations for this model are given by:⁶

$\left[\begin{array}{cc}X\prime X& X\prime {\rm Z}\\ {\rm Z}\prime X& {\rm Z}\prime {\rm Z}+{{\rm A}}^{-1}\alpha \end{array}\right]\left[\begin{array}{c}{b}_i\\ a_i\end{array}\right]=\left[\begin{array}{c}X\prime y\\ Z\prime y\end{array}\right]$ >

Where 𝛼 in this model is a ratio of the residual variance to the additive variance:

$\alpha =\frac{{\sigma }_e^2}{{\sigma }_a^2}$

Where 𝜎_𝑎² is additive genetic variance

Genetic parameters

Using mixed models, it is possible to estimate the variance components, and from them calculate the hereability, which is given by:⁵

$h^2=\frac{{\sigma }_a^2}{{\sigma }_p^2}$

Where ℎ² is heritability, and 𝜎_𝑝² is phenotypic variance, therefore, heritability is defined as a quotient between the additive variance and the phenotypic variance. The additive component (additive variance) of the numerator of the formula of ℎ² can be estimated using several procedures, a well-known one is to use a reproductive model to estimate the variance between breeders, which is ¼ of the additive variance, therefore, a formula to estimate 𝜎_𝑎² is:

${\sigma }_a^2=4{\sigma }_s^2$

Where 4𝜎_𝑠² is four times the variance among breeders, therefore, heritability can be calculated as:⁷

$h^2=\frac{4_s^2}{{\sigma }_p^2}$

If the heritability is known, the heritability component can be 𝜎_𝑎² component can be calculated using the following formula:⁸

${\sigma }_a^2=h^2{\sigma }_p^2$

Another parameter of interest is the environmental proportion coefficient, which indicates how much of the differences observed in the phenotype (data) of the animals are due to non-genetic (environmental) factors, this coefficient has the following mathematical formula:⁹

$e^2=\frac{{\sigma }_{en}^2}{{\sigma }_p^2}$

Where 𝜎_𝑒𝑛² is the environmental variance. The variance component 𝜎_𝑒𝑛² is calculated using the difference between 𝜎_𝑝² 𝜎_𝑎² therefore, the formula for 𝜎_𝑒𝑛² is:⁸

${\sigma }_{en}^2={\sigma }_p^2-{\sigma }_a^2$

Finally, the variance 𝜎_𝑝² is the sum of the variance components:

${\sigma }_p^2={\sigma }_a^2+{\sigma }_{en}^2={\sigma }_s^2+{\sigma }_e^2$

Variance component estimation using a reproductive model analysis of variance

There are several classical methods for estimating the variance components needed to compute ℎ² y 𝑒²including analysis of variance (ANOVA), maximum likelihood (ML) and restricted maximum likelihood (REML).

ANOVA is a technique that attempts to separate out different sources of variability. 𝜎_𝑝² into different sources of variability, this involves the separation of sums of squares (SC), degrees of freedom (GL) and mean squares (MS) for each source of variation. Variance components estimated using ANOVA are calculated by equating the expected values of the CM (E (CM)) for each source of variation, with their respective CM and solving the resulting system of equations.¹⁰ CMs are a ratio of SC to GLs for each source of variation:¹⁰

$CM=\frac{SC}{GL}$

In the case of a fixed factor and a random factor, without interaction, the reproductive model, in elementary algebra, is given by:

$y_{ijk}=\mu +s_i+b_j+e_{ijk}$

And the ANOVA square for the above model is presented in Table 1.

Where 𝑘 is the number of replicates of the design, 𝑛_{𝑓𝑖𝑗𝑜} is the number of levels of the fixed effect and 𝑛_𝑠 is the number of levels of the random factor. The variance components are calculated by equating the CM to their E (CM):

$CMs={\sigma }_e^2+n_{fijo}k{\sigma }_s^2$

$CMe={\sigma }_e^2$

And the unique solution of this system of equations is:

${\sigma }_s^2=\frac{CMs-CMe}{n_{fijo}k}$

${\sigma }_e^2=C{M}_e$

In balanced data, the CS can be estimated directly without the need for adjustment, for a model with two non-interacting factors:

$SCs=\frac{\sum y_{s.}^2}{k}-\frac{{\left(\sum y\right)}^2}{n}$

$SCb=\frac{\sum y_{b.}^2}{w}-\frac{{\left(\sum y\right)}^2}{n}$

Where ∑ 𝑦_𝑠²_. is the sum of the sum of the sum of the data for each player squared, ∑ 𝑦_𝑏²_. is the sum of the sum of the sum of the data for each level of the fixed effect and w is the number of replicates for the fixed effect. For the unbalanced case, the SCs have to be calculated using the type III SCs for the random factor (sire), since type III calculates the SCs of an effect by correcting them with respect to any other effect that does not contain it and orthogonal to any effect (if it exists) that contains it. Type III CS can be expressed as:¹¹

$SCs=SC\left(\mu ,s,b\right)-SC\left(\mu ,b\right)$

The 𝑆𝐶𝑠 is corrected for the effects of 𝜇 𝑦 𝑏where 𝜇 is the intercept or herd mean effect. In order to find the values of 𝑆𝐶𝑠 it is necessary to fit a complete model and calculate (𝜇, 𝑠, 𝑏, ) and subtract 𝑆𝐶(𝜇, 𝑏) a reduced model.

Maximum likelihood  

The maximum likelihood (ML) method is a classical method of parameter estimation proposed by Fisher,¹² but it was not until Hartley and Rao,¹³ that it was used for mixed models in general. Knowing the likelihood function as a function of the parameters of a statistical model given some data, in ML we try to obtain estimators of the variance components that maximize the likelihood function, that is, that have the maximum probability of representing the population parameters.

The likelihood function is defined as the product of the likelihood function of the data, but in practice, the natural logarithm of the likelihood function is used because it is more manageable, if the distribution of the data is normal, in matrix algebra the natural logarithm of the likelihood function is defined as:¹¹

$Ln\left(L\right)=-0.5\left(n\right).In\left(2\pi \right)-0.5In\left|V\right|-0.5\left(y-Xb\right)\prime V^{-1}\left(y-Xb\right)$

Where (𝐿) is the natural logarithm of the likelihood function and 𝑉 = 𝑍𝐺𝑍^′ + 𝑅 is the variance and phenotypic covariance matrix of the model. To find the estimators that maximize the likelihood, we need to find the maximum of equation (𝐿)This is achieved with different methodologies, for example, if the data structure is balanced and we have a mixed model, with a random effect and a fixed one with no interaction, the derivative of 𝐿𝑛(𝐿)with respect to the parameters to be estimated σ²_s y σ²_e will lead us to a system of equations whose solution is:

${\sigma }_s^2=\frac{\frac{SCs}{n_s}-{\sigma }_e^2}{n_{fijo}k}$

${\sigma }_e^2=\left[1-\frac{n_{fijo}-1}{n_s\left(n_{fijo}k-1\right)}\right]CMe$

An important point of ML estimation, for this model, is that even with balanced data, it is possible to find estimators different from the ones presented above, since these solutions will be valid if the inequality 𝐶𝑀𝑠 > 𝐶𝑀𝑒is met, but on the other hand, if the inequality is 𝐶𝑀𝑠 < 𝐶𝑀𝑒 ML estimates for this model and balanced data are given by:¹¹

${\sigma }_e^2=\frac{\stackrel{{\sigma }_s^2=0}{SCtotal}}{\left(n_s\right)\left(n_{fijo}\right)\left(k\right)}$

That is all phenotypic variability is residual, which may indicate that the model used is incorrect or that the number of data is insufficient, thus increasing the variability of the error. The variance σ²_p is the sum of the variance components σ²_e y σ²_s whose sum gives an estimate of σ²_p given mathematically by:

${\sigma }_p^2={\sigma }_s^2+{\sigma }_e^2=\frac{\left(y-Xb\right)\prime \left(y-Xb\right)}{n}$

Which is biased, since it is associated with n degrees of freedom. If the structure of the daros is unbalanced, the partial derivatives of (𝐿) lead to nonlinear maximum likelihood equations for the parameters to be estimated, therefore, the system of equations cannot be solved with direct methods. Faced with this problem, iterative number methods are used to try to approximate the maximum of (𝐿) which are applied to the logarithmic likelihood itself and not to the equations resulting from its first derivative, in order to be able to simultaneously calculate the variance components and (𝐿)which we can use to find fit criteria for our model, such as the Akaike information criterion (AIC) and the Bayesian information criterion (BIC).

Restricted maximum likelihood

The restricted maximum likelihood method (REML) is a method proposed by Paterson and Thompson,¹³ which takes into account the loss of degrees of freedom by including fixed effects in the statistical model, therefore, the estimation of variability components are unbiased, since they are associated to degrees of freedom, which leads to an estimation of variance of the model. 𝑛 − (𝑋) degrees of freedom, which leads to an estimate of the variance, which is defined as σ²_pwhich is defined as:

${\sigma }_p^2=\frac{\left(y-Xb\right)\prime \left(y-Xb\right)}{n-Rango\left(X\right)}$

Where (𝑋) is the rank of the incidence matrix for the fixed effects of the model. For the case where the only fixed effect is 𝜇the variance σ²_p is associated with 𝑛 − 1 degrees of freedom.

As in ML, in REML, the objective is to maximize the logarithm of a function of the parameters, but in this case restricted, which is known as restricted likelihood function, which in matrix algebra is defined as:¹⁴

$Ln\left(Lr\right)=-0.5\left(n-p\right).In\left(2\pi \right)-0.5In\left|V\right|-0.5\left|X\prime V^{-1}X\right|-0.5\left(y-Xb\right)\prime V^{-1}\left(y-Xb\right)$

Where (𝐿_𝑟) is the logarithm of the restricted likelihood function. If we have a balanced data structure and by deriving (𝐿_𝑟) as a function of the variability components of the model (model above), we can solve a system of equations that give rise to estimates given by:¹¹

${\sigma }_s^2=\frac{CMs-CMe}{n_{fijo}k}$

${\sigma }_e^2=CMe$

These are identical to estimates using an ANOVA, since a property of REML is that in a balanced data structure, REML estimates = ANOVA as long as the inequality is satisfied. 𝐶𝑀𝑠 > 𝐶𝑀𝑒Otherwise, estimates via ANOVA would be negative and in REML all phenotypic variability is residual. In unbalanced data structure, the derivative of (𝐿_𝑟) with respect to the variance components, gives rise to nonlinear equations, which cannot be solved directly, therefore, in these cases, as in ML, iterative numerical methods are used to approximate the value of the variance components.

REML estimates using kinship information in an animal model

In the case of a simple animal model, where each animal has only one data (and there are animals without data), the ANOVA method cannot be applied, since it is not possible to estimate the variation within groups using this methodology, because the classification variable is each animal that has a unique record, but the ML and REML estimations are applicable since they allow introducing kinship information in the matrix. 𝐴. In a mixed model, maximizing (𝐿_𝑟) is equivalent to minimize −2 (𝐿_𝑟) Therefore, the objective function to be minimized, in matrix algebra, can be defined as:¹⁴

$-2Ln\left(Lr\right)=\left(n-p\right).In\left(2\pi \right)+Ln\left|R\right|+Ln\left|G\right|+Ln\left|C\right|+y\prime Py$

Where 𝐿𝑛|𝐶| is the natural logarithm of the determinant of the coefficient matrix of the normal Henderson equations and ′𝑃𝑦 is the generalized residual sum of squares. Obviously to minimize −2𝐿𝑛(𝐿_𝑟) iterative numerical methods are needed, but it has the advantage that it is easier than maximizing 𝐿𝑛(𝐿_𝑟) Therefore, most specialized REML programs use sparse matrix algorithms and numerical methods to try to find estimators resulting from the minimization of −2𝐿𝑛(𝐿_𝑟).

An example of birth weight (BW) in cattle was used to calculate the VG, ℎ² y 𝑒² using a mixed model, with a fixed and a random factor. The database is presented in Table 2. In this problem, we want to eliminate the variability that exists between the sexes, therefore, the sex factor is considered as fixed and the father factor as random, which leads us to the statistical model for this problem:

$PN=media+padre+sexo+error$

ANOVA, ML, and REML methods were used to calculate ℎ² , 𝑒² and GVs using the data in Table 2, first using all the information and then assuming the last missing data. For the animal model, a similar model was used:

$PN=media+animal+sexo+error$

Where all the kinship information and the value of the variance components found in the previous model were used to solve the Henderson normal equations.

Balanced data in a reproductive model

To calculate the CM, it is necessary to calculate the SC and GL for each source of variation, for this model and our data structure, we can calculate them using the formulas in Table 1.

$\begin{array}{l}G{L}_s=3-1=2\\ G{L}_{sexo}=2-1=1\\ G{L}_{total}=12-1=11\\ G{L}_e=11-1-2=11-3=8\end{array}$

And since the design is balanced, the SCs are:

$\begin{array}{l}S{C}_{total}={36}^2+{35}^2+{33}^2+\dots +{35}^2-\frac{{395}^2}{12}=152.916\\ S{C}_s=\frac{{\left(132\right)}^2+{\left(124\right)}^2+{\left(139\right)}^2}{4}-\frac{{395}^2}{12}=28.166\\ S{C}_{sexo}=\frac{{\left(185\right)}^2+{\left(210\right)}^2}{6}-\frac{{395}^2}{12}=52.083\\ S{C}_e=152.916-\left(28.166+52.083\right)=72.667\end{array}$

And the CMs come from:

$\begin{array}{l}C{M}_S=\frac{28.166}{2}=14.083\\ C{M}_e=\frac{72.667}{8}=9.083\end{array}$

And from the CM we can calculate the variance components:

${\sigma }_s^2= \frac{14.083-9.083}{2\left(2\right)}=1.25$

Therefore, ℎ² using ANOVA is:

$\begin{array}{l}{h}^2=\frac{4\left(1.25\right)}{1.25+9.083}=0.483\\ e^2=\frac{10.33-5}{10.33}=0.515\end{array}$

And these ANOVA estimates, too, are REML, since the data structure is balanced and the 𝐶𝑀𝑠 > 𝐶𝑀𝑒. Now the calculation of the GVs, using the REML estimates, comes from the solutions of the normal Henderson equations, using the estimated value of the variance components, and calculating the value of 𝛼 we have that:

$\alpha =\frac{9.083}{1.25}=7.2664$

Therefore, the equations are:

$\left(\begin{array}{ccc}6& 0& \begin{array}{ccc}2& 2& 2\end{array}\\ 0& 6& \begin{array}{ccc}2& 2& 2\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}4+7.266& 0& 0\end{array}\\ \begin{array}{ccc}0& 4+7.266& 0\end{array}\\ \begin{array}{ccc}0& 0& 4+7.266\end{array}\end{array}\end{array}\right)\left[\begin{array}{c}\begin{array}{c}{b}_1\\ b_2\end{array}\\ s_1\\ s_2\\ s_3\end{array}\right]=\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]\to \left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2& 2& 2\end{array}\\ 0& 6& \begin{array}{ccc}2& 2& 2\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}11.266& 0& 0\end{array}\\ \begin{array}{ccc}0& 11.266& 0\end{array}\\ \begin{array}{ccc}0& 0& 11.266\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]=\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]$

In the previous equations the value of 𝜇 was forced to be zero in order to break the linear dependence between the rows and columns of the coefficient matrix. The solution of this system of equations is given by:

$\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]={\left[\left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2& 2& 2\end{array}\\ 0& 6& \begin{array}{ccc}2& 2& 2\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}11.266& 0& 0\end{array}\\ \begin{array}{ccc}0& 11.266& 0\end{array}\\ \begin{array}{ccc}0& 0& 11.266\end{array}\end{array}\end{array}\right]\right]}^{-1}\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]=\left[\begin{array}{c}30.83\\ 35\\ 0.029\\ -0.680\\ 0.650\end{array}\right]$

ML estimates are:

$\begin{array}{l}{\sigma }_s^2=\frac{\frac{28.166}{3}-8.073}{2\left(2\right)}=0.328\\ {\sigma }_e^2= \left[1-\frac{2-1}{3\left(2\left(2\right)-1\right)}\right]9.083=8.073\end{array}$

Y ℎ²y 𝑒²using ML is:

$\begin{array}{l}{h}^2=\frac{4\left(0.328\right)}{0.328+8.073}=0.156\\ e^2=\frac{8.401-1.3212}{8.401}=0.842\end{array}$

And the equations using ML are:

$\left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2& 2& 2\end{array}\\ 0& 6& \begin{array}{ccc}2& 2& 2\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}28.612& 0& 0\end{array}\\ \begin{array}{ccc}0& 28.612& 0\end{array}\\ \begin{array}{ccc}0& 0& 28.612\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]=\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]$

Therefore the solution is:

$\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]={\left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2& 2& 2\end{array}\\ 0& 6& \begin{array}{ccc}2& 2& 2\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}28.612& 0& 0\end{array}\\ \begin{array}{ccc}0& 28.612& 0\end{array}\\ \begin{array}{ccc}0& 0& 28.612\end{array}\end{array}\end{array}\right]}^{-1}\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]=\left[\begin{array}{c}30.83\\ 35\\ 0.011\\ -0.267\\ 0.256\end{array}\right]$

Introduction of the parentage matrix in a reproductive model

Now it is assumed that animal 1 is the father of animal 2, therefore, the equations take into account all the genealogy between males. First we have to calculate 𝐴⁻¹. Applying Henderson's rules,⁵ we have:

$A^{-1}=\left[\begin{array}{ccc}1+1/3& -2/3& 0\\ -2/3& 1/4& 0\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}1/4& -2/3& 0\\ -2/3& 1/4& 0\\ 0& 0& 1\end{array}\right]$

Therefore,

$A^{-1}\alpha =\left[\begin{array}{ccc}\frac{1}{4}& -\frac{2}{3}& 0\\ -\frac{2}{3}& \frac{1}{4}& 0\\ 0& 0& 1\end{array}\right]\left(7.2664\right)=\left[\begin{array}{ccc}1.8166& -4.8442& 0\\ -4.8442& 1.8166& 0\\ 0& 0& 7.2664\end{array}\right]$

And adding the Z'Z matrix:

$Z\prime Z+A^{-1}\alpha =\left[\begin{array}{ccc}4& 0& 0\\ 0& 4& 0\\ 0& 0& 4\end{array}\right]+\left[\begin{array}{ccc}1.8166& -4.8442& 0\\ -4.8442& 1.8166& 0\\ 0& 0& 7.2664\end{array}\right]=\left[\begin{array}{ccc}5.8166& -0.8442& 0\\ -0.8442& 5.8166& 0\\ 0& 0& 11.2664\end{array}\right]$

Therefore, the Henderson normal equations are:

$\left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2.0000& 2.0000& 2.0000\end{array}\\ 0& 6& \begin{array}{ccc}2.0000& 2.0000& 2.0000\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}5.8166& -0.8442& 0.0000\end{array}\\ \begin{array}{ccc}-.08442& 5.8166& 0.0000\end{array}\\ \begin{array}{ccc}0.0000& 0.0000& 11.2664\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]=\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]$

And the solution of these equations is:

$\left[\begin{array}{c}{b}_1\\ b_2\\ s_1\\ s_2\\ s_3\end{array}\right]={\left[\begin{array}{ccc}6& 0& \begin{array}{ccc}2.0000& 2.0000& 2.0000\end{array}\\ 0& 6& \begin{array}{ccc}2.0000& 2.0000& 2.0000\end{array}\\ \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}2\\ 2\\ 2\end{array}& \begin{array}{c}\begin{array}{ccc}5.8166& -0.8442& 0.0000\end{array}\\ \begin{array}{ccc}-0.8442& 5.8166& 0.0000\end{array}\\ \begin{array}{ccc}0.0000& 0.0000& 11.2664\end{array}\end{array}\end{array}\right]}^{-1}\left[\begin{array}{c}185\\ 210\\ 132\\ 124\\ 139\end{array}\right]=\left[\begin{array}{c}31.6285\\ 35.7952\\ -0.7765\\ -1.9776\\ 0.3685\end{array}\right]$

In this solution, the fixed effect (sex) is adjusted for the random example, and the random effect is adjusted for the fixed effect.

Unbalanced data in a reproductive model

Assuming the last missing data, the 𝑆𝐶𝑡𝑜𝑡𝑎𝑙 is:

$SCtotal=S{C}_{total}={36}^2+{35}^2+{33}^2+\dots +{29}^2-\frac{{360}^2}{11}=148.1818182$

The GLs are:

$\begin{array}{l}G{L}_s=3-1=2\\ G{L}_{sexo}=2-1=1\\ G{L}_{total}=11-1=10\\ G{L}_e=10-1-2=11-3=7\end{array}$

To calculate the ordinary least squares solutions for the reduced model (without the sire factor) must be estimated:

$b={\left[\begin{array}{ccc}11& 5& 6\\ 5& 5& 0\\ 6& 0& 6\end{array}\right]}^{-1}\left[\begin{array}{c}360\\ 150\\ 0\end{array}\right]=\left[\begin{array}{ccc}0.1666& -0.1666& 0\\ -0.1666& 0.3666& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}360\\ 150\\ 0\end{array}\right]=\left[\begin{array}{c}35\\ -5\\ 0\end{array}\right]$

Therefore, the 𝑆𝐶𝑠𝑒𝑥𝑜 for the reduced model is:

$S{C}_{sexo\left(reducido\right)}=\left[35-5 0\right]\left[\begin{array}{c}360\\ 150\\ 210\end{array}\right]-\frac{{360}^2}{11}=\left(35\right)\left(360\right)+\left(-5\right)\left(150\right)-\frac{{360}^2}{11}=68.1818$

To obtain the adjusted SCs, we calculate the SC for the full model and subtract from it, respectively (𝑟𝑒𝑑𝑢𝑐𝑖𝑑𝑜):

$SCs=SC\left(\mu ,s,sexo,\right)-SC\left(\mu ,sexo\right)=83.6818182-68.1818=15.5$

And the 𝑆𝐶_𝑒 is:

$S{C}_e=148.1818182-\left(15.5+68.1818\right)=64.5$

Therefore, the CMs are:

$\begin{array}{l}C{M}_s=\frac{15.5}{2}=7.75\\ C{M}_e=\frac{64.5}{7}=9.214286\end{array}$

Finally, the variance components are:

$\begin{array}{l}{\sigma }_s^2= \frac{7.75-9.214286}{3.6}=-0.4067\\ {\sigma }_e^2= 9.214286\end{array}$

The component σ²_s component is a negative estimate of the variance, because the 𝐶𝑀𝑠 < 𝐶𝑀𝑒is negative, therefore, the ML and REML estimates for σ²_s are:

σ²_s = 0

In other words, all the total variability is residual. Table 3 shows the ML and REML iterations for the calculation of the residual variance with the Newton-Rapson method: As shown in Table 3, when an unbalanced database is used, the estimates for the variance components are different in ANOVA and REML, since with ANOVA we have σ²_e = 9.2142 and with REML σ²_e = 8.88. For this particular case, ℎ² = 0 because the variance of the numerator is zero, obviously to find a more credible estimate, one should increase the number of data used in the genetic evaluation or try another model and compare the AIC.

Animal model

The solutions of the Henderson normal equations, using the values of σ²_e = 9.083 y σ²_a = 5 are presented in Table 4 Henderson's normal equations are not presented for this case due to its large dimensions.

Genetic values and parameters can be estimated using ANOVA, ML and REML, but with the risk of estimating negative variance components using ANOVA or zero (or overestimated) heritabilities with likelihood-based methods when the data structure or model is not correct. When the data structure is unbalanced, mathematical calculations with ANOVA, ML and REML are more complex and require computational algorithms with higher performance.

None.

The authors declared that there are no conflicts of interest.

1. Caballero J, Pablo E, Martinez C. Restricted maximum likelihood estimation of variance and covariance components of multiple traits under designs i and ii of North Carolina. Rev Fitotec Mex. 2003;26(1):53–66.
2. Castejón, Osiris. Design and analysis of experiments with statistix. Venezuela: Rafael Urdaneta University; 2008.
3. Elzo M, Garay O. Evolution of genetic improvement practices in domestic animal populations. Gainesville: University of Florida; 2012.
4. Román R, Aranguren A. Genetic evaluation of breeding stock: achievements and challenges From: Achievements and challenges in dual purpose cattle breeding; 2014.
5. Gutiérrez P. Introduction to animal genetic assessment methodology adapted to the EHEA. Madrid: Complutense; 2010.
6. Aranguren A, Román R. The simple animal model: a methodology for geneticists from: achievements and challenges of dual purpose cattle breeding; 2014.
7. Perez J, Montiel N. Heritability of the IBMI index of Italian buffaloes used in artificial insemination. Rev Cientif FCV-LUZ XXXIII. 2023;205–207.
8. Román R, Aranguren J, Garcidueñas R, et al. Association between reproductive characteristics and milk production in crossbred heifers. Revista Espamciencia. 2023;14(2):63–70.
9. Montgomery D. Experiments with random factors. From: Design and analysis of experiments. 2nd edn. Limusa-wiley; 2015.  
10. Searle R, Casella G, Mcculloch C. Variance Components. Wiley series; 1992.
11. Aldrich J. R. A. Fisher and the making of maximum likelihood 1912 - 1922. Statist Sci. 1997;12(3):162–176.
12. Hartley HO, Rao JNK. "Maximum likelihood estimation for the mixed analysis of variance model". Biometrika. 1967;54(1):93–108.
13. Patterson HD, Thompson R. "Recovery of inter-block information when block sizes are unequal". Biometrika. 1971;58:545–554.
14. Meyer K. Restricted maximum likelihood to estimate variance components for animal models with several random effects using a derivative-free algorithm. Genet Sel Evol. 1989;21:317–340.
15. Román R, Aranguren J, Villasmil Y, et al. Analysis of fertility at first service in dual-purpose heifers under an animal model. Revista Científica. 2010;20(4):383–389.