Translating lme4 models to sommer
Giovanny Covarrubias-Pazaran
2026-10-06
Source:vignettes/sommer.vs.lme4.Rmd
sommer.vs.lme4.RmdThe 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.
- Random slopes with same intercept
- Random slopes and random intercepts (without correlation)
- Random slopes and random intercepts (with correlation)
- Random slopes with a different intercept
- Other models not available in lme4
1) Random slopes
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.
## Loading required package: Matrix## Loading required package: MASS## Loading required package: crayon## Loading required package: enhancer
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.9628 505.6592 2.725082
## 2 vsm(ism(units)) sigma2 960.4982 107.0553 8.9719852) Random slopes and random intercepts (without correlation)
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.55262 283.63263 2.212554
## 2 vsm(ism(Days), ism(Subject)) sigma2 35.85476 14.52757 2.468050
## 3 vsm(ism(units)) sigma2 653.58363 76.65735 8.526040Notice that Days is a numerical (not factor) variable.
3) Random slopes and random intercepts (with correlation)
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[] 610.776208 288.11618
## 2 vsm(usm(cbind(1, Days)), ism(Subject)) covariance[Days,] 9.640222 46.64122
## 3 vsm(usm(cbind(1, Days)), ism(Subject)) variance[Days] 35.093150 14.79332
## 4 vsm(ism(units)) sigma2 654.975406 77.19162
## Zratio
## 1 2.1198956
## 2 0.2066889
## 3 2.3722301
## 4 8.4850587
cov2cor(fm2$theta[[1]])## [,1] [,2]
## [1,] 1.00000000 0.06584682
## [2,] 0.06584682 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.
4) Random slopes with a different intercept
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.70874 19.09541 2.760284
## 2 vsm(ism(units)) sigma2 842.09182 93.85935 8.9718484) Other models available in sommer but not in lme4
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.8655 74.68294 6.7333385
## 2 vsm(dsm(Daysf), ism(Subject)) variance[1] 558.3574 503.32816 1.1093307
## 3 vsm(dsm(Daysf), ism(Subject)) variance[2] 361.3899 439.53385 0.8222118
## 4 vsm(dsm(Daysf), ism(Subject)) variance[3] 917.5182 607.46547 1.5104038
## 5 vsm(dsm(Daysf), ism(Subject)) variance[4] 1216.6107 695.59696 1.7490167
## 6 vsm(dsm(Daysf), ism(Subject)) variance[5] 2060.8449 953.76616 2.1607444
## 7 vsm(dsm(Daysf), ism(Subject)) variance[6] 3272.2788 1339.70886 2.4425298
## 8 vsm(dsm(Daysf), ism(Subject)) variance[7] 1901.0460 906.31736 2.0975500
## 9 vsm(dsm(Daysf), ism(Subject)) variance[8] 2958.9116 1243.29915 2.3798871
## 10 vsm(dsm(Daysf), ism(Subject)) variance[9] 3834.0375 1530.74664 2.5046846
## 11 vsm(ism(units)) sigma2 518.0607 345.48785 1.4995048
## 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
## 1 vsm(usm(Daysf), ism(Subject)) variance[0] 813.26915 248.268843
## 2 vsm(usm(Daysf), ism(Subject)) covariance[1,0] 696.13240 264.507108
## 3 vsm(usm(Daysf), ism(Subject)) covariance[2,0] 393.66573 230.272132
## 4 vsm(usm(Daysf), ism(Subject)) covariance[3,0] 599.99763 338.094848
## 5 vsm(usm(Daysf), ism(Subject)) covariance[4,0] 602.56858 368.518591
## 6 vsm(usm(Daysf), ism(Subject)) covariance[5,0] 735.88291 463.761458
## 7 vsm(usm(Daysf), ism(Subject)) covariance[6,0] 576.93467 516.791239
## 8 vsm(usm(Daysf), ism(Subject)) covariance[7,0] 726.51241 424.276809
## 9 vsm(usm(Daysf), ism(Subject)) covariance[8,0] 759.00751 534.840892
## 10 vsm(usm(Daysf), ism(Subject)) covariance[9,0] 1101.81778 580.698304
## 11 vsm(usm(Daysf), ism(Subject)) variance[1] 1094.90618 379.120151
## 12 vsm(usm(Daysf), ism(Subject)) covariance[2,1] 811.05242 328.864405
## 13 vsm(usm(Daysf), ism(Subject)) covariance[3,1] 1143.25846 482.394412
## 14 vsm(usm(Daysf), ism(Subject)) covariance[4,1] 1106.21307 515.371025
## 15 vsm(usm(Daysf), ism(Subject)) covariance[5,1] 1251.17647 637.517153
## 16 vsm(usm(Daysf), ism(Subject)) covariance[6,1] 1029.32522 690.170378
## 17 vsm(usm(Daysf), ism(Subject)) covariance[7,1] 947.33195 539.809559
## 18 vsm(usm(Daysf), ism(Subject)) covariance[8,1] 1176.93714 713.737050
## 19 vsm(usm(Daysf), ism(Subject)) covariance[9,1] 1495.44833 761.461145
## 20 vsm(usm(Daysf), ism(Subject)) variance[2] 952.89796 345.359515
## 21 vsm(usm(Daysf), ism(Subject)) covariance[3,2] 1230.35117 479.804405
## 22 vsm(usm(Daysf), ism(Subject)) covariance[4,2] 1153.04191 502.171855
## 23 vsm(usm(Daysf), ism(Subject)) covariance[5,2] 1142.24632 606.089905
## 24 vsm(usm(Daysf), ism(Subject)) covariance[6,2] 1290.09025 687.592858
## 25 vsm(usm(Daysf), ism(Subject)) covariance[7,2] 1146.66492 536.105578
## 26 vsm(usm(Daysf), ism(Subject)) covariance[8,2] 1186.77270 687.932486
## 27 vsm(usm(Daysf), ism(Subject)) covariance[9,2] 1220.32738 709.947241
## 28 vsm(usm(Daysf), ism(Subject)) variance[3] 1936.51707 762.088845
## 29 vsm(usm(Daysf), ism(Subject)) covariance[4,3] 2016.88391 826.492564
## 30 vsm(usm(Daysf), ism(Subject)) covariance[5,3] 2191.55305 996.440052
## 31 vsm(usm(Daysf), ism(Subject)) covariance[6,3] 2439.57295 1132.667562
## 32 vsm(usm(Daysf), ism(Subject)) covariance[7,3] 1680.73820 828.627064
## 33 vsm(usm(Daysf), ism(Subject)) covariance[8,3] 2251.60184 1120.215250
## 34 vsm(usm(Daysf), ism(Subject)) covariance[9,3] 2239.32204 1132.197416
## 35 vsm(usm(Daysf), ism(Subject)) variance[4] 2361.71439 950.755365
## 36 vsm(usm(Daysf), ism(Subject)) covariance[5,4] 2678.23675 1146.488257
## 37 vsm(usm(Daysf), ism(Subject)) covariance[6,4] 2895.61635 1296.195149
## 38 vsm(usm(Daysf), ism(Subject)) covariance[7,4] 2081.48667 964.914774
## 39 vsm(usm(Daysf), ism(Subject)) covariance[8,4] 2852.34245 1297.034922
## 40 vsm(usm(Daysf), ism(Subject)) covariance[9,4] 2867.10709 1306.796839
## 41 vsm(usm(Daysf), ism(Subject)) variance[5] 3769.00723 1507.850194
## 42 vsm(usm(Daysf), ism(Subject)) covariance[6,5] 3617.02502 1629.131158
## 43 vsm(usm(Daysf), ism(Subject)) covariance[7,5] 2605.11633 1224.376240
## 44 vsm(usm(Daysf), ism(Subject)) covariance[8,5] 4096.54340 1710.063033
## 45 vsm(usm(Daysf), ism(Subject)) covariance[9,5] 4092.86377 1717.156698
## 46 vsm(usm(Daysf), ism(Subject)) variance[6] 5178.44322 2046.223653
## 47 vsm(usm(Daysf), ism(Subject)) covariance[7,6] 3131.57428 1429.049709
## 48 vsm(usm(Daysf), ism(Subject)) covariance[8,6] 4148.86148 1892.084039
## 49 vsm(usm(Daysf), ism(Subject)) covariance[9,6] 3260.99049 1785.576480
## 50 vsm(usm(Daysf), ism(Subject)) variance[7] 3058.80838 1193.824166
## 51 vsm(usm(Daysf), ism(Subject)) covariance[8,7] 3267.11660 1460.611643
## 52 vsm(usm(Daysf), ism(Subject)) covariance[9,7] 3051.94729 1441.676772
## 53 vsm(usm(Daysf), ism(Subject)) variance[8] 5078.98313 2047.284253
## 54 vsm(usm(Daysf), ism(Subject)) covariance[9,8] 4930.60559 2027.462275
## 55 vsm(usm(Daysf), ism(Subject)) variance[9] 5637.13534 2154.429391
## 56 vsm(ism(units)) sigma2 27.84182 3.891448
## Zratio
## 1 3.275760
## 2 2.631810
## 3 1.709567
## 4 1.774643
## 5 1.635110
## 6 1.586770
## 7 1.116379
## 8 1.712355
## 9 1.419128
## 10 1.897401
## 11 2.888019
## 12 2.466221
## 13 2.369966
## 14 2.146440
## 15 1.962577
## 16 1.491407
## 17 1.754937
## 18 1.648979
## 19 1.963919
## 20 2.759148
## 21 2.564277
## 22 2.296110
## 23 1.884615
## 24 1.876241
## 25 2.138879
## 26 1.725130
## 27 1.718899
## 28 2.541065
## 29 2.440293
## 30 2.199383
## 31 2.153830
## 32 2.028341
## 33 2.009972
## 34 1.977855
## 35 2.484040
## 36 2.336035
## 37 2.233935
## 38 2.157172
## 39 2.199125
## 40 2.193996
## 41 2.499590
## 42 2.220217
## 43 2.127709
## 44 2.395551
## 45 2.383512
## 46 2.530732
## 47 2.191368
## 48 2.192747
## 49 1.826296
## 50 2.562193
## 51 2.236814
## 52 2.116943
## 53 2.480839
## 54 2.431910
## 55 2.616533
## 56 7.154617
## 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] 2813.9872 1009.64062
## 2 vsm(dsm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.7831 199.71947
## 3 vsm(ism(units)) sigma2 654.9355 77.18456
## Zratio
## 1 2.787118
## 2 2.372243
## 3 8.485317
## 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] 2814.1252
## 2 vsm(usm(leg(Days, 1)), ism(Subject)) covariance[leg1,leg0] 868.9816
## 3 vsm(usm(leg(Days, 1)), ism(Subject)) variance[leg1] 473.3727
## 4 vsm(ism(units)) sigma2 654.9173
## StdError Zratio
## 1 1009.70965 2.787064
## 2 380.62715 2.283026
## 3 199.50722 2.372710
## 4 77.18134 8.485435Literature
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.