The sommer package was developed to provide R users with a flexible univariate and multivariate linear mixed-model solver. Multi-environment trial (MET) analyses often need a genetic covariance structure among environments. A fully unstructured covariance (usm()) captures every pairwise environment relationship but requires \(q(q+1)/2-1\) parameters for \(q\) environments, which quickly becomes difficult to estimate reliably as the number of environments grows. Factor-analytic (FA) and reduced-rank (RR) models approximate that same covariance with far fewer parameters by assuming that genotype-by-environment interaction is driven by a small number of latent factors.
This vignette focuses on the fam() and rrm() covariance-shaping factors used inside vsm(), and on the loadings_mmes()/scores_mmes() helper functions used to extract and visualize their fitted quantities.
SECTION 1: Theory
- Why reduce the cross-environment covariance
- The factor-analytic (FA) model
- The reduced-rank (RR) model
- Choosing k and comparing to other structures
SECTION 2: Fitting factor-analytic and reduced-rank models
- Data
- Factor-analytic model with
fam() - Reduced-rank model with
rrm() - Comparing FA and RR
SECTION 3: Extracting loadings, scores, and diagnostic plots
- Loadings and specific variances
- Genotype scores
- Why and how to rotate loadings and scores
- Percentage of variance explained
- Loadings plot
- Genotype scores plot
- Latent regression plots
- Heatmap of the fitted genetic correlation matrix
SECTION 1: Theory
1) Why reduce the cross-environment covariance
Let \(q\) be the number of environments and \(G\) the \(q\times q\) genetic covariance matrix among them. An unstructured model estimates every variance and covariance directly, \(q(q+1)/2-1\) working parameters after removing the single overall scale owned by vsm(). As \(q\) grows, this saturated model becomes weakly identified relative to the available genotype replication, and its estimates become unstable. Diagonal (dsm()) and compound-symmetry (csm()) models are far more parsimonious but assume, respectively, no genetic correlation among environments or one common correlation everywhere. Factor-analytic and reduced-rank models sit between these extremes: they estimate genuine environment-specific covariance patterns while controlling the number of parameters through a rank \(k \ll q\).
2) The factor-analytic (FA) model
The FA model represents the environment covariance shape as
$$ K = \Lambda\Lambda^{\mathsf T} + \Psi, $$
where \(\Lambda\) is a \(q\times k\) loading matrix and \(\Psi\) is a diagonal matrix of environment-specific variances. Each environment’s genetic variance is split into a part explained by the \(k\) common latent factors (\(\Lambda\Lambda^{\mathsf T}\)) and a part unique to that environment (\(\Psi\)). vsm() still owns the single overall variance scale, so fam(x, k) reports a normalized shape with \(K_{11}=1\); the number of estimated working parameters for \(k\) factors and \(q\) environments is
$$ kq-\frac{k(k-1)}{2} + (q-1), $$
the first term for the (rotationally-constrained, lower-triangular) loadings and the second for the environment-specific variance ratios.
3) The reduced-rank (RR) model
The RR model uses the same loading structure but fixes every specific variance to zero:
$$ K = \Lambda\Lambda^{\mathsf T},\qquad \Lambda_{11}=1. $$
rrm(x, k) therefore has
$$ kq-\frac{k(k-1)}{2}-1 $$
free shape parameters, plus the overall scale supplied by vsm(): exactly \(q\) fewer parameters than heterogeneous fam(x, k). The covariance has rank \(k\), with no diagonal remainder, matching ASReml’s default rr() covariance family. The Henderson engine works with latent factors rather than inverting the singular environment covariance; direct inversion is also supported.
For an estimated common nonzero specific variance, use fam(x, k, variance="homogeneous"). This is the covariance family previously fitted by rrm(); the default fam() remains heterogeneous.
4) Choosing k and comparing to other structures
Increasing \(k\) enlarges the latent subspace, but RR remains singular while \(k<q\). In practice \(k\) is chosen small enough to remain identifiable and interpretable (often 1-3 factors). Zero-specific RR lies on the boundary of the FA family, so ordinary chi-squared likelihood-ratio reference distributions need not apply.
SECTION 2: Fitting factor-analytic and reduced-rank models
1) Data
We use the DT_h2 multi-environment potato yield dataset, which has 15 environments (Env, combining location and year) and 41 genotypes (Name).
library(sommer)
data(DT_h2, package="enhancer")
DT <- DT_h2
DT <- DT[with(DT, order(Env)), ]
length(unique(DT$Env))
## [1] 15
length(unique(DT$Name))
## [1] 41
head(DT)
## Name Env Loc Year Block y
## 67 MSL007-B CA.2011 CA 2011 CA.2011.2 5
## 105 MSL007-B CA.2011 CA 2011 CA.2011.1 6
## 308 MSK061-4 CA.2011 CA 2011 CA.2011.2 9
## 393 MSK061-4 CA.2011 CA 2011 CA.2011.1 10
## 469 MSR169-8Y CA.2011 CA 2011 CA.2011.1 11
## 471 NY148 CA.2011 CA 2011 CA.2011.1 11
2) Factor-analytic model with fam()
fitFA <- mmes(y ~ Env,
random = ~ vsm(fam(Env, 2), ism(Name)),
rcov = ~ units,
data = DT,
control = mmes_control(
# convergence = list(maxIters = 150),
diagnostics = list(verbose = FALSE)))
summary(fitFA)$varcomp
## term parameter estimate
## 1 vsm(fam(Env, 2), ism(Name)) loading[CA.2011,F1] 3.003812e+00
## 2 vsm(fam(Env, 2), ism(Name)) loading[CA.2012,F1] 1.505602e+00
## 3 vsm(fam(Env, 2), ism(Name)) loading[CA.2013,F1] 1.977072e+00
## 4 vsm(fam(Env, 2), ism(Name)) loading[FL.2011,F1] 4.316637e-01
## 5 vsm(fam(Env, 2), ism(Name)) loading[FL.2012,F1] 1.763942e-01
## 6 vsm(fam(Env, 2), ism(Name)) loading[FL.2013,F1] 9.708521e-01
## 7 vsm(fam(Env, 2), ism(Name)) loading[MI.2011,F1] 2.136873e+00
## 8 vsm(fam(Env, 2), ism(Name)) loading[MI.2012,F1] 1.441421e+00
## 9 vsm(fam(Env, 2), ism(Name)) loading[MI.2013,F1] 2.697790e+00
## 10 vsm(fam(Env, 2), ism(Name)) loading[MO.2011,F1] 8.359788e-01
## 11 vsm(fam(Env, 2), ism(Name)) loading[MO.2012,F1] 2.298807e+00
## 12 vsm(fam(Env, 2), ism(Name)) loading[MO.2013,F1] 1.177742e+00
## 13 vsm(fam(Env, 2), ism(Name)) loading[NY.2011,F1] 9.778560e-01
## 14 vsm(fam(Env, 2), ism(Name)) loading[NY.2012,F1] 1.369488e+00
## 15 vsm(fam(Env, 2), ism(Name)) loading[NY.2013,F1] 8.897588e-01
## 16 vsm(fam(Env, 2), ism(Name)) loading[CA.2012,F2] 1.734993e-01
## 17 vsm(fam(Env, 2), ism(Name)) loading[CA.2013,F2] -5.764190e-01
## 18 vsm(fam(Env, 2), ism(Name)) loading[FL.2011,F2] 1.197546e+00
## 19 vsm(fam(Env, 2), ism(Name)) loading[FL.2012,F2] 5.521793e-01
## 20 vsm(fam(Env, 2), ism(Name)) loading[FL.2013,F2] 1.153011e-01
## 21 vsm(fam(Env, 2), ism(Name)) loading[MI.2011,F2] 1.158865e+00
## 22 vsm(fam(Env, 2), ism(Name)) loading[MI.2012,F2] 8.069280e-01
## 23 vsm(fam(Env, 2), ism(Name)) loading[MI.2013,F2] 9.415421e-01
## 24 vsm(fam(Env, 2), ism(Name)) loading[MO.2011,F2] -2.449628e-01
## 25 vsm(fam(Env, 2), ism(Name)) loading[MO.2012,F2] 1.484394e+00
## 26 vsm(fam(Env, 2), ism(Name)) loading[MO.2013,F2] 6.425122e-01
## 27 vsm(fam(Env, 2), ism(Name)) loading[NY.2011,F2] 2.064666e+00
## 28 vsm(fam(Env, 2), ism(Name)) loading[NY.2012,F2] 4.524148e-01
## 29 vsm(fam(Env, 2), ism(Name)) loading[NY.2013,F2] 2.806597e+00
## 30 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2011] 6.678659e+00
## 31 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2012] 3.048041e+00
## 32 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2013] 3.609733e+00
## 33 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2011] 4.942937e-02
## 34 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2012] 6.678659e-06
## 35 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2013] 6.678659e-06
## 36 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2011] 3.051821e+00
## 37 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2012] 2.316374e+00
## 38 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2013] 1.126774e+01
## 39 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2011] 6.678659e-06
## 40 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2012] 9.013260e+00
## 41 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2013] 2.712617e-02
## 42 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2011] 6.678659e-06
## 43 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2012] 6.678659e-06
## 44 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2013] 3.876432e-01
## 45 vsm(ism(units)) sigma2 3.968801e+00
## StdError Zratio
## 1 0.81947698 3.665523e+00
## 2 0.45669326 3.296748e+00
## 3 0.57092898 3.462904e+00
## 4 0.51065302 8.453170e-01
## 5 0.30859879 5.715973e-01
## 6 0.00000000 NA
## 7 0.76673577 2.786974e+00
## 8 0.49171676 2.931405e+00
## 9 0.87054119 3.098981e+00
## 10 0.29853230 2.800296e+00
## 11 0.88273930 2.604175e+00
## 12 0.37042561 3.179430e+00
## 13 0.75725190 1.291322e+00
## 14 0.00000000 NA
## 15 0.96093484 9.259304e-01
## 16 0.64303666 2.698124e-01
## 17 0.81860003 -7.041522e-01
## 18 0.47310686 2.531237e+00
## 19 0.29008269 1.903524e+00
## 20 0.50234825 2.295243e-01
## 21 0.81555379 1.420954e+00
## 22 0.54432809 1.482429e+00
## 23 1.13285715 8.311216e-01
## 24 0.37122106 -6.598840e-01
## 25 1.00438814 1.477908e+00
## 26 0.46767450 1.373845e+00
## 27 0.59355627 3.478468e+00
## 28 0.46515991 9.726005e-01
## 29 0.63598977 4.412960e+00
## 30 3.47078776 1.924249e+00
## 31 1.32905103 2.293397e+00
## 32 1.87403484 1.926182e+00
## 33 0.77799207 6.353454e-02
## 34 0.00000000 NA
## 35 0.00000000 NA
## 36 2.00400446 1.522861e+00
## 37 1.09940707 2.106930e+00
## 38 3.87667802 2.906544e+00
## 39 0.00000000 NA
## 40 3.28758564 2.741605e+00
## 41 0.58599217 4.629101e-02
## 42 1.22718043 5.442279e-06
## 43 0.08649756 7.721210e-05
## 44 1.93500199 2.003322e-01
## 45 0.24799917 1.600328e+01
3) Reduced-rank model with rrm()
fitRR <- mmes(y ~ Env,
random = ~ vsm(rrm(Env, 2), ism(Name)),
rcov = ~ units,
data = DT,
control = mmes_control(
# convergence = list(maxIters = 150),
diagnostics = list(verbose = FALSE)))
summary(fitRR)$varcomp
## term parameter estimate
## 1 vsm(rrm(Env, 2), ism(Name)) loading[CA.2011,F1] 3.426629389
## 2 vsm(rrm(Env, 2), ism(Name)) loading[CA.2012,F1] 1.526552388
## 3 vsm(rrm(Env, 2), ism(Name)) loading[CA.2013,F1] 1.953275332
## 4 vsm(rrm(Env, 2), ism(Name)) loading[FL.2011,F1] 0.546451896
## 5 vsm(rrm(Env, 2), ism(Name)) loading[FL.2012,F1] 0.251931289
## 6 vsm(rrm(Env, 2), ism(Name)) loading[FL.2013,F1] 0.741325547
## 7 vsm(rrm(Env, 2), ism(Name)) loading[MI.2011,F1] 2.425919747
## 8 vsm(rrm(Env, 2), ism(Name)) loading[MI.2012,F1] 1.519441403
## 9 vsm(rrm(Env, 2), ism(Name)) loading[MI.2013,F1] 3.841674401
## 10 vsm(rrm(Env, 2), ism(Name)) loading[MO.2011,F1] 0.856346428
## 11 vsm(rrm(Env, 2), ism(Name)) loading[MO.2012,F1] 2.914235717
## 12 vsm(rrm(Env, 2), ism(Name)) loading[MO.2013,F1] 1.284434347
## 13 vsm(rrm(Env, 2), ism(Name)) loading[NY.2011,F1] 1.780391579
## 14 vsm(rrm(Env, 2), ism(Name)) loading[NY.2012,F1] 1.240314199
## 15 vsm(rrm(Env, 2), ism(Name)) loading[NY.2013,F1] 1.400196358
## 16 vsm(rrm(Env, 2), ism(Name)) loading[CA.2012,F2] 0.515297404
## 17 vsm(rrm(Env, 2), ism(Name)) loading[CA.2013,F2] -1.453686127
## 18 vsm(rrm(Env, 2), ism(Name)) loading[FL.2011,F2] 0.837850110
## 19 vsm(rrm(Env, 2), ism(Name)) loading[FL.2012,F2] 0.307601042
## 20 vsm(rrm(Env, 2), ism(Name)) loading[FL.2013,F2] 0.329087672
## 21 vsm(rrm(Env, 2), ism(Name)) loading[MI.2011,F2] 1.534871002
## 22 vsm(rrm(Env, 2), ism(Name)) loading[MI.2012,F2] 0.777469166
## 23 vsm(rrm(Env, 2), ism(Name)) loading[MI.2013,F2] -0.935909371
## 24 vsm(rrm(Env, 2), ism(Name)) loading[MO.2011,F2] -0.081619919
## 25 vsm(rrm(Env, 2), ism(Name)) loading[MO.2012,F2] 2.815485404
## 26 vsm(rrm(Env, 2), ism(Name)) loading[MO.2013,F2] -0.003401437
## 27 vsm(rrm(Env, 2), ism(Name)) loading[NY.2011,F2] 2.094595786
## 28 vsm(rrm(Env, 2), ism(Name)) loading[NY.2012,F2] 0.408297756
## 29 vsm(rrm(Env, 2), ism(Name)) loading[NY.2013,F2] 1.798032047
## 30 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2011] 0.000000000
## 31 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2012] 0.000000000
## 32 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2013] 0.000000000
## 33 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2011] 0.000000000
## 34 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2012] 0.000000000
## 35 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2013] 0.000000000
## 36 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2011] 0.000000000
## 37 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2012] 0.000000000
## 38 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2013] 0.000000000
## 39 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2011] 0.000000000
## 40 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2012] 0.000000000
## 41 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2013] 0.000000000
## 42 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2011] 0.000000000
## 43 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2012] 0.000000000
## 44 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2013] 0.000000000
## 45 vsm(ism(units)) sigma2 5.765728207
## StdError Zratio
## 1 0.5181774 6.612849802
## 2 0.3345536 4.562953341
## 3 0.4655191 4.195908133
## 4 0.4555556 1.199528545
## 5 0.3320777 0.758651663
## 6 0.4351673 1.703541300
## 7 0.5839836 4.154088596
## 8 0.3542538 4.289132574
## 9 0.5635509 6.816907698
## 10 0.4691616 1.825269493
## 11 0.7340443 3.970108585
## 12 0.3701975 3.469592076
## 13 0.7022109 2.535408507
## 14 0.3445412 3.599901223
## 15 0.4924655 2.843237407
## 16 0.3599574 1.431551108
## 17 0.4471123 -3.251277126
## 18 0.4781954 1.752108347
## 19 0.3596013 0.855394615
## 20 0.5560246 0.591858107
## 21 0.4771067 3.217039317
## 22 0.3280407 2.370038700
## 23 0.5031489 -1.860104176
## 24 0.5021126 -0.162553029
## 25 0.5913268 4.761301831
## 26 0.4217239 -0.008065555
## 27 0.5944575 3.523541838
## 28 0.3648052 1.119221371
## 29 0.4102523 4.382747164
## 30 0.0000000 NA
## 31 0.0000000 NA
## 32 0.0000000 NA
## 33 0.0000000 NA
## 34 0.0000000 NA
## 35 0.0000000 NA
## 36 0.0000000 NA
## 37 0.0000000 NA
## 38 0.0000000 NA
## 39 0.0000000 NA
## 40 0.0000000 NA
## 41 0.0000000 NA
## 42 0.0000000 NA
## 43 0.0000000 NA
## 44 0.0000000 NA
## 45 0.2929790 19.679664485
4) Comparing FA and RR
Both models use rank 2, but RR fixes all specific variances at zero. Compare their fitted likelihoods:
c(logLik_FA = tail(fitFA$llik, 1), logLik_RR = tail(fitRR$llik, 1))
## logLik_FA1 logLik_FA2 logLik_FA3 logLik_FA4 logLik_FA5 logLik_FA6
## -29.932612 25.601513 45.949417 58.736232 64.257852 68.370377
## logLik_FA7 logLik_FA8 logLik_FA9 logLik_FA10 logLik_FA11 logLik_FA12
## 77.945279 81.125024 82.256633 85.330110 86.026199 86.553244
## logLik_FA13 logLik_FA14 logLik_FA15 logLik_FA16 logLik_FA17 logLik_FA18
## 86.842560 87.003094 87.137520 87.157912 87.233575 87.270951
## logLik_FA19 logLik_FA20 logLik_FA21 logLik_FA22 logLik_FA23 logLik_FA24
## 87.288354 87.296045 87.299083 87.299855 87.302477 87.302907
## logLik_FA25 logLik_FA26 logLik_RR1 logLik_RR2 logLik_RR3 logLik_RR4
## 87.303056 87.303063 -108.865994 -39.971428 -31.808583 -15.413570
## logLik_RR5 logLik_RR6 logLik_RR7 logLik_RR8 logLik_RR9 logLik_RR10
## 2.984253 12.932587 21.669742 27.588034 27.615913 32.487973
## logLik_RR11 logLik_RR12 logLik_RR13 logLik_RR14 logLik_RR15 logLik_RR16
## 34.058806 34.395555 34.541674 34.614326 34.667450 34.687818
## logLik_RR17 logLik_RR18 logLik_RR19 logLik_RR20 logLik_RR21 logLik_RR22
## 34.708731 35.012824 35.052648 35.082935 35.109266 35.132370
## logLik_RR23 logLik_RR24 logLik_RR25 logLik_RR26 logLik_RR27 logLik_RR28
## 35.152689 35.170610 35.186411 35.200371 35.212687 35.223572
## logLik_RR29 logLik_RR30
## 35.233175 35.241658
The less restrictive FA model can capture genetic variation outside the common-factor subspace. A formal likelihood-ratio test requires accounting for the zero-variance boundary; do not interpret a naive chi-squared test as an ordinary interior-parameter comparison.
SECTION 3: Extracting loadings, scores, and diagnostic plots
1) Loadings and specific variances
loadings_mmes() reconstructs the fitted loadings and specific variances. Here we request varianceScale = FALSE explicitly to obtain the normalized loadings \(\Lambda\) and specific variances \(\Psi\), so that \(\sigma^2(\Lambda\Lambda^{\mathsf T}+\Psi)\) reproduces the fitted covariance matrix. The default rotation = TRUE returns rotated loadings; rotation does not change this covariance.
faInfo <- loadings_mmes(fitFA, varianceScale = FALSE, rotation = TRUE)
round(faInfo$loadings, 3)
## F1 F2
## CA.2011 0.673 0.349
## CA.2012 0.357 0.136
## CA.2013 0.376 0.359
## FL.2011 0.236 -0.218
## FL.2012 0.104 -0.103
## FL.2013 0.231 0.087
## MI.2011 0.613 -0.011
## MI.2012 0.417 -0.013
## MI.2013 0.714 0.103
## MO.2011 0.159 0.152
## MO.2012 0.688 -0.065
## MO.2013 0.339 -0.007
## NY.2011 0.459 -0.349
## NY.2012 0.359 0.058
## NY.2013 0.526 -0.525
round(faInfo$specific, 3)
## CA.2011 CA.2012 CA.2013 FL.2011 FL.2012 FL.2013 MI.2011 MI.2012 MI.2013 MO.2011
## 0.425 0.194 0.230 0.003 0.000 0.000 0.194 0.148 0.718 0.000
## MO.2012 MO.2013 NY.2011 NY.2012 NY.2013
## 0.574 0.002 0.000 0.000 0.025
faInfo$sigma2
## [1] 15.70155
rrInfo <- loadings_mmes(fitRR, varianceScale = FALSE, rotation = TRUE)
round(rrInfo$specific, 3)
## CA.2011 CA.2012 CA.2013 FL.2011 FL.2012 FL.2013 MI.2011 MI.2012 MI.2013 MO.2011
## 0 0 0 0 0 0 0 0 0 0
## MO.2012 MO.2013 NY.2011 NY.2012 NY.2013
## 0 0 0 0 0
Notice that rrInfo$specific is exactly zero in every environment, unlike the estimated heterogeneous values in faInfo$specific. With fam(..., variance="homogeneous"), the extracted specifics are equal but generally nonzero.
2) Genotype scores
scores_mmes() predicts each genotype’s position on the \(k\) latent factors from its BLUPs, the fitted loadings, and the fitted covariance. We use the default regression method and varianceScale = TRUE: scores are then on the standardized latent-factor scale, with the overall genetic scale absorbed into the corresponding loadings. The term “regression” here describes the score estimator, not a new regression fitted to the plotted points.
faScores <- scores_mmes(fitFA, method = "regression",
varianceScale = TRUE, rotation = TRUE)
head(faScores)
## F1 F2
## A01143-3C 1.18406048 1.3806693
## AC00206-2W 0.00846617 -0.3384607
## AC01151-5W 0.54826010 -0.1007708
## AC03433-1W -1.02576139 -0.1051194
## AC03452-2W 2.00158652 -0.1861375
## AC05153-1W -1.17569818 -0.3737137
3) Why and how to rotate loadings and scores
The fitted covariance identifies a latent subspace, not a unique set of factor axes. For any orthogonal \(k\times k\) matrix \(Q\), with \(Q Q^{\mathsf T}=I\),
$$ (\Lambda Q)(\Lambda Q)^{\mathsf T}=\Lambda\Lambda^{\mathsf T}. $$
Consequently, many loading matrices describe exactly the same model. The lower-triangular constraints used during fitting identify one parameterization, but its axes depend on the ordering of environments and need not be convenient for interpretation. Rotation at extraction is a change of coordinates, not a refit: it leaves the likelihood, specific variances, genetic covariance, and fitted genotype-by-environment effects unchanged.
Both loadings_mmes() and scores_mmes() accept the logical argument rotation:
rotation = FALSEretains the fitted, lower-triangular loading orientation.rotation = TRUE(the default) rotates using the right singular vectors of the loading matrix. If \(\Lambda=A D V^{\mathsf T}\) is its singular value decomposition, sommer uses \(Q=-V\) and returns \(\Lambda_* = \Lambda Q\). The loading columns are then orthogonal, ordered by decreasing squared length. This is an SVD rotation, not a varimax or oblique rotation; no rotation name or custom matrix is accepted by this argument.
The rotated first factor captures the largest share of the common genetic variance; subsequent factors capture decreasing shares. When first-factor loadings have the same sign across environments, that factor can describe broad performance, while later factors describe contrasting responses or crossover interaction. This interpretation is data-dependent: the first factor is not necessarily an overall-performance axis, and later factors are not automatically particular environmental variables.
Always rotate loadings and scores together. Writing \(L=\sqrt{\sigma^2}\Lambda\) for the loadings in trait units and \(F\) for the genotype-by-factor score matrix, the common genetic prediction is
$$ \widehat U_{\mathrm{common}} = F L^{\mathsf T}. $$
With \(L_* = LQ\) and \(F_* = FQ\), this product is unchanged. Mixing rotated scores with unrotated loadings gives incorrect reconstructions and misleading plots. Keep term, method, and varianceScale consistent across extraction calls as well. Here the normalized loadings in faInfo are converted to trait units explicitly to match faScores:
faUnrotated <- loadings_mmes(fitFA, varianceScale = FALSE,
rotation = FALSE)
scoresUnrotated <- scores_mmes(fitFA, method = "regression",
varianceScale = TRUE, rotation = FALSE)
L <- sqrt(faInfo$sigma2) * faInfo$loadings
LUnrotated <- sqrt(faUnrotated$sigma2) * faUnrotated$loadings
common <- faScores %*% t(L)
commonUnrotated <- scoresUnrotated %*% t(LUnrotated)
stopifnot(
isTRUE(all.equal(tcrossprod(L), tcrossprod(LUnrotated))),
isTRUE(all.equal(common, commonUnrotated)),
isTRUE(all.equal(faInfo$specific, faUnrotated$specific))
)
round(crossprod(L), 3)
## F1 F2
## F1 49.134 0.000
## F2 0.000 12.157
The same arguments and reasoning apply to fitRR. For one factor there is no multi-axis rotation, only a possible sign change. A factor’s sign is arbitrary: reversing both a loading column and its score column changes neither the covariance nor the common predictions. Thus a negative score is not intrinsically “poor performance”; interpret it together with the loadings. Across separate fits, axes can also swap or become unstable when singular values are nearly equal, so raw score signs and labels should not be compared without aligning the axes.
This extraction argument is distinct from vsm(..., Gu = Gu, rotation = TRUE), which requests Lee–van der Werf rotation of a relationship term during fitting, not the post-fit orientation of FA/RR factors.
4) Percentage of variance explained
A standard factor-analytic diagnostic is the proportion of total genetic variance attributed to each latent factor:
varPerFactor <- colSums(faInfo$loadings^2)
totalVar <- sum(faInfo$loadings^2) + sum(faInfo$specific)
propExplained <- varPerFactor / totalVar
round(100 * propExplained, 1)
## F1 F2
## 48.8 12.1
These percentages use the normalized loadings and specific variances, so the common overall scale cancels. Their sum is the percentage of total genetic variance explained by the common factors, not necessarily 100%. Rotation redistributes variance between factors but preserves that sum. Under the default SVD rotation, factors are ordered by decreasing contribution; the percentages refer to variance across environments, not phenotypic variance or genotype-specific prediction accuracy.
5) Loadings plot
Plotting the loadings by environment shows which environments are most associated with each latent factor.
barplot(t(faInfo$loadings), beside = TRUE,
col = c("steelblue4", "tomato"),
las = 2, cex.names = 0.7,
ylab = "Loading",
main = "Factor-analytic loadings by environment")
legend("topright", legend = colnames(faInfo$loadings),
fill = c("steelblue4", "tomato"), bty = "n")
6) Genotype scores plot
This is a plot of genotype scores only, rather than a biplot containing both genotype scores and environment loadings. Genotypes close together have similar common-factor responses, although their environment-specific effects may differ. The signs and directions must be read alongside the loading plot.
plot(faScores[,2] ~ faScores[,1],
xlab = "Factor 1 score", ylab = "Factor 2 score",
main = "Genotype scores")
text(faScores[,2] ~ faScores[,1], labels = rownames(faScores),
cex = 0.6, pos = 1)
abline(h = 0, v = 0, lty = 3)
7) Latent regression plots
In the FA literature, a latent regression plot displays genotype performance against an environmental coordinate estimated by the model: a factor loading. Unlike regression on measured rainfall or temperature, the horizontal coordinate is latent and learned from the genetic covariance. Each genotype has its own regression slope, given by its factor score. Such plots help distinguish broad adaptation from sensitivity to the environmental contrasts captured by a factor; see Smith and Cullis (2018) for FA-based selection tools and interpretation of performance and stability.
For genotype \(g\) in environment \(e\), the genetic effect is decomposed as
$$ \widehat u_{ge} = \sum_{j=1}^k \widehat f_{gj} L_{ej} + \widehat\delta_{ge}, $$
where \(\widehat\delta_{ge}\) is the predicted environment-specific remainder. In a one-factor model the common response is the straight line \(\widehat f_{g1} L_{e1}\) through the origin. With multiple factors, plotting the full genetic effect against only one loading also includes contributions from the other factors. To isolate factor \(j\), subtract those contributions first:
$$ \widehat u_{ge}^{(j)} = \widehat u_{ge}
- \sum_{h\ne j}\widehat f_{gh} L_{eh} = \widehat f_{gj} L_{ej} + \widehat\delta_{ge}. $$
The following base-R plots show these adjusted genetic BLUPs as points and the common-factor response as a line. We select four genotypes spanning the rotated first-factor scores to keep the display readable. Their selection is illustrative, not a ranking rule.
U <- fitFA$uList[[faInfo$term]]
U <- U[rownames(faScores), rownames(L), drop = FALSE]
specificBLUP <- U - common
scoreOrder <- order(faScores[, 1])
selected <- scoreOrder[unique(round(seq(1, length(scoreOrder),
length.out = 4)))]
genotypeColors <- c("steelblue4", "tomato3", "darkgreen", "goldenrod3")
local({
oldPar <- par(mfrow = c(1, ncol(L)), mar = c(4.5, 4.5, 3, 1))
on.exit(par(oldPar))
for (factorIndex in seq_len(ncol(L))) {
loading <- L[, factorIndex]
factorResponse <- outer(faScores[selected, factorIndex], loading)
adjusted <- factorResponse + specificBLUP[selected, , drop = FALSE]
plot(range(loading), range(c(adjusted, factorResponse, 0)),
type = "n", xlab = paste("Environment loading:", colnames(L)[factorIndex]),
ylab = "Adjusted genetic effect", main = paste("Latent regression", factorIndex))
abline(h = 0, v = 0, lty = 3, col = "grey70")
for (genotypeIndex in seq_along(selected)) {
points(loading, adjusted[genotypeIndex, ], pch = 16,
col = genotypeColors[genotypeIndex])
abline(a = 0, b = faScores[selected[genotypeIndex], factorIndex],
col = genotypeColors[genotypeIndex], lwd = 2)
}
legend("topleft", legend = rownames(faScores)[selected],
col = genotypeColors[seq_along(selected)], pch = 16, lty = 1,
cex = 0.7, bty = "n")
}
})
Each point represents an environment, located at its loading on that factor. A steep slope indicates a strong response to that latent environmental contrast; a nearly flat slope indicates little response along that axis. Opposite slopes indicate contrasting responses. Points away from a genotype’s line reflect the environment-specific remainder after removing the other factors, not raw observation residuals. The lines use model-derived scores; they are not ordinary least-squares regressions fitted to these points.
These are plots of genetic deviations, not absolute yield: the fixed environment effects are omitted, so all common-factor lines pass through zero and their crossings alone do not establish differences in total yield. For selection, consider the full predicted genetic effects (or predictions including fixed effects), all retained factors, specific effects, and prediction uncertainty. The points are shrunken BLUPs, and uncertainty in both loadings and scores is not shown. The default regression scores and Bartlett scores need not give the same slopes; use one method consistently. To draw analogous RR plots, extract the RR loadings and scores with the same rotation and scale choices and use the RR term’s BLUP matrix.
8) Heatmap of the fitted genetic correlation matrix
The normalized loadings and specific variances returned by loadings_mmes() can be combined directly into the fitted genetic covariance, and then converted to a correlation matrix for visualization.
Sigma <- faInfo$sigma2 * (
faInfo$loadings %*% t(faInfo$loadings) + diag(faInfo$specific)
)
corMat <- cov2cor(Sigma)
stopifnot(isTRUE(all.equal(unname(Sigma),
unname(fitFA$theta[[faInfo$term]]))))
heatmap(corMat, symm = TRUE,
main = "Fitted genetic correlation among environments")
Environments with a strong positive fitted correlation cluster together in the heatmap, while environments explained by different latent factors, or with a large specific variance, appear less correlated with the rest.
Literature
Covarrubias-Pazaran G. 2016. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6):1-15.
Smith AB, Cullis BR, Thompson R. 2001. Analyzing variety by environment data using multiplicative mixed models and expectations of trials. Biometrics 57(4).
Thompson R, Cullis B, Smith A, Gilmour A. 2003. A sparse implementation of the average information algorithm for factor analytic and reduced rank variance models. Australian & New Zealand Journal of Statistics 45(4).
Smith AB, Cullis BR. 2018. Plant breeding selection tools built on factor analytic mixed models for multi-environment trial data. Euphytica 214:143. https://doi.org/10.1007/s10681-018-2220-5.
Gilmour AR, Thompson R, Cullis BR. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51:1440-1450.