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 (formulation=“direct”) 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,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError Zratio
## 1 vsm(ism(Subject)) sigma2 1378.2022 505.7744 2.724935
## 2 vsm(ism(units)) sigma2 960.3449 107.0296 8.972701This 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,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError Zratio
## 1 vsm(ism(Subject)) sigma2 627.59783 283.65666 2.212526
## 2 vsm(ism(Days), ism(Subject)) sigma2 35.85785 14.52919 2.467987
## 3 vsm(ism(units)) sigma2 653.58081 76.65691 8.526052Notice 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,
# control = mmes_control(formulation="henderson"),
random= ~ vsm(usm(cbind(1,Days)), ism(Subject)),
data=DT,
control = mmes_control(convergence = list(maxIters = 200),
execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError
## 1 vsm(usm(cbind(1, Days)), ism(Subject)) variance[] 612.075419 288.77522
## 2 vsm(usm(cbind(1, Days)), ism(Subject)) covariance[Days,] 9.600784 46.67787
## 3 vsm(usm(cbind(1, Days)), ism(Subject)) variance[Days] 35.071429 14.78194
## 4 vsm(ism(units)) sigma2 654.941027 77.18554
## Zratio
## 1 2.1195566
## 2 0.2056817
## 3 2.3725859
## 4 8.4852814## [,1] [,2]
## [1,] 1.00000000 0.06552809
## [2,] 0.06552809 1.00000000Notice 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,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError Zratio
## 1 vsm(ism(Days), ism(Subject)) sigma2 52.70784 19.09481 2.760322
## 2 vsm(ism(units)) sigma2 842.02116 93.84754 8.972224One 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,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError Zratio
## 1 vsm(dsm(Daysf), ism(Subject)) variance[0] 656.6223 114.7957 5.719918
## 2 vsm(dsm(Daysf), ism(Subject)) variance[1] 714.5969 446.0912 1.601908
## 3 vsm(dsm(Daysf), ism(Subject)) variance[2] 518.3159 384.8090 1.346943
## 4 vsm(dsm(Daysf), ism(Subject)) variance[3] 1073.6579 557.7301 1.925049
## 5 vsm(dsm(Daysf), ism(Subject)) variance[4] 1372.7595 652.1531 2.104965
## 6 vsm(dsm(Daysf), ism(Subject)) variance[5] 2216.7852 922.5747 2.402825
## 7 vsm(dsm(Daysf), ism(Subject)) variance[6] 3428.2196 1317.9499 2.601176
## 8 vsm(dsm(Daysf), ism(Subject)) variance[7] 2056.9049 873.4375 2.354954
## 9 vsm(dsm(Daysf), ism(Subject)) variance[8] 3114.8998 1219.6650 2.553898
## 10 vsm(dsm(Daysf), ism(Subject)) variance[9] 3990.1124 1511.7008 2.639485
## 11 vsm(ism(units)) sigma2 363.4387 256.2290 1.418413## unstructured model
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(usm(Daysf), ism(Subject)),
data=DT,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError
## 1 vsm(usm(Daysf), ism(Subject)) variance[0] 1264.36579 451.426562
## 2 vsm(usm(Daysf), ism(Subject)) covariance[1,0] 1021.07525 410.055185
## 3 vsm(usm(Daysf), ism(Subject)) covariance[2,0] 538.67488 301.274843
## 4 vsm(usm(Daysf), ism(Subject)) covariance[3,0] 796.19323 416.166401
## 5 vsm(usm(Daysf), ism(Subject)) covariance[4,0] 755.53781 428.110812
## 6 vsm(usm(Daysf), ism(Subject)) covariance[5,0] 929.53916 541.129616
## 7 vsm(usm(Daysf), ism(Subject)) covariance[6,0] 662.81699 585.404176
## 8 vsm(usm(Daysf), ism(Subject)) covariance[7,0] 928.73204 513.795947
## 9 vsm(usm(Daysf), ism(Subject)) covariance[8,0] 915.96704 606.708346
## 10 vsm(usm(Daysf), ism(Subject)) covariance[9,0] 1435.77806 721.190980
## 11 vsm(usm(Daysf), ism(Subject)) variance[1] 1279.10126 454.752887
## 12 vsm(usm(Daysf), ism(Subject)) covariance[2,1] 826.95699 335.747397
## 13 vsm(usm(Daysf), ism(Subject)) covariance[3,1] 1135.01993 459.290989
## 14 vsm(usm(Daysf), ism(Subject)) covariance[4,1] 1037.13917 460.539007
## 15 vsm(usm(Daysf), ism(Subject)) covariance[5,1] 1176.14002 569.804981
## 16 vsm(usm(Daysf), ism(Subject)) covariance[6,1] 847.48780 600.354698
## 17 vsm(usm(Daysf), ism(Subject)) covariance[7,1] 923.74428 512.739680
## 18 vsm(usm(Daysf), ism(Subject)) covariance[8,1] 1040.45401 619.232958
## 19 vsm(usm(Daysf), ism(Subject)) covariance[9,1] 1509.90053 731.924824
## 20 vsm(usm(Daysf), ism(Subject)) variance[2] 862.15773 304.509328
## 21 vsm(usm(Daysf), ism(Subject)) covariance[3,2] 1055.31726 393.402487
## 22 vsm(usm(Daysf), ism(Subject)) covariance[4,2] 909.53160 382.460250
## 23 vsm(usm(Daysf), ism(Subject)) covariance[5,2] 856.73615 454.157810
## 24 vsm(usm(Daysf), ism(Subject)) covariance[6,2] 914.07015 509.755770
## 25 vsm(usm(Daysf), ism(Subject)) covariance[7,2] 922.60757 432.426478
## 26 vsm(usm(Daysf), ism(Subject)) covariance[8,2] 828.63962 500.122114
## 27 vsm(usm(Daysf), ism(Subject)) covariance[9,2] 964.49806 566.413839
## 28 vsm(usm(Daysf), ism(Subject)) variance[3] 1620.50093 571.455378
## 29 vsm(usm(Daysf), ism(Subject)) covariance[4,3] 1587.82979 575.755371
## 30 vsm(usm(Daysf), ism(Subject)) covariance[5,3] 1668.86239 684.508749
## 31 vsm(usm(Daysf), ism(Subject)) covariance[6,3] 1781.99755 762.453839
## 32 vsm(usm(Daysf), ism(Subject)) covariance[7,3] 1279.10305 599.598873
## 33 vsm(usm(Daysf), ism(Subject)) covariance[8,3] 1602.00789 741.343417
## 34 vsm(usm(Daysf), ism(Subject)) covariance[9,3] 1736.93532 826.406717
## 35 vsm(usm(Daysf), ism(Subject)) variance[4] 1817.69214 630.906634
## 36 vsm(usm(Daysf), ism(Subject)) covariance[5,4] 1998.93237 752.208367
## 37 vsm(usm(Daysf), ism(Subject)) covariance[6,4] 2073.99885 827.330131
## 38 vsm(usm(Daysf), ism(Subject)) covariance[7,4] 1546.78126 651.454084
## 39 vsm(usm(Daysf), ism(Subject)) covariance[8,4] 2024.50880 819.223350
## 40 vsm(usm(Daysf), ism(Subject)) covariance[9,4] 2193.04658 910.546602
## 41 vsm(usm(Daysf), ism(Subject)) variance[5] 2925.09152 1024.793000
## 42 vsm(usm(Daysf), ism(Subject)) covariance[6,5] 2597.76647 1050.869248
## 43 vsm(usm(Daysf), ism(Subject)) covariance[7,5] 1936.35997 833.204292
## 44 vsm(usm(Daysf), ism(Subject)) covariance[8,5] 3051.03781 1120.627007
## 45 vsm(usm(Daysf), ism(Subject)) covariance[9,5] 3228.67085 1230.348379
## 46 vsm(usm(Daysf), ism(Subject)) variance[6] 3981.62872 1360.761703
## 47 vsm(usm(Daysf), ism(Subject)) covariance[7,6] 2301.72051 958.764170
## 48 vsm(usm(Daysf), ism(Subject)) covariance[8,6] 2921.42474 1191.860370
## 49 vsm(usm(Daysf), ism(Subject)) covariance[9,6] 2203.46316 1209.083294
## 50 vsm(usm(Daysf), ism(Subject)) variance[7] 2528.76668 880.129570
## 51 vsm(usm(Daysf), ism(Subject)) covariance[8,7] 2427.56796 976.466848
## 52 vsm(usm(Daysf), ism(Subject)) covariance[9,7] 2392.74649 1054.491779
## 53 vsm(usm(Daysf), ism(Subject)) variance[8] 3802.35473 1331.886895
## 54 vsm(usm(Daysf), ism(Subject)) covariance[9,8] 3837.72642 1426.439075
## 55 vsm(usm(Daysf), ism(Subject)) variance[9] 4798.60167 1685.875458
## 56 vsm(ism(units)) sigma2 20.09991 6.638579
## Zratio
## 1 2.800823
## 2 2.490092
## 3 1.787985
## 4 1.913161
## 5 1.764818
## 6 1.717775
## 7 1.132238
## 8 1.807589
## 9 1.509732
## 10 1.990843
## 11 2.812739
## 12 2.463033
## 13 2.471244
## 14 2.252012
## 15 2.064110
## 16 1.411645
## 17 1.801585
## 18 1.680230
## 19 2.062917
## 20 2.831302
## 21 2.682538
## 22 2.378108
## 23 1.886428
## 24 1.793153
## 25 2.133559
## 26 1.656875
## 27 1.702815
## 28 2.835744
## 29 2.757820
## 30 2.438044
## 31 2.337187
## 32 2.133265
## 33 2.160952
## 34 2.101792
## 35 2.881079
## 36 2.657418
## 37 2.506858
## 38 2.374352
## 39 2.471254
## 40 2.408495
## 41 2.854324
## 42 2.472017
## 43 2.323992
## 44 2.722617
## 45 2.624192
## 46 2.926029
## 47 2.400716
## 48 2.451147
## 49 1.822425
## 50 2.873175
## 51 2.486073
## 52 2.269099
## 53 2.854863
## 54 2.690424
## 55 2.846356
## 56 3.027743## random regression (legendre polynomials)
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(dsm(leg(Days,1)), ism(Subject)),
data=DT,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate StdError
## 1 vsm(dsm(leg(Days, 1)), ism(Subject)) variance[leg0] 2817.4417 1011.41776
## 2 vsm(dsm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.4674 199.55780
## 3 vsm(ism(units)) sigma2 654.9410 77.18554
## Zratio
## 1 2.785636
## 2 2.372583
## 3 8.485282## unstructured random regression (legendre)
fm2 <- mmes(Reaction ~ Days,
random= ~ vsm(usm(leg(Days,1)), ism(Subject)),
data=DT,
control = mmes_control(execution = list(inverseTol = 1e-6),
diagnostics = list(verbose = FALSE)))
summary(fm2)$varcomp## term parameter estimate
## 1 vsm(usm(leg(Days, 1)), ism(Subject)) variance[leg0] 2817.4277
## 2 vsm(usm(leg(Days, 1)), ism(Subject)) covariance[leg1,leg0] 869.9652
## 3 vsm(usm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.4643
## 4 vsm(ism(units)) sigma2 654.9410
## StdError Zratio
## 1 1011.41057 2.785642
## 2 381.08010 2.282893
## 3 199.55625 2.372586
## 4 77.18554 8.485281Covarrubias-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.