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
- Percentage of variance explained
- Loadings plot
- Genotype scores biplot
- 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 the environment-specific remainder to be homogeneous:
$$ K = \Lambda\Lambda^{\mathsf T} + I. $$
rrm(x, k) therefore has
$$ kq-\frac{k(k-1)}{2} $$
working parameters: exactly \(q-1\) fewer than fam(x, k), because it does not estimate separate environment-specific variances. RR is a restricted (nested) special case of FA in which every environment is assumed to have the same residual/specific variance after accounting for the \(k\) common factors. It is useful when there is not enough replication to support \(q\) separate specific variances, or as a more parsimonious first approximation before fitting the full FA model.
4) Choosing k and comparing to other structures
Increasing \(k\) moves the model continuously from compound symmetry-like behavior (\(k=0\), not directly supported, but conceptually the limit) toward the saturated unstructured model (\(k=q-1\)). In practice \(k\) is chosen small enough to remain identifiable and interpretable (often 1-3 factors), and increased only if it meaningfully improves the likelihood. Because RR is nested inside FA with the same \(k\), anova.mmes() can be used to test whether the extra \(q-1\) FA specific-variance parameters are supported by the data.
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)
## Loading required package: Matrix
## Loading required package: MASS
## Loading required package: crayon
## Loading required package: enhancer
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)
## Solver selected: ldlt
summary(fitFA)$varcomp
## term parameter estimate
## 1 vsm(fam(Env, 2), ism(Name)) loading[CA.2011,F1] 3.0048776208
## 2 vsm(fam(Env, 2), ism(Name)) loading[CA.2012,F1] 1.5069866827
## 3 vsm(fam(Env, 2), ism(Name)) loading[CA.2013,F1] 1.9780068616
## 4 vsm(fam(Env, 2), ism(Name)) loading[FL.2011,F1] 0.4320595180
## 5 vsm(fam(Env, 2), ism(Name)) loading[FL.2012,F1] 0.1769448480
## 6 vsm(fam(Env, 2), ism(Name)) loading[FL.2013,F1] 0.9697303822
## 7 vsm(fam(Env, 2), ism(Name)) loading[MI.2011,F1] 2.1365415969
## 8 vsm(fam(Env, 2), ism(Name)) loading[MI.2012,F1] 1.4410923427
## 9 vsm(fam(Env, 2), ism(Name)) loading[MI.2013,F1] 2.6954117697
## 10 vsm(fam(Env, 2), ism(Name)) loading[MO.2011,F1] 0.8366710461
## 11 vsm(fam(Env, 2), ism(Name)) loading[MO.2012,F1] 2.2985581637
## 12 vsm(fam(Env, 2), ism(Name)) loading[MO.2013,F1] 1.1770823749
## 13 vsm(fam(Env, 2), ism(Name)) loading[NY.2011,F1] 0.9772044183
## 14 vsm(fam(Env, 2), ism(Name)) loading[NY.2012,F1] 1.3696291976
## 15 vsm(fam(Env, 2), ism(Name)) loading[NY.2013,F1] 0.8890274338
## 16 vsm(fam(Env, 2), ism(Name)) loading[CA.2012,F2] 0.1722044759
## 17 vsm(fam(Env, 2), ism(Name)) loading[CA.2013,F2] -0.5769071594
## 18 vsm(fam(Env, 2), ism(Name)) loading[FL.2011,F2] 1.1970074325
## 19 vsm(fam(Env, 2), ism(Name)) loading[FL.2012,F2] 0.5512164979
## 20 vsm(fam(Env, 2), ism(Name)) loading[FL.2013,F2] 0.1170903499
## 21 vsm(fam(Env, 2), ism(Name)) loading[MI.2011,F2] 1.1594271596
## 22 vsm(fam(Env, 2), ism(Name)) loading[MI.2012,F2] 0.8073523687
## 23 vsm(fam(Env, 2), ism(Name)) loading[MI.2013,F2] 0.9449444436
## 24 vsm(fam(Env, 2), ism(Name)) loading[MO.2011,F2] -0.2454739294
## 25 vsm(fam(Env, 2), ism(Name)) loading[MO.2012,F2] 1.4846023009
## 26 vsm(fam(Env, 2), ism(Name)) loading[MO.2013,F2] 0.6438739193
## 27 vsm(fam(Env, 2), ism(Name)) loading[NY.2011,F2] 2.0645100911
## 28 vsm(fam(Env, 2), ism(Name)) loading[NY.2012,F2] 0.4523827527
## 29 vsm(fam(Env, 2), ism(Name)) loading[NY.2013,F2] 2.8078628941
## 30 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2011] 6.6633457278
## 31 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2012] 3.0463630516
## 32 vsm(fam(Env, 2), ism(Name)) specific_variance[CA.2013] 3.6090565572
## 33 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2011] 0.0501792853
## 34 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2012] 0.0001869543
## 35 vsm(fam(Env, 2), ism(Name)) specific_variance[FL.2013] 0.0001585274
## 36 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2011] 3.0528782578
## 37 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2012] 2.3169108493
## 38 vsm(fam(Env, 2), ism(Name)) specific_variance[MI.2013] 11.2727450504
## 39 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2011] 0.0001995198
## 40 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2012] 9.0139439913
## 41 vsm(fam(Env, 2), ism(Name)) specific_variance[MO.2013] 0.0271409700
## 42 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2011] 0.0002702601
## 43 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2012] 0.0003522949
## 44 vsm(fam(Env, 2), ism(Name)) specific_variance[NY.2013] 0.3846473158
## 45 vsm(ism(units)) sigma2 3.9688245476
## StdError Zratio
## 1 0.8220594 3.6553045198
## 2 0.4635645 3.2508670004
## 3 0.5739203 3.4464837813
## 4 0.5192200 0.8321319261
## 5 0.3261388 0.5425445689
## 6 0.3677144 2.6371834410
## 7 0.7774825 2.7480254329
## 8 0.4946987 2.9130709233
## 9 0.9137577 2.9498103398
## 10 0.3921379 2.1336141187
## 11 0.8826952 2.6040225347
## 12 0.3734346 3.1520442450
## 13 0.7887810 1.2388792069
## 14 0.3416226 4.0091879156
## 15 0.9711351 0.9154518344
## 16 0.6429618 0.2678300440
## 17 0.8191669 -0.7042608519
## 18 0.4728781 2.5313232760
## 19 0.2963932 1.8597475191
## 20 0.5016201 0.2334243484
## 21 0.8156593 1.4214600361
## 22 0.5440171 1.4840569610
## 23 1.1337752 0.8334495705
## 24 0.4912119 -0.4997312657
## 25 1.0037833 1.4790067647
## 26 0.4670677 1.3785452271
## 27 0.5967803 3.4594138034
## 28 0.4764476 0.9494911077
## 29 0.6534564 4.2969400219
## 30 3.4890290 1.9097994786
## 31 1.3314603 2.2879864725
## 32 1.8765349 1.9232557874
## 33 0.7823119 0.0641423020
## 34 0.5338621 0.0003501921
## 35 0.7536168 0.0002103555
## 36 2.0093099 1.5193665644
## 37 1.1006971 2.1049486168
## 38 3.8942898 2.8946857132
## 39 0.9373965 0.0002128446
## 40 3.2886427 2.7409313930
## 41 0.5877766 0.0461756522
## 42 1.2768260 0.0002116656
## 43 0.5766149 0.0006109710
## 44 2.0163236 0.1907666569
## 45 0.2652474 14.9627247529
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)
## Solver selected: ldlt
summary(fitRR)$varcomp
## term parameter estimate StdError
## 1 vsm(rrm(Env, 2), ism(Name)) loading[CA.2011,F1] 3.38934720 0.5734663
## 2 vsm(rrm(Env, 2), ism(Name)) loading[CA.2012,F1] 1.50785453 0.3882794
## 3 vsm(rrm(Env, 2), ism(Name)) loading[CA.2013,F1] 1.73866657 0.4968970
## 4 vsm(rrm(Env, 2), ism(Name)) loading[FL.2011,F1] 0.59655414 0.5456171
## 5 vsm(rrm(Env, 2), ism(Name)) loading[FL.2012,F1] 0.27676260 0.4045255
## 6 vsm(rrm(Env, 2), ism(Name)) loading[FL.2013,F1] 0.72260685 0.5130813
## 7 vsm(rrm(Env, 2), ism(Name)) loading[MI.2011,F1] 2.53246286 0.6444075
## 8 vsm(rrm(Env, 2), ism(Name)) loading[MI.2012,F1] 1.58460992 0.3991016
## 9 vsm(rrm(Env, 2), ism(Name)) loading[MI.2013,F1] 3.59836849 0.6445868
## 10 vsm(rrm(Env, 2), ism(Name)) loading[MO.2011,F1] 0.84308408 0.5670675
## 11 vsm(rrm(Env, 2), ism(Name)) loading[MO.2012,F1] 3.00524303 0.7631330
## 12 vsm(rrm(Env, 2), ism(Name)) loading[MO.2013,F1] 1.25710408 0.4418251
## 13 vsm(rrm(Env, 2), ism(Name)) loading[NY.2011,F1] 1.82205728 0.7628243
## 14 vsm(rrm(Env, 2), ism(Name)) loading[NY.2012,F1] 1.24861155 0.4086312
## 15 vsm(rrm(Env, 2), ism(Name)) loading[NY.2013,F1] 1.60427308 0.4921529
## 16 vsm(rrm(Env, 2), ism(Name)) loading[CA.2012,F2] 0.48668487 0.4619068
## 17 vsm(rrm(Env, 2), ism(Name)) loading[CA.2013,F2] -1.32753190 0.5251877
## 18 vsm(rrm(Env, 2), ism(Name)) loading[FL.2011,F2] 0.74601232 0.5814341
## 19 vsm(rrm(Env, 2), ism(Name)) loading[FL.2012,F2] 0.24127276 0.4656463
## 20 vsm(rrm(Env, 2), ism(Name)) loading[FL.2013,F2] 0.37541144 0.7062174
## 21 vsm(rrm(Env, 2), ism(Name)) loading[MI.2011,F2] 1.12552899 0.5800393
## 22 vsm(rrm(Env, 2), ism(Name)) loading[MI.2012,F2] 0.48288018 0.4247066
## 23 vsm(rrm(Env, 2), ism(Name)) loading[MI.2013,F2] -1.50447942 0.6400994
## 24 vsm(rrm(Env, 2), ism(Name)) loading[MO.2011,F2] -0.04044858 0.6087592
## 25 vsm(rrm(Env, 2), ism(Name)) loading[MO.2012,F2] 2.50244983 0.6653206
## 26 vsm(rrm(Env, 2), ism(Name)) loading[MO.2013,F2] -0.15701037 0.5290051
## 27 vsm(rrm(Env, 2), ism(Name)) loading[NY.2011,F2] 1.66698281 0.6550612
## 28 vsm(rrm(Env, 2), ism(Name)) loading[NY.2012,F2] 0.30101764 0.4700564
## 29 vsm(rrm(Env, 2), ism(Name)) loading[NY.2013,F2] 1.29832973 0.4804319
## 30 vsm(rrm(Env, 2), ism(Name)) common_specific_variance 1.76775123 0.3366846
## 31 vsm(ism(units)) sigma2 4.24057879 0.2923669
## Zratio
## 1 5.9102810
## 2 3.8834265
## 3 3.4990485
## 4 1.0933568
## 5 0.6841660
## 6 1.4083670
## 7 3.9299089
## 8 3.9704429
## 9 5.5824423
## 10 1.4867438
## 11 3.9380331
## 12 2.8452528
## 13 2.3885673
## 14 3.0555954
## 15 3.2597043
## 16 1.0536429
## 17 -2.5277282
## 18 1.2830557
## 19 0.5181460
## 20 0.5315806
## 21 1.9404359
## 22 1.1369735
## 23 -2.3503841
## 24 -0.0664443
## 25 3.7612689
## 26 -0.2968031
## 27 2.5447743
## 28 0.6403862
## 29 2.7024221
## 30 5.2504666
## 31 14.5043034
4) Comparing FA and RR
Because rrm(Env, 2) is nested inside fam(Env, 2) (both use rank 2, but RR fixes the specific variances to be equal), we can compare them with a likelihood ratio test:
c(AIC_FA = fitFA$AIC, AIC_RR = fitRR$AIC)
## AIC_FA AIC_RR
## -144.60344 -72.58723
anova.mmes(fitFA, fitRR)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 60 -144.60344 -73.041603 87.30172
## mod2 46 -72.58723 -1.025386 51.29361 72.01622 14 8.3066901414204e-10 ***
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 60 -144.60344 -73.041603 87.30172
## mod2 46 -72.58723 -1.025386 51.29361 72.01622 14 8.3066901414204e-10 ***
A lower AIC/BIC and a significant likelihood ratio test favor the less restrictive FA model whenever the data support heterogeneous environment-specific variances.
SECTION 3: Extracting loadings, scores, and diagnostic plots
1) Loadings and specific variances
loadings_mmes() reconstructs the fitted, normalized loadings \(\Lambda\) and specific variances \(\Psi\) so that \(\sigma^2(\Lambda\Lambda^{\mathsf T}+\Psi)\) reproduces the fitted covariance matrix exactly.
faInfo <- loadings_mmes(fitFA)
round(faInfo$loadings, 3)
## F1 F2
## CA.2011 2.666 1.386
## CA.2012 1.416 0.542
## CA.2013 1.489 1.424
## FL.2011 0.936 -0.863
## FL.2012 0.411 -0.407
## FL.2013 0.914 0.343
## MI.2011 2.430 -0.043
## MI.2012 1.651 -0.052
## MI.2013 2.827 0.405
## MO.2011 0.629 0.604
## MO.2012 2.724 -0.257
## MO.2013 1.341 -0.028
## NY.2011 1.819 -1.381
## NY.2012 1.424 0.230
## NY.2013 2.084 -2.081
round(faInfo$specific, 3)
## CA.2011 CA.2012 CA.2013 FL.2011 FL.2012 FL.2013 MI.2011 MI.2012 MI.2013 MO.2011
## 1.682 0.769 0.911 0.013 0.000 0.000 0.771 0.585 2.846 0.000
## MO.2012 MO.2013 NY.2011 NY.2012 NY.2013
## 2.275 0.007 0.000 0.000 0.097
faInfo$sigma2
## [1] 15.69264
rrInfo <- loadings_mmes(fitRR)
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.486 0.486 0.486 0.486 0.486 0.486 0.486 0.486 0.486 0.486
## MO.2012 MO.2013 NY.2011 NY.2012 NY.2013
## 0.486 0.486 0.486 0.486 0.486
Notice that rrInfo$specific is constant across environments: this is the direct consequence of rrm() fixing a homogeneous specific variance, unlike the heterogeneous values in faInfo$specific.
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.
faScores <- scores_mmes(fitFA)
head(faScores)
## F1 F2
## A01143-3C 1.183505909 1.37997480
## AC00206-2W 0.008557919 -0.33899056
## AC01151-5W 0.548555060 -0.09987134
## AC03433-1W -1.025464223 -0.10407867
## AC03452-2W 2.001549121 -0.18598226
## AC05153-1W -1.175657607 -0.37370696
3) 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
## 69.0 17.1
4) 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")
5) Genotype scores biplot
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)
6) 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)
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).
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.