The sommer package was developed to provide R users with a powerful and reliable multivariate mixed model solver for different genetic and non-genetic analyses in diploid and polyploid organisms. This package allows the user to estimate variance components for a mixed model with the advantages of specifying the variance-covariance structure of the random effects, specifying heterogeneous variances, and obtaining other parameters such as BLUPs, BLUEs, residuals, fitted values, variances for fixed and random effects, etc. The core algorithms of the package are coded in C++ using the Armadillo library to optimize dense matrix operations common in the derect-inversion algorithms. Although the vignette shows examples using the mmes function with the default direct inversion algorithm (henderson=FALSE) the Henderson’s approach can be faster when the number of records surpasses the number of coefficients to estimate and setting the henderson argument to TRUE can bring significant speed ups.
The purpose of this vignette is to show how to translate the syntax formula from lme4 models to sommer models. Feel free to remove the comment marks from the lme4 code so you can compare the results.
This is the simplest model people use when a random effect is desired and the levels of the random effect are considered to have the same intercept.
# install.packages("lme4")
# library(lme4)
library(sommer)
data(DT_sleepstudy, package="enhancer")
DT <- DT_sleepstudy
###########
## lme4
###########
# fm1 <- lmer(Reaction ~ Days + (1 | Subject), data=DT)
# summary(fm1) # or vc <- VarCorr(fm1); print(vc,comp=c("Variance"))
# Random effects:
# Groups Name Variance Std.Dev.
# Subject (Intercept) 1378.2 37.12
# Residual 960.5 30.99
# Number of obs: 180, groups: Subject, 18
###########
## sommer
###########
fm2 <- mmes(Reaction ~ Days,
random= ~ Subject,
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(ism(Subject)) sigma2 1377.9699 505.6628 2.725077
## 2 vsm(ism(units)) sigma2 960.4965 107.0550 8.971993
This is the a model where you assume that the random effect has different intercepts based on the levels of another variable. In addition the || in lme4 assumes that slopes and intercepts have no correlation.
###########
## lme4
###########
# fm1 <- lmer(Reaction ~ Days + (Days || Subject), data=DT)
# summary(fm1) # or vc <- VarCorr(fm1); print(vc,comp=c("Variance"))
# Random effects:
# Groups Name Variance Std.Dev.
# Subject (Intercept) 627.57 25.051
# Subject.1 Days 35.86 5.988
# Residual 653.58 25.565
# Number of obs: 180, groups: Subject, 18
###########
## sommer
###########
fm2 <- mmes(Reaction ~ Days,
random= ~ Subject + vsm(ism(Days), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(ism(Subject)) sigma2 627.47974 283.61898 2.212404
## 2 vsm(ism(Days), ism(Subject)) sigma2 35.86523 14.53347 2.467768
## 3 vsm(ism(units)) sigma2 653.72685 76.68247 8.525115
Notice that Days is a numerical (not factor) variable.
This is the a model where you assume that the random effect has different intercepts based on the levels of another variable. In addition a single | in lme4 assumes that slopes and intercepts have a correlation to be estimated.
###########
## lme4
###########
# fm1 <- lmer(Reaction ~ Days + (Days | Subject), data=DT)
# summary(fm1) # or # vc <- VarCorr(fm1); print(vc,comp=c("Variance"))
# Random effects:
# Groups Name Variance Std.Dev. Corr
# Subject (Intercept) 612.10 24.741
# Days 35.07 5.922 0.07
# Residual 654.94 25.592
# Number of obs: 180, groups: Subject, 18
###########
## sommer
###########
fm2 <- mmes(Reaction ~ Days, # henderson=TRUE,
random= ~ vsm(usm(cbind(1,Days)), ism(Subject)) ,
nIters = 200, data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError
## 1 vsm(usm(cbind(1, Days)), ism(Subject)) variance[] 611.679210 288.57422
## 2 vsm(usm(cbind(1, Days)), ism(Subject)) covariance[Days,] 9.615479 46.66683
## 3 vsm(usm(cbind(1, Days)), ism(Subject)) variance[Days] 35.078371 14.78558
## 4 vsm(ism(units)) sigma2 654.951742 77.18743
## Zratio
## 1 2.1196599
## 2 0.2060453
## 3 2.3724724
## 4 8.4852120
cov2cor(fm2$theta[[1]])
## [,1] [,2]
## [1,] 1.00000000 0.06564314
## [2,] 0.06564314 1.00000000
Notice that this last model require a new function called covm() which creates the two random effects as before but now they have to be encapsulated in covm() instead of just added.
This is the a model where you assume that the random effect has different intercepts based on the levels of another variable but there’s not a main effect. The 0 in the intercept in lme4 assumes that random slopes interact with an intercept but without a main effect.
###########
## lme4
###########
# fm1 <- lmer(Reaction ~ Days + (0 + Days | Subject), data=DT)
# summary(fm1) # or vc <- VarCorr(fm1); print(vc,comp=c("Variance"))
# Random effects:
# Groups Name Variance Std.Dev.
# Subject Days 52.71 7.26
# Residual 842.03 29.02
# Number of obs: 180, groups: Subject, 18
###########
## sommer
###########
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(ism(Days), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(ism(Days), ism(Subject)) sigma2 52.71141 19.09684 2.760217
## 2 vsm(ism(units)) sigma2 842.12287 93.86454 8.971683
One of the strengths of sommer is the availability of other variance covariance structures. In this section we show 4 models available in sommer that are not available in lme4 and might be useful.
library(orthopolynom)
## diagonal model
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(dsm(Daysf), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(dsm(Daysf), ism(Subject)) variance[0] 502.6988 74.47413 6.7499788
## 2 vsm(dsm(Daysf), ism(Subject)) variance[1] 559.2660 503.31690 1.1111608
## 3 vsm(dsm(Daysf), ism(Subject)) variance[2] 362.3340 439.69453 0.8240585
## 4 vsm(dsm(Daysf), ism(Subject)) variance[3] 918.4486 607.42693 1.5120314
## 5 vsm(dsm(Daysf), ism(Subject)) variance[4] 1217.5541 695.56641 1.7504499
## 6 vsm(dsm(Daysf), ism(Subject)) variance[5] 2061.7874 953.75270 2.1617631
## 7 vsm(dsm(Daysf), ism(Subject)) variance[6] 3273.2376 1339.72125 2.4432229
## 8 vsm(dsm(Daysf), ism(Subject)) variance[7] 1901.9830 906.29992 2.0986243
## 9 vsm(dsm(Daysf), ism(Subject)) variance[8] 2959.8719 1243.30782 2.3806429
## 10 vsm(dsm(Daysf), ism(Subject)) variance[9] 3835.0209 1530.77404 2.5052822
## 11 vsm(ism(units)) sigma2 517.2536 345.32710 1.4978656
## unstructured model
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(usm(Daysf), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(usm(Daysf), ism(Subject)) variance[0] 405.57603 112.77359 3.596374
## 2 vsm(usm(Daysf), ism(Subject)) covariance[1,0] 419.54898 149.60096 2.804454
## 3 vsm(usm(Daysf), ism(Subject)) covariance[2,0] 269.94635 154.77620 1.744108
## 4 vsm(usm(Daysf), ism(Subject)) covariance[3,0] 396.56888 224.89559 1.763347
## 5 vsm(usm(Daysf), ism(Subject)) covariance[4,0] 437.22604 260.40675 1.679012
## 6 vsm(usm(Daysf), ism(Subject)) covariance[5,0] 519.16370 331.44832 1.566349
## 7 vsm(usm(Daysf), ism(Subject)) covariance[6,0] 444.15322 378.45746 1.173588
## 8 vsm(usm(Daysf), ism(Subject)) covariance[7,0] 531.80093 308.08739 1.726137
## 9 vsm(usm(Daysf), ism(Subject)) covariance[8,0] 571.59031 406.90707 1.404720
## 10 vsm(usm(Daysf), ism(Subject)) covariance[9,0] 774.18737 411.63528 1.880761
## 11 vsm(usm(Daysf), ism(Subject)) variance[1] 842.93172 286.50404 2.942129
## 12 vsm(usm(Daysf), ism(Subject)) covariance[2,1] 753.02600 254.71659 2.956329
## 13 vsm(usm(Daysf), ism(Subject)) covariance[3,1] 1014.18375 363.07725 2.793300
## 14 vsm(usm(Daysf), ism(Subject)) covariance[4,1] 1011.38456 380.20322 2.660116
## 15 vsm(usm(Daysf), ism(Subject)) covariance[5,1] 1106.77204 497.95570 2.222632
## 16 vsm(usm(Daysf), ism(Subject)) covariance[6,1] 976.74160 540.36119 1.807572
## 17 vsm(usm(Daysf), ism(Subject)) covariance[7,1] 889.23135 449.88349 1.976581
## 18 vsm(usm(Daysf), ism(Subject)) covariance[8,1] 1116.15712 595.53479 1.874210
## 19 vsm(usm(Daysf), ism(Subject)) covariance[9,1] 1320.83826 638.38385 2.069035
## 20 vsm(usm(Daysf), ism(Subject)) variance[2] 978.22845 304.93364 3.208004
## 21 vsm(usm(Daysf), ism(Subject)) covariance[3,2] 1253.82005 401.35567 3.123962
## 22 vsm(usm(Daysf), ism(Subject)) covariance[4,2] 1279.81144 433.98172 2.948999
## 23 vsm(usm(Daysf), ism(Subject)) covariance[5,2] 1289.74905 560.51727 2.300998
## 24 vsm(usm(Daysf), ism(Subject)) covariance[6,2] 1490.62884 664.18520 2.244297
## 25 vsm(usm(Daysf), ism(Subject)) covariance[7,2] 1319.86722 551.93759 2.391334
## 26 vsm(usm(Daysf), ism(Subject)) covariance[8,2] 1468.18756 705.19623 2.081956
## 27 vsm(usm(Daysf), ism(Subject)) covariance[9,2] 1403.69388 712.59805 1.969826
## 28 vsm(usm(Daysf), ism(Subject)) variance[3] 1916.91118 673.49316 2.846222
## 29 vsm(usm(Daysf), ism(Subject)) covariance[4,3] 2117.49170 750.70371 2.820676
## 30 vsm(usm(Daysf), ism(Subject)) covariance[5,3] 2331.60879 960.06810 2.428587
## 31 vsm(usm(Daysf), ism(Subject)) covariance[6,3] 2665.24651 1121.75369 2.375964
## 32 vsm(usm(Daysf), ism(Subject)) covariance[7,3] 1899.64099 852.86867 2.227355
## 33 vsm(usm(Daysf), ism(Subject)) covariance[8,3] 2552.46704 1165.80535 2.189445
## 34 vsm(usm(Daysf), ism(Subject)) covariance[9,3] 2438.70404 1156.17321 2.109290
## 35 vsm(usm(Daysf), ism(Subject)) variance[4] 2536.35767 931.08330 2.724093
## 36 vsm(usm(Daysf), ism(Subject)) covariance[5,4] 2968.06557 1191.87407 2.490251
## 37 vsm(usm(Daysf), ism(Subject)) covariance[6,4] 3283.10388 1378.42518 2.381779
## 38 vsm(usm(Daysf), ism(Subject)) covariance[7,4] 2432.16306 1063.19059 2.287608
## 39 vsm(usm(Daysf), ism(Subject)) covariance[8,4] 3356.86980 1458.41214 2.301729
## 40 vsm(usm(Daysf), ism(Subject)) covariance[9,4] 3250.12813 1442.77951 2.252685
## 41 vsm(usm(Daysf), ism(Subject)) variance[5] 4132.83284 1677.97235 2.462992
## 42 vsm(usm(Daysf), ism(Subject)) covariance[6,5] 4140.83371 1841.98376 2.248029
## 43 vsm(usm(Daysf), ism(Subject)) covariance[7,5] 3130.01495 1448.87371 2.160309
## 44 vsm(usm(Daysf), ism(Subject)) covariance[8,5] 4824.14737 2064.51791 2.336694
## 45 vsm(usm(Daysf), ism(Subject)) covariance[9,5] 4671.98066 2025.63569 2.306427
## 46 vsm(usm(Daysf), ism(Subject)) variance[6] 5843.09484 2355.66296 2.480446
## 47 vsm(usm(Daysf), ism(Subject)) covariance[7,6] 3820.45430 1733.85978 2.203439
## 48 vsm(usm(Daysf), ism(Subject)) covariance[8,6] 5083.37605 2341.09439 2.171367
## 49 vsm(usm(Daysf), ism(Subject)) covariance[9,6] 4028.45638 2159.88736 1.865123
## 50 vsm(usm(Daysf), ism(Subject)) variance[7] 3602.47750 1494.52761 2.410446
## 51 vsm(usm(Daysf), ism(Subject)) covariance[8,7] 4163.51193 1903.74792 2.187008
## 52 vsm(usm(Daysf), ism(Subject)) covariance[9,7] 3749.93743 1814.18180 2.067013
## 53 vsm(usm(Daysf), ism(Subject)) variance[8] 6262.42482 2700.09447 2.319335
## 54 vsm(usm(Daysf), ism(Subject)) covariance[9,8] 5978.41689 2613.24451 2.287737
## 55 vsm(usm(Daysf), ism(Subject)) variance[9] 6404.46192 2665.80185 2.402452
## 56 vsm(ism(units)) sigma2 68.32424 30.65416 2.228873
## random regression (legendre polynomials)
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(dsm(leg(Days,1)), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate StdError
## 1 vsm(dsm(leg(Days, 1)), ism(Subject)) variance[leg0] 2815.6001 1010.47021
## 2 vsm(dsm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.6353 199.64378
## 3 vsm(ism(units)) sigma2 654.9382 77.18503
## Zratio
## 1 2.786426
## 2 2.372402
## 3 8.485300
## unstructured random regression (legendre)
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(usm(leg(Days,1)), ism(Subject)),
data=DT, tolParInv = 1e-6, verbose = FALSE)
## Solver selected: ldlt
summary(fm2)$varcomp
## term parameter estimate
## 1 vsm(usm(leg(Days, 1)), ism(Subject)) variance[leg0] 2816.8784
## 2 vsm(usm(leg(Days, 1)), ism(Subject)) covariance[leg1,leg0] 869.8019
## 3 vsm(usm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.4508
## 4 vsm(ism(units)) sigma2 654.9370
## StdError Zratio
## 1 1011.12756 2.785878
## 2 381.00537 2.282913
## 3 199.54894 2.372605
## 4 77.18482 8.485308
Covarrubias-Pazaran G. 2016. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6):1-15.
Covarrubias-Pazaran G. 2018. Software update: Moving the R package sommer to multivariate mixed models for genome-assisted prediction. doi: https://doi.org/10.1101/354639
Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp.
Gilmour et al. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51(4):1440-1450.
Henderson C.R. 1975. Best Linear Unbiased Estimation and Prediction under a Selection Model. Biometrics vol. 31(2):423-447.
Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.
Lee, D.-J., Durban, M., and Eilers, P.H.C. (2013). Efficient two-dimensional smoothing with P-spline ANOVA mixed models and nested bases. Computational Statistics and Data Analysis, 61, 22 - 37.
Lee et al. 2015. MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Cold Spring Harbor. doi: http://dx.doi.org/10.1101/027201.
Maier et al. 2015. Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder. Am J Hum Genet; 96(2):283-294.
Rodriguez-Alvarez, Maria Xose, et al. Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spatial Statistics 23 (2018): 52-71.
Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco.
Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208.
Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.