Factor analytic and reduced rank models in sommer

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

  1. Why reduce the cross-environment covariance
  2. The factor-analytic (FA) model
  3. The reduced-rank (RR) model
  4. Choosing k and comparing to other structures

SECTION 2: Fitting factor-analytic and reduced-rank models

  1. Data
  2. Factor-analytic model with fam()
  3. Reduced-rank model with rrm()
  4. Comparing FA and RR

SECTION 3: Extracting loadings, scores, and diagnostic plots

  1. Loadings and specific variances
  2. Genotype scores
  3. Why and how to rotate loadings and scores
  4. Percentage of variance explained
  5. Loadings plot
  6. Genotype scores plot
  7. Latent regression plots
  8. 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,
              nIters = 150, verbose = FALSE,
              data = DT)
## Warning in mmes(y ~ Env, random = ~vsm(fam(Env, 2), ism(Name)), rcov = ~units,
## : AI-first optimization stalled at the last accepted parameters
## (recovery-stalled). Consider different starting values; convergence was not
## certified.
summary(fitFA)$varcomp
##                           term                  parameter      estimate
## 1  vsm(fam(Env, 2), ism(Name))        loading[CA.2011,F1]  3.001750e+00
## 2  vsm(fam(Env, 2), ism(Name))        loading[CA.2012,F1]  1.505817e+00
## 3  vsm(fam(Env, 2), ism(Name))        loading[CA.2013,F1]  1.972337e+00
## 4  vsm(fam(Env, 2), ism(Name))        loading[FL.2011,F1]  4.351413e-01
## 5  vsm(fam(Env, 2), ism(Name))        loading[FL.2012,F1]  1.759811e-01
## 6  vsm(fam(Env, 2), ism(Name))        loading[FL.2013,F1]  9.730630e-01
## 7  vsm(fam(Env, 2), ism(Name))        loading[MI.2011,F1]  2.138747e+00
## 8  vsm(fam(Env, 2), ism(Name))        loading[MI.2012,F1]  1.443908e+00
## 9  vsm(fam(Env, 2), ism(Name))        loading[MI.2013,F1]  2.701687e+00
## 10 vsm(fam(Env, 2), ism(Name))        loading[MO.2011,F1]  8.350447e-01
## 11 vsm(fam(Env, 2), ism(Name))        loading[MO.2012,F1]  2.304051e+00
## 12 vsm(fam(Env, 2), ism(Name))        loading[MO.2013,F1]  1.179139e+00
## 13 vsm(fam(Env, 2), ism(Name))        loading[NY.2011,F1]  9.810314e-01
## 14 vsm(fam(Env, 2), ism(Name))        loading[NY.2012,F1]  1.370377e+00
## 15 vsm(fam(Env, 2), ism(Name))        loading[NY.2013,F1]  8.943842e-01
## 16 vsm(fam(Env, 2), ism(Name))        loading[CA.2012,F2]  1.708272e-01
## 17 vsm(fam(Env, 2), ism(Name))        loading[CA.2013,F2] -5.747152e-01
## 18 vsm(fam(Env, 2), ism(Name))        loading[FL.2011,F2]  1.196518e+00
## 19 vsm(fam(Env, 2), ism(Name))        loading[FL.2012,F2]  5.541338e-01
## 20 vsm(fam(Env, 2), ism(Name))        loading[FL.2013,F2]  1.105096e-01
## 21 vsm(fam(Env, 2), ism(Name))        loading[MI.2011,F2]  1.155608e+00
## 22 vsm(fam(Env, 2), ism(Name))        loading[MI.2012,F2]  8.024194e-01
## 23 vsm(fam(Env, 2), ism(Name))        loading[MI.2013,F2]  9.323091e-01
## 24 vsm(fam(Env, 2), ism(Name))        loading[MO.2011,F2] -2.448994e-01
## 25 vsm(fam(Env, 2), ism(Name))        loading[MO.2012,F2]  1.474594e+00
## 26 vsm(fam(Env, 2), ism(Name))        loading[MO.2013,F2]  6.403232e-01
## 27 vsm(fam(Env, 2), ism(Name))        loading[NY.2011,F2]  2.061535e+00
## 28 vsm(fam(Env, 2), ism(Name))        loading[NY.2012,F2]  4.494236e-01
## 29 vsm(fam(Env, 2), ism(Name))        loading[NY.2013,F2]  2.808992e+00
## 30 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2011]  6.688974e+00
## 31 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2012]  3.046572e+00
## 32 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2013]  3.622008e+00
## 33 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2011]  4.733796e-02
## 34 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2012]  1.493871e-05
## 35 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2013]  8.439076e-06
## 36 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2011]  3.051284e+00
## 37 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2012]  2.315665e+00
## 38 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2013]  1.126748e+01
## 39 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2011]  1.340884e-05
## 40 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2012]  9.013077e+00
## 41 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2013]  2.733208e-02
## 42 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2011]  1.620564e-04
## 43 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2012]  3.153882e-04
## 44 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2013]  3.729244e-01
## 45             vsm(ism(units))                     sigma2  3.968861e+00
##     StdError        Zratio
## 1  0.8227142  3.648594e+00
## 2  0.4632823  3.250323e+00
## 3  0.5734481  3.439434e+00
## 4  0.5192763  8.379763e-01
## 5  0.3268001  5.384977e-01
## 6  0.3709634  2.623070e+00
## 7  0.7771504  2.752038e+00
## 8  0.4941013  2.922292e+00
## 9  0.9125495  2.960592e+00
## 10 0.3912243  2.134440e+00
## 11 0.8814928  2.613806e+00
## 12 0.3730415  3.160878e+00
## 13 0.7889077  1.243531e+00
## 14 0.3412853  4.015342e+00
## 15 0.9726583  9.195257e-01
## 16 0.6435562  2.654425e-01
## 17 0.8193255 -7.014492e-01
## 18 0.4728979  2.530183e+00
## 19 0.2961967  1.870830e+00
## 20 0.5044144  2.190848e-01
## 21 0.8168394  1.414731e+00
## 22 0.5449630  1.472429e+00
## 23 1.1356595  8.209407e-01
## 24 0.4929317 -4.968221e-01
## 25 1.0050254  1.467220e+00
## 26 0.4678678  1.368599e+00
## 27 0.5963306  3.457034e+00
## 28 0.4773092  9.415775e-01
## 29 0.6540804  4.294568e+00
## 30 3.4973130  1.912604e+00
## 31 1.3314269  2.288201e+00
## 32 1.8766161  1.930074e+00
## 33 0.7816303  6.056310e-02
## 34 0.5206526  2.869229e-05
## 35 0.8313128  1.015150e-05
## 36 2.0088326  1.518934e+00
## 37 1.1003868  2.104410e+00
## 38 3.8950417  2.892775e+00
## 39 0.9308603  1.440478e-05
## 40 3.2881191  2.741104e+00
## 41 0.5878004  4.649892e-02
## 42 1.2734167  1.272611e-04
## 43 0.5768502  5.467419e-04
## 44 2.0163249  1.849526e-01
## 45 0.2644497  1.500800e+01

3) Reduced-rank model with rrm()

fitRR <- mmes(y ~ Env,
              random = ~ vsm(rrm(Env, 2), ism(Name)),
              rcov = ~ units,
              nIters = 150, verbose = FALSE,
              data = DT)
summary(fitRR)$varcomp
##                           term                  parameter    estimate  StdError
## 1  vsm(rrm(Env, 2), ism(Name))        loading[CA.2011,F1]  3.51395448 0.5290545
## 2  vsm(rrm(Env, 2), ism(Name))        loading[CA.2012,F1]  1.55652058 0.3325198
## 3  vsm(rrm(Env, 2), ism(Name))        loading[CA.2013,F1]  1.85826420 0.4618439
## 4  vsm(rrm(Env, 2), ism(Name))        loading[FL.2011,F1]  0.60225228 0.4621690
## 5  vsm(rrm(Env, 2), ism(Name))        loading[FL.2012,F1]  0.26384843 0.3310584
## 6  vsm(rrm(Env, 2), ism(Name))        loading[FL.2013,F1]  0.73549270 0.4278361
## 7  vsm(rrm(Env, 2), ism(Name))        loading[MI.2011,F1]  2.50303356 0.5757270
## 8  vsm(rrm(Env, 2), ism(Name))        loading[MI.2012,F1]  1.55771594 0.3460649
## 9  vsm(rrm(Env, 2), ism(Name))        loading[MI.2013,F1]  3.78708101 0.5949514
## 10 vsm(rrm(Env, 2), ism(Name))        loading[MO.2011,F1]  0.87430884 0.4771665
## 11 vsm(rrm(Env, 2), ism(Name))        loading[MO.2012,F1]  3.07358283 0.7202287
## 12 vsm(rrm(Env, 2), ism(Name))        loading[MO.2013,F1]  1.28460775 0.3699165
## 13 vsm(rrm(Env, 2), ism(Name))        loading[NY.2011,F1]  1.88915895 0.6937754
## 14 vsm(rrm(Env, 2), ism(Name))        loading[NY.2012,F1]  1.25485884 0.3419116
## 15 vsm(rrm(Env, 2), ism(Name))        loading[NY.2013,F1]  1.53407576 0.4674285
## 16 vsm(rrm(Env, 2), ism(Name))        loading[CA.2012,F2]  0.45186512 0.3585174
## 17 vsm(rrm(Env, 2), ism(Name))        loading[CA.2013,F2] -1.48782374 0.4416179
## 18 vsm(rrm(Env, 2), ism(Name))        loading[FL.2011,F2]  0.83690449 0.4748106
## 19 vsm(rrm(Env, 2), ism(Name))        loading[FL.2012,F2]  0.29585454 0.3601858
## 20 vsm(rrm(Env, 2), ism(Name))        loading[FL.2013,F2]  0.35171587 0.5593048
## 21 vsm(rrm(Env, 2), ism(Name))        loading[MI.2011,F2]  1.38385183 0.4646787
## 22 vsm(rrm(Env, 2), ism(Name))        loading[MI.2012,F2]  0.65446597 0.3244937
## 23 vsm(rrm(Env, 2), ism(Name))        loading[MI.2013,F2] -1.37340771 0.5156478
## 24 vsm(rrm(Env, 2), ism(Name))        loading[MO.2011,F2] -0.07736326 0.4956230
## 25 vsm(rrm(Env, 2), ism(Name))        loading[MO.2012,F2]  2.65048421 0.5807362
## 26 vsm(rrm(Env, 2), ism(Name))        loading[MO.2013,F2] -0.10370748 0.4210157
## 27 vsm(rrm(Env, 2), ism(Name))        loading[NY.2011,F2]  1.94601737 0.5683801
## 28 vsm(rrm(Env, 2), ism(Name))        loading[NY.2012,F2]  0.34548312 0.3623471
## 29 vsm(rrm(Env, 2), ism(Name))        loading[NY.2013,F2]  1.60781749 0.4002131
## 30 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2011]  0.00000000 0.0000000
## 31 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2012]  0.00000000 0.0000000
## 32 vsm(rrm(Env, 2), ism(Name)) specific_variance[CA.2013]  0.00000000 0.0000000
## 33 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2011]  0.00000000 0.0000000
## 34 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2012]  0.00000000 0.0000000
## 35 vsm(rrm(Env, 2), ism(Name)) specific_variance[FL.2013]  0.00000000 0.0000000
## 36 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2011]  0.00000000 0.0000000
## 37 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2012]  0.00000000 0.0000000
## 38 vsm(rrm(Env, 2), ism(Name)) specific_variance[MI.2013]  0.00000000 0.0000000
## 39 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2011]  0.00000000 0.0000000
## 40 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2012]  0.00000000 0.0000000
## 41 vsm(rrm(Env, 2), ism(Name)) specific_variance[MO.2013]  0.00000000 0.0000000
## 42 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2011]  0.00000000 0.0000000
## 43 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2012]  0.00000000 0.0000000
## 44 vsm(rrm(Env, 2), ism(Name)) specific_variance[NY.2013]  0.00000000 0.0000000
## 45             vsm(ism(units))                     sigma2  5.76619309 0.2930191
##        Zratio
## 1   6.6419522
## 2   4.6809856
## 3   4.0235763
## 4   1.3030996
## 5   0.7969845
## 6   1.7190992
## 7   4.3476050
## 8   4.5012247
## 9   6.3653618
## 10  1.8322929
## 11  4.2675097
## 12  3.4726964
## 13  2.7230122
## 14  3.6701266
## 15  3.2819472
## 16  1.2603714
## 17 -3.3690295
## 18  1.7626070
## 19  0.8213942
## 20  0.6288448
## 21  2.9780831
## 22  2.0168833
## 23 -2.6634609
## 24 -0.1560929
## 25  4.5640069
## 26 -0.2463268
## 27  3.4237956
## 28  0.9534591
## 29  4.0174034
## 30         NA
## 31         NA
## 32         NA
## 33         NA
## 34         NA
## 35         NA
## 36         NA
## 37         NA
## 38         NA
## 39         NA
## 40         NA
## 41         NA
## 42         NA
## 43         NA
## 44         NA
## 45 19.6785567

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   66.479270 
##  logLik_FA7  logLik_FA8  logLik_FA9 logLik_FA10 logLik_FA11 logLik_FA12 
##   67.326222   67.643139   67.760623   67.803979   79.269472   82.519873 
## logLik_FA13 logLik_FA14 logLik_FA15 logLik_FA16 logLik_FA17 logLik_FA18 
##   83.697757   84.301934   84.670037   84.928786   85.430834   85.992197 
## logLik_FA19 logLik_FA20 logLik_FA21 logLik_FA22 logLik_FA23 logLik_FA24 
##   86.281680   86.463984   86.535274   86.805749   86.915006   86.984039 
## logLik_FA25 logLik_FA26 logLik_FA27 logLik_FA28 logLik_FA29 logLik_FA30 
##   87.022103   87.029279   87.121833   87.147329   87.165486   87.234728 
## logLik_FA31 logLik_FA32 logLik_FA33 logLik_FA34 logLik_FA35 logLik_FA36 
##   87.249526   87.269693   87.277306   87.281706   87.284290   87.285567 
## logLik_FA37 logLik_FA38 logLik_FA39 logLik_FA40 logLik_FA41 logLik_FA42 
##   87.285682   87.290406   87.291552   87.297759   87.299271   87.299883 
## logLik_FA43 logLik_FA44 logLik_FA45 logLik_FA46 logLik_FA47 logLik_FA48 
##   87.300029   87.301129   87.301419   87.301624   87.301764   87.301844 
## logLik_FA49 logLik_FA50 logLik_FA51 logLik_FA52 logLik_FA53 logLik_FA54 
##   87.302532   87.302880   87.303000   87.303062   87.303096   87.303110 
## logLik_FA55 logLik_FA56 logLik_FA57 logLik_FA58 logLik_FA59 logLik_FA60 
##   87.303120   87.303122   87.303188   87.303204   87.303205   87.303218 
## logLik_FA61 logLik_FA62  logLik_RR1  logLik_RR2  logLik_RR3  logLik_RR4 
##   87.303227   87.303234 -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 logLik_RR31 logLik_RR32 logLik_RR33 logLik_RR34 
##   35.233175   35.241658   35.249137   35.255734   35.261540   35.266649 
## logLik_RR35 logLik_RR36 logLik_RR37 logLik_RR38 logLik_RR39 logLik_RR40 
##   35.271132   35.275057   35.278485   35.281462   35.284038   35.286243 
## logLik_RR41 logLik_RR42 logLik_RR43 logLik_RR44 logLik_RR45 logLik_RR46 
##   35.288122   35.289686   35.290980   35.291999   35.292787   35.293322 
## logLik_RR47 logLik_RR48 logLik_RR49 logLik_RR50 logLik_RR51 logLik_RR52 
##   35.293653   35.293739   35.302627   35.303451   35.304178   35.304822 
## logLik_RR53 logLik_RR54 logLik_RR55 logLik_RR56 logLik_RR57 logLik_RR58 
##   35.305392   35.305898   35.306346   35.306743   35.307094   35.307405 
## logLik_RR59 logLik_RR60 logLik_RR61 logLik_RR62 logLik_RR63 logLik_RR64 
##   35.307681   35.307924   35.308139   35.308330   35.308498   35.308646 
## logLik_RR65 logLik_RR66 logLik_RR67 logLik_RR68 
##   35.308776   35.308890   35.308991   35.309078

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.348
## CA.2012 0.357  0.136
## CA.2013 0.376  0.358
## FL.2011 0.236 -0.218
## FL.2012 0.104 -0.104
## FL.2013 0.231  0.088
## MI.2011 0.613 -0.011
## MI.2012 0.417 -0.012
## MI.2013 0.714  0.104
## MO.2011 0.159  0.152
## MO.2012 0.687 -0.063
## MO.2013 0.339 -0.007
## NY.2011 0.459 -0.348
## NY.2012 0.359  0.058
## NY.2013 0.526 -0.526
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.426   0.194   0.231   0.003   0.000   0.000   0.194   0.147   0.718   0.000 
## MO.2012 MO.2013 NY.2011 NY.2012 NY.2013 
##   0.574   0.002   0.000   0.000   0.024
faInfo$sigma2
## [1] 15.69948
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.183497889  1.3795638
## AC00206-2W  0.008665097 -0.3382699
## AC01151-5W  0.548157714 -0.1036613
## AC03433-1W -1.025464860 -0.1048459
## AC03452-2W  2.001615758 -0.1857154
## AC05153-1W -1.175652952 -0.3726147

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 = FALSE retains 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.138  0.000
## F2  0.000 12.138

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.