Structural Equation Modelling of Additive Genetic and Residual Covariance Matrices in Beef Cattle

1. Introduction

In animal breeding, standard multiple-trait mixed models (SMTM) are often applied to analyze multiple traits [1], estimating unstructured covariance matrices. The number of (co)variance parameters involved grows quadratically with the number of traits, and when phenotypes are highly correlated, the residual and genetic covariance matrices can be singular, making the convergence of estimation algorithms complex [1,2] and inferences unreliable.

Structural equation models (SEM) are often used to obtain more parsimonious (and sometimes more interpretable) parametrization of multivariate models. This approach has been used in quantitative genetics to structure the genetic and environmental covariance matrices of multivariate linear mixed models [2,3,4,5,6], commonly used for genetic evaluations. Primarily due to its ability to consider causal relationships between latent variables, allowing for the modeling of complex phenomena while simultaneously reducing the dimensionality of the data [7,8,9]. SEM methodology comprises many different techniques and procedures used together to model covariance structures. These models can be viewed as an extension of the multivariate linear mixed model [10] that can express functional networks among traits [2]. Gianola and Sorensen [3] discussed the use of recursive (REC) and simultaneous equation models (both special cases of SEM) acting on phenotypes. Valente et al. [11] proposed a methodology to search for recursive causal structures among phenotypes within a multivariate linear mixed model, which has been applied to various species and traits, including carcass and meat quality traits in beef cattle [12]. Some authors discussed the use of REC models operating on genetic covariance matrices [6,8,13]. In certain cases, modeling the covariance matrices helps infer or analyze causality between the breeding values of different traits or examine the ratio of two correlated traits, thereby reducing the model’s dimension with more parsimonious and efficient parameterizations [8,14].

Alternatively, factor analysis (FA), another specific case of SEM, can be used to reduce the dimension of the estimated genetic covariance matrix [15,16] and achieve a more parsimonious model of genetic effects within a multivariate linear mixed model, without decreasing the original number of traits and records. In the FA framework, the observed traits are modeled as linear combinations of fewer unobservable latent factors and trait-specific residual components, effectively capturing the common sources of variation underlying correlated traits while preserving individual trait identities. This approach has been successfully applied in quantitative genetics to structure genetic and environmental covariance matrices [16,17], proving particularly useful when analyzing highly correlated traits where traditional unstructured covariance matrices may suffer from convergence issues or poor estimability. The reduction in the number of parameters not only improves computational efficiency but can also lead to more stable parameter estimates and better model convergence, especially in datasets with limited sample sizes relative to the number of traits being analyzed [7]. Moreover, FA models can provide biological insights by identifying latent factors that may represent underlying biological processes or developmental pathways common to multiple phenotypes.

While previous studies have demonstrated the utility of SEM in various contexts, their application to structuring both genetic and residual covariance matrices independently in beef cattle breeding remains underexplored. In the current omics era, breeding programs are increasingly collecting larger numbers of correlated phenotypes, making computationally efficient, statistically robust methods for multi-trait genetic evaluation even more critical. These more parsimonious models offer practical advantages for routine genetic evaluations in breeding programs, enabling the simultaneous analysis of multiple traits across large populations with improved computational efficiency and numerical stability compared to fully parameterized models. Therefore, the objective of this study was to evaluate the performance of factor analysis and recursive models applied independently to additive genetic and residual covariance matrices for growth and ultrasound carcass traits in Nellore cattle, comparing them with standard multiple-trait mixed models in terms of parameter estimation, breeding value prediction, and model fit criteria.

2. Materials and Methods

2.1. Dataset Description

The present study used data from 2,942 Nellore animals collected over 3 years from an existing database of 10 herds participating in the Brazilian Nellore Breeding Program across 6 States. On-farm data collection was conducted according to routine management procedures and in accordance with the [18] guidelines. To collect carcass data without slaughtering the animals, the ultrasound methodology was applied [19]. Real-time ultrasound measurements of carcass traits included the longissimus muscle area (LMA, cm²) and backfat thickness (BF, mm), both obtained from a cross-sectional image of the longissimus dorsi muscle between the 12th and 13th ribs. Rump fat thickness (RF, mm) was measured at the intersection of the gluteus medius and biceps femoris muscles, located between the hooks and pin bones. Backfat thickness was measured at the 3/4 position from the chine bone end of the longissimus muscle, using the cross-sectional ribeye image. Other traits recorded included body weight (BW, kg) and hip height (HH, cm), both measured at the date of ultrasound scanning, and 450-day age-standardized scrotal circumference (SC, mm). Age-standardized SC was obtained as described in the Guidelines for Uniform Beef Improvement Programs [20]. The traits were assessed in animals aged between 480 and 629 days.

Contemporary groups (CG) were defined as animals of the same sex, born in the same herd, year, and season (dry or rainy), and reared within the same management group. A total of 302 contemporary groups were formed, with an average of 9.7 animals per group (range: 3 to 45 animals). Additional effects considered in the analyses included age of animal at scanning (linear covariate for all traits except SC) and age of dam at calving (linear and quadratic covariates for BF, RF, HH, and BW; classes for LMA and SC). A description of the traits and data structure for Brazilian Nellore cattle is presented in Table 1.

2.2. Standard Multi-Trait Mixed Model

Genetic parameters and model solutions were estimated using Bayesian inference implemented via Gibbs sampling. In the following sections, we first describe the standard multi-trait mixed model commonly used in genetic evaluation, then present extensions based on factor analytic (FA) and recursive (REC) model structures.

In the standard multi-trait mixed model (SMTM), a sire model was fitted for p traits measured on n individuals. The mixed-model equation can be described as:

$$\mathbf{y}=\mathbf{X}\beta +\mathbf{Z}\mathbf{u}+ϵ,$$

where $\mathbf{y}$ is the vector of phenotypic observations, $\beta$ is the vector of fixed effects, $\mathbf{u}$ is the vector of random additive sire effects, $ϵ$ is the vector of residual effects, and $\mathbf{X}$ and $\mathbf{Z}$ are the corresponding incidence matrices. Under the Bayesian framework, the likelihood function was specified as:

$$\mathbf{y}|\beta ,\mathbf{u},{\mathbf{G}}_0,{\mathbf{R}}_0\sim N(\mathbf{X}\beta +\mathbf{Z}\mathbf{u},{\mathbf{I}}_n\otimes {\mathbf{R}}_0),$$

with prior distributions $\mathbf{u}|\mathbf{A},{\mathbf{G}}_0\sim N(\mathbf{0},\mathbf{A}\otimes {\mathbf{G}}_0)$ and $ϵ|{\mathbf{R}}_0\sim N(\mathbf{0},{\mathbf{I}}_n\otimes {\mathbf{R}}_0)$, where ${\mathbf{G}}_0$ is the $p\times p$ additive genetic (co)variance matrix among sire effects, ${\mathbf{R}}_0$ is the $p\times p$ residual (co)variance matrix, $\mathbf{A}$ is Wright’s numerator relationship matrix among sires, ${\mathbf{I}}_n$ is an identity matrix of order n, and ⊗ denotes the Kronecker product operator.

Prior distributions for model parameters were specified as follows. Systematic effects included contemporary groups, animal age at scanning (linear effect, except for SC), and dam age (linear and quadratic effects for BF, RF, HH, and BW). For the vector of systematic effects $\beta$ with vague information [21], prior distributions were assigned as:

$$\beta \sim N(0,\mathbf{I}\otimes {100}^2),$$

For the (co)variance matrices, scaled inverse Wishart distributions were assigned as prior distributions:

$${\mathbf{G}}_0\sim \mathrm{IW}({\mathbf{S}}_g,{\nu }_g)\mathrm{and}{\mathbf{R}}_0\sim \mathrm{IW}({\mathbf{S}}_e,{\nu }_e),$$

where ${\mathbf{S}}_g$ and ${\mathbf{S}}_e$ are $p\times p$ scale matrices and ${\nu }_g$ and ${\nu }_e$ are degrees of freedom parameters. The joint posterior distribution of all model parameters is proportional to the product of the likelihood and prior distributions:

$$p(\beta ,\mathbf{u},{\mathbf{G}}_0,{\mathbf{R}}_0|\mathbf{y})\propto p\left(\mathbf{y}\right|\beta ,\mathbf{u},{\mathbf{G}}_0,{\mathbf{R}}_0)\times p\left(\mathbf{u}\right|\mathbf{A},{\mathbf{G}}_0)\times p\left(\beta \right)\times p\left({\mathbf{G}}_0\right)\times p\left({\mathbf{R}}_0\right).$$

2.3. Factor Analytic Model

Factor analysis is a multivariate statistical technique that simplifies complex and correlated data structures [22,23]. It can be viewed as an extension of principal component analysis, where the factor-analytic variance-covariance structure serves as a parsimonious approximation to a fully unstructured variance-covariance matrix, enabling more efficient parameterization of complex correlation patterns. In a standard FA, a vector of random variables (say $\mathbf{u}$) can be decomposed as a linear combination of fewer unobservable random variables called common factors ($\mathbf{f}$), with an unobservable incidence matrix ($\Lambda$) of factor loadings, plus a vector of trait-specific errors ($\delta$) to each individual i that represent the lack of fit of the model. In compact notation:

$${\mathbf{u}}_i={\Lambda }_u{\mathbf{f}}_i+{\delta }_{u,i},$$

In vector notation, the factor analytic model for these effects can be written as:

$$\mathbf{u}=({\mathbf{I}}_n\otimes \Lambda )\mathbf{f}+\delta ,$$

where $\mathbf{u}={({\mathbf{u}}_1^{\prime },\dots ,{\mathbf{u}}_n^{\prime })}^{\prime }$, $\mathbf{f}={({\mathbf{f}}_1^{\prime },\dots ,{\mathbf{f}}_n^{\prime })}^{\prime }$, and $\delta ={({\delta }_1^{\prime },\dots ,{\delta }_n^{\prime })}^{\prime }$. This results in the following covariance matrix of $\mathbf{u}$:

$$\mathrm{Cov}\left({\mathbf{u}}_i\right)={\mathbf{G}}_0={\Lambda }_u{\Lambda }_u^{\prime }+{\Psi }_u,$$

where ${\Psi }_u=\mathrm{Diag}\{{\psi }_{u,1},\dots ,{\psi }_{u,p}\}$ is a diagonal matrix of specific variances. Likewise, the residual vector ($ϵ$) can be decomposed into common and trait-specific factors. Consequently, the marginal distributions of ${\mathbf{u}}_i$ and ${ϵ}_i$ are:

$${\mathbf{u}}_i\stackrel{\mathrm{iid}}{\sim }N(\mathbf{0},{\Lambda }_u{\Lambda }_u^{\prime }+{\Psi }_u)\mathrm{and}{ϵ}_i\stackrel{\mathrm{iid}}{\sim }N(\mathbf{0},{\Lambda }_{ϵ}{\Lambda }_{ϵ}^{\prime }+{\Psi }_{ϵ}),$$

where “iid” stands for “independent and identically distributed”.

To model the additive genetic and residual (co)variance matrices using FA, it was assumed that Equation 2 holds for the vector of random additive genetic effects ($\mathbf{u}$) in Equation 1 and, likewise, for the vector of random residual effects ($ϵ$) such that

$$\mathbf{y}=\mathbf{X}\beta +\mathbf{Z}\left[({\mathbf{I}}_u\otimes {\Lambda }_u){\mathbf{f}}_u\right]+\mathbf{Z}{\delta }_u+\left[({\mathbf{I}}_{ϵ}\otimes {\Lambda }_{ϵ}){\mathbf{f}}_{ϵ}\right]+{\delta }_{ϵ},$$

where $\Lambda$, $\mathbf{X}$, $\mathbf{Z}$ and $\beta$ are as before and $\mathbf{f}$ and $\delta$ are interpreted as vectors of common and specific additive genetic effects, respectively, similarly for the random residual effects ($ϵ$). Combining the assumptions of the FA model described above with those of SMTM leads to the following random effects joint distribution,

$$p(ϵ,\mathbf{u})=N\left[ϵ\right|\mathbf{0},{\mathbf{I}}_{ϵ}\otimes ({\Lambda }_{ϵ}{\Lambda }_{ϵ}^{\prime }+{\Psi }_{ϵ})]N\left[\mathbf{u}\right|\mathbf{0},\mathbf{A}\otimes ({\Lambda }_u{\Lambda }_u^{\prime }+{\Psi }_u)],$$

where ${\Lambda }_u$ and ${\Lambda }_{ϵ}$ are the matrices of additive genetic and residual factor loads, respectively, ${\Psi }_u$ and ${\Psi }_{ϵ}$ are the diagonal matrices of specific additive genetic and residual variances, respectively. More details on the implementation of this model can be found in [17].

2.4. Recursive Models

Recursive models are a category of SEMs that postulate causal and unidirectional relationships between latent variables, unlike simultaneous models, which admit mutual feedback effects between traits. For more information on the procedure for REC effects decomposition, see [17,21]. Using REC, the vector $\mathbf{u}$ can be decomposed as

$${\mathbf{u}}_i={\Pi }_u{\mathbf{u}}_i+{\pounds }_{u,i},$$

where ${\mathbf{u}}_i$ is the random effect of the ith sire, ${\Pi }_u=\left\{{\pi }_{u,ij}\right\}$ is a strictly lower-triangular matrix; i.e., ${\pi }_{u,ij}=0$ for all $i\le j$, whose non-zero entries define recursive effects, and ${\pounds }_{u,i}$ is a vector of random effects whose (co)variance matrix is: $\mathrm{Cov}\left({\pounds }_u\right)=\mathbf{A}\otimes {\Gamma }_u$ where ${\Gamma }_u=\mathrm{Cov}\left({\pounds }_{u,i}\right)$ is a diagonal matrix. The reduced form model is ${\mathbf{u}}_i={(\mathbf{I}-{\Pi }_u)}^{-1}{\pounds }_{u,i}$ (Equation 4).

Similarly to the FA model, to implement the REC to model the additive genetic and the residual (co)variance matrices, it was assumed that Equation 4 holds the vector of random additive genetic effects ($\mathbf{u}$) in Equation 1, likewise for the vector of random residual effects ($ϵ$). Therefore, the additive genetic covariance matrix is equal to ${\mathbf{G}}_0$, as $\mathrm{Cov}\left({\mathbf{u}}_i\right)={\mathbf{G}}_0={(\mathbf{I}-{\Pi }_u)}^{-1}{\Gamma }_u{(\mathbf{I}-{\Pi }_u)}^{-1\prime }$. Similarly, using REC for residuals, ${\mathbf{R}}_0$ can be represented as $\mathrm{Cov}\left({ϵ}_i\right)={\mathbf{R}}_0={(\mathbf{I}-{\Pi }_{ϵ})}^{-1}{\Gamma }_{ϵ}{(\mathbf{I}-{\Pi }_{ϵ})}^{-1\prime }$, where ${\Pi }_{ϵ}$ is a strictly lower-triangular matrix defining recursive effects between model residuals, and ${\Gamma }_{ϵ}$ is diagonal. A fully recursive model (FRM) occurs when all lower-triangular entries of ${\Pi }_u$ and of ${\Pi }_{ϵ}$ are non-zero, i.e., all ${\pi }_{u,ij}$ and ${\pi }_{ϵ,ij}$ with $i>j$ are parameters to be estimated. This model has as many (co)dispersion parameters as the SMTM, and there is a one to one map from the unknowns in $\{{\Pi }_u,{\Gamma }_u,{\Pi }_{ϵ},{\Gamma }_{ϵ}\}$ to those in $\{{\mathbf{G}}_0,{\mathbf{R}}_0\}$. Therefore, the FRM is just a re-parameterization of SMTM; however, various degrees of parsimony can be obtained by setting some of the lower-triangular entries of ${\Pi }_u$ and ${\Pi }_{ϵ}$ to zero. Models in this study differed in how ${\Pi }_u$ and ${\Pi }_{ϵ}$ were structured.

2.5. Statistical Analyses

Bayesian inference and MCMC implementation: All models were implemented in a Bayesian framework using Gibbs sampling via the MTM package [24] in R software [25]. An initial chain of 1,500,000 iterations was run to assess convergence. The final chain length, burn-in period, and thinning interval were determined using the criteria of [26], implemented in the Bayesian Output Analysis (BOA) package [27]. Convergence was further evaluated through visual inspection of trace plots for variance components. Based on these diagnostics, posterior inference was based on 146,000 samples obtained after discarding the first 40,000 iterations as burn-in and retaining every 10th iteration thereafter.

Model comparison: Multiple alternative models were fitted to assess the effectiveness of parsimonious covariance structures. The standard multi-trait mixed model (SMTM) with unstructured (co)variance matrices served as the baseline, requiring estimation of 42 (co)variance parameters (21 for ${\mathbf{G}}_0$ and 21 for ${\mathbf{R}}_0$).

Factor analytic models: Four FA models were evaluated. The FA2F model specified two common factors for both additive genetic and residual (co)variance matrices, reducing the total number of parameters to 24 (12 genetic + 12 residual). The FA3F model extended this to three factors, increasing the parameter count to 36 (18 genetic + 18 residual). To evaluate the benefit of parsimony in only one (co)variance component, two additional models were fitted: FA2G applied two factors to ${\mathbf{G}}_0$ while maintaining the unstructured ${\mathbf{R}}_0$ (33 parameters total), whereas FA2R applied two factors to ${\mathbf{R}}_0$ while maintaining the unstructured ${\mathbf{G}}_0$ (33 parameters total).

Recursive models: Three recursive models were evaluated. First, a fully recursive model (FRM) was fitted, which is equivalent to SMTM in terms of number of parameters but reparameterizes the covariance structures using recursive paths. Based on FRM results, two reduced models were constructed. The REC1 model set the $ij$th entry of ${\Pi }_{ϵ}$ to zero if $|{\widehat{\pi }}_{ij}|/{\mathrm{sd}}_{ij}<1.96$ or $i\le j$, where ${\widehat{\pi }}_{ij}$ and ${\mathrm{sd}}_{ij}$ denote the posterior mean and standard deviation of the $ij$th entry of ${\Pi }_{ϵ}$ from FRM, respectively. This criterion resulted in six recursive effects being removed from ${\Pi }_{ϵ}$. The REC2 model removed any effect in ${\Pi }_u$ or ${\Pi }_{ϵ}$ whose posterior mean from FRM had absolute value less than 0.15 (an arbitrary threshold chosen to evaluate a more parsimonious structure), resulting in removal of 15 recursive effects (10 from ${\Pi }_{ϵ}$ and 5 from ${\Pi }_u$).

Model selection criteria: Models were compared using the deviance information criterion [DIC; [28,29], the posterior mean of the log-likelihood [Mean($\mathcal{L}$)], and the effective number of parameters [$p_D$; [28,29]. Lower DIC values indicate better model fit penalized for complexity, while $p_D$ quantifies model complexity. Additionally, Spearman and Pearson correlations between predicted sire breeding values across models were calculated to assess ranking stability.

3. Results

3.1. Genetic Parameters

Table 2, Table 3, and Table 4 present estimates of heritabilities (diagonal), genetic (above diagonal), and residual (below diagonal) correlations, obtained from the models SMTM, FA2G, and REC1, respectively. Standard errors obtained by these three models ranged from 0.05 to 0.12 for estimates of heritabilities. The Spearman and Pearson correlations of the estimated breeding values of all traits, comparing the same trait across models, ranged from 0.94 to 1.00.

For the genetic correlations, the standard errors obtained by these three models presented in Table 2, Table 3, and Table 4 ranged from 0.07 to 0.23. Considering the overlap within $\pm 1.0$ standard errors, differences between genetic correlation estimates obtained by SMTM (Table 2) and the REC1 (Table 4) model were not important. Differences between the genetic correlations estimates obtained by SMTM (Table 2) and FA2G (Table 3) models were observed for LMA and BW, BF and RF, BF and HH, RF and HH, and BW and HH.

For the residual correlations, the standard errors obtained by these three models presented in Table 2, Table 3, and Table 4 ranged from 0.01 to 0.07. Considering these standard error ranges, differences in residual correlation estimates obtained with the SMTM (Table 2) and REC1 (Table 4) models were observed only for SC and BF. The differences between the residual correlation estimates obtained by the SMTM (Table 2) and FA2G (Table 3) models were not significant.

3.2. Recursive Effects

Estimates of genetic (above diagonal) and residual (below diagonal) recursive effects from the FRM in Nellore cattle are shown in Table 5, assuming that the causation flows from the left (trait of the first column) to the right. In other words, a consequent variable never exerts causal influence (either directly or indirectly) on an antecedent variable that first exerts causal influence on it.

3.3. Model comparison

The deviance information criterion (DIC), the posterior mean of the log-likelihood (Mean(L)), and the estimated numbers of effective parameters ($p_D$) used to compare the proposed models are shown in Table 6. Although the SMTM model is considered fully parameterized for the additive genetic and residual matrices, the $p_D$ were largest in models that have many (co)dispersion parameters, as in the FA3F and FA2R, compared to the SMTM model; thus, the model with the smallest $p_D$ was FA2G.

The models with the best (smallest) DIC and Mean(L) were REC1 and SMTM, respectively. Conversely, FA3F was the worst model for these three comparison criteria (DIC, Mean(L), and $p_D$).

4. Discussion

This study aimed to evaluate the potential advantages of using structured equation models, such as recursive and factor-analytic models, for sire modeling in beef cattle genetic evaluations of ultrasound carcass and growth traits.

Given the magnitude of standard errors, differences between heritability estimates obtained by SMTM and the two best models (REC1 and FA2G) were only observed for BF and RF traits in the FA2G model, which yielded slightly lower estimates, except for SC in the REC1 model, but with the same pattern and within one standard deviation. The correlations (Spearman and Pearson) between estimated breeding values for each trait, across models were high and close to the unit, showing that selections decision would be basically the same using any of the methodologies considered here, despite the lower genetic variability (heritability estimates) of some traits in some models, because when estimating genetic breeding values the variances are generally fixed.

The FA2G model estimated genetic correlations with the same pattern as the SMTM model, but with a lower magnitude. It is suggested that only when genetic correlations are of lower magnitude can FA models properly estimate the true covariance structure; therefore, FA models are not a good alternative for estimating high correlations, as discussed by [30]. On the other hand, most of the genetic correlations estimated with the REC1 model were of slightly larger magnitude than those estimated by the SMTM, but yet very similar, indicating a good alternative for modeling additive genetic covariance matrices in beef cattle.

As expected, REC1 and FA2G modeled residual covariance matrices of beef cattle data very similarly to the SMTM model, except for the association between SC and BF in the REC1 model. Nevertheless, the residual correlation in the REC1 model was basically null $(0.09\pm 0.01)$. Considering the reasonable agreement between covariance parameters and breeding values estimated by the REC1 and FA2G models with the traditional model (SMTM), it could be suggested that these structural equations can be successfully used to reduce dimensionality of genetic evaluations, optimizing processing by a single analysis including all traits of interest and, consequently, providing faster and reliable estimation of breeding values compared to the traditional univariate analyses for each selection criteria. Thus, beef cattle breeding programs could run multiple-trait mixed-model analyses of large datasets nearly every day, such that as soon as the farmers send in new field data, they can get prompt results for selection decisions on their herds.

Recursive models act separately on $\mathbf{u}$ and $ϵ$, allowing specific patterns to be captured in each component (Table 3). The different patterns of genetic and residual correlations can be interpreted using estimates of recursive effects, and based on this information, irrelevant effects can be suppressed from the covariance matrices. The six recursive effects zeroed out in ${\Pi }_{ϵ}$ for the REC1 model were between RF and LMA, SC and LMA, BW and RF, HH and RF, SC and RF, and SC and HH. Therefore, the REC1 model saves these six recursive effects, which were virtually null, out of the total 21 parameters to be estimated, suggesting that this type of structural equation model can facilitate the estimation of breeding values for multiple traits in a single analysis.

In terms of the $p_D$ criterion, all recursive models (FRM, REC1, and REC2) had smaller $p_D$ values than the SMTM model; even for FRM, which has the same dimension as the SMTM, some benefits were observed when using this kind of structural equation model. For the FA models, $p_D$ favored the FA2G (two factors) model over all other FA and REC models. However, when we tested a model with three factors (FA3F), the $p_D$ was the largest. This was expected due to the number of parameters (p) of $\mathrm{Cov}\left(\mathbf{u}\right)={\mathbf{G}}_0$ in the FA, which can be calculated as $p=q+mq-m(m-1)/2$ parameters, where $q\times m$ are the dimensions of the $\Lambda$ matrix of factor loads (Equation 2). FA3F has 21 p in each matrix (additive genetic and residual), while FA2F has 17 p. In the SMTM, the p of ${\mathbf{G}}_0$ can be calculated as $p=q(q+1)/2$ parameters, where $q\times q$ are the dimensions of the ${\mathbf{G}}_0$ matrix, thus yielding 21 p in each matrix (additive genetic and residual) in the SMTM model. Depending on the data set, i.e., the number of records, a greater number of traits, and the most correlated traits, the FA models may be used as a special case of SEM to reduce the covariance matrices’ rank in the model and provide a more parsimonious estimation of genetic parameters compared to SMTM. In this paper, the FA3F model had a poorer fit to the data than the other models because this dataset has only six traits. Nevertheless, if the dataset had more traits, a model with more factors would probably work better at reducing the $p_D$. Generally, when more factors are included in the model, it better explains the relationships between traits; however, with fewer traits, fewer factors can be precisely estimated. Because we have only six traits and the number of records of SC used is limited, just because it is collected only in males ($n=1,340$), this could at least partially explain the worst performance in terms of DIC, $p_D$, and Mean(L) of the FA3F model. As discussed by [31] and [30], the FA model usually reduces the rank of covariance matrices. Still, the $\Psi$ matrix has to be close to zero $\mathrm{Cov}\left(\mathbf{u}\right)={\mathbf{G}}_0=\Lambda {\Lambda }^{\prime }+\Psi$, when specific effects are assumed absent. This results in a mixed-model formulation with covariance matrices that are not of full rank.

In the same sense, the Mean(L) and DIC demonstrated considerably better fits to the data for the models SMTM and REC1 than for all the other FA and REC models. Since DIC also accounts for dimensionality through a penalty ($p_D$), the REC1 model might be an alternative to the SMTM in beef cattle genetic evaluation schemes that involve multiple-trait estimates.

In general, even though FA2G was effective in terms of $p_D$, comparing all the criteria and parameters, the REC1 model showed a good alternative for this kind of dataset (e.g., six traits) because it had a smaller $p_D$ compared to the SMTM model, better DIC, and a reasonable Mean(L) compared with all other models tested. This model, when applied to the residual effect, can be used to obtain more parsimonious representations of the residual covariance matrix. This reduces the number of parameters to be estimated, and may yield useful interpretations of the underlying patterns of (co)variability. Exploration of genetic and residual relationships between traits can be detected using the REC1 model, helping to understand how selecting certain criteria can affect another trait. Moreover, unlike recursive models acting on phenotypes, these models do not impose cross-restrictions between the genetic and residual covariance matrices. Results from this study indicate that beef cattle growth and ultrasound carcass traits can be successfully analyzed using recursive models to estimate breeding values and identify a possible causal effect between latent variables, allowing for the modeling of complex phenomena while simultaneously reducing the dimensionality of the data.

5. Conclusions

Results from this study indicate that beef cattle growth and ultrasound carcass traits can be successfully analyzed using factor analysis or recursive models to estimate breeding values for multiple traits in a single analysis. Depending on the size of the data set, namely the number of records and traits, the factor analyses and recursive models applied independently to genetic and residual effects can reduce covariance matrices’ ranks and provide a more parsimonious estimation of genetic dispersion parameters compared to standard multiple-trait mixed models, especially if the covariances between the traits are of low magnitude. Exploration of genetic and residual relationships between traits via recursive models provides information to better understand how selection based on specific criteria can affect another trait through a causal effect.