I have observed in the literature that the sommer package brings possibilities associated with the inference of the parameters of a mixed model, in comparison with classic packages: nlme or lme4.
Due of this, I'm trying to adjust some models via the mmer() function, which were previously adjusted via the lme() function.
Considering the nlme package, the following models have been adjusted:
data(Orthodont)
library(nlme)
m1 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ 1|Subject,
method = "REML")
# Different variances - separate variance per factor level
m2 <- lme(distance ~ Sex * age, data = Orthodont, random = list(Subject = pdDiag(~ age)),
method = "REML")
m3 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ age | Subject,
method = "REML")
In the sommer package, which can be accessed through the link: https://cran.r-project.org/web/packages/sommer/sommer.pdf, the following information is said at the end of page 9, beginning of page 10 and of page 47 and
In this link some additional information can also be found: https://rpubs.com/samuelkn/CovarianceStructuresInR:
#####===========================================####
#####Univariateheterogeneousvariancemodels####
#####===========================================####
### Compound simmetry(CS)+ Diagonal(DIAG) model
ans3<-mmer(Yield~Env, random=~Name+vsr(dsr(Env),Name),
rcov=~vsr(dsr(Env),units), data=DT)
summary(ans3)
ans4<-mmec(Yield~Env, random=~Name+vsc(dsc(Env),isc(Name)),
rcov=~vsc(dsc(Env),isc(units)), data=DT)
summary(ans4)
On page 38 there is information about the dsc() function and page 92 isc() function for the diagonal covariance structure. I tried to adjust the models and they looked like this:
data(Orthodont)
library(sommer)
m11<- mmer(distance ~ Sex * age, random = ~ vsr(Subject),
data = Orthodont, method = "REML")
m22 <- mmer(distance ~ Sex * age, random = ~vsr(dsc(age), isc(Subject)),
rcov =~vsr(dsc(age), units),
data = Orthodont, method = "REML")
m33 <- mmer(distance ~ Sex * age, random = ~ Subject + vsr(age, Subject),
data = Orthodont, method = "REML")
One fact is that:
- For the m1 and m11 I obtained satisfactory results:
> summary(m1)
Linear mixed-effects model fit by REML
Data: Orthodont
AIC BIC logLik
445.7572 461.6236 -216.8786
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 1.816214 1.386382
Fixed effects: distance ~ Sex * age
Value Std.Error DF t-value p-value
(Intercept) 16.340625 0.9813122 79 16.651810 0.0000
SexFemale 1.032102 1.5374208 25 0.671321 0.5082
age 0.784375 0.0775011 79 10.120823 0.0000
SexFemale:age -0.304830 0.1214209 79 -2.510520 0.0141
Correlation:
(Intr) SexFml age
SexFemale -0.638
age -0.869 0.555
SexFemale:age 0.555 -0.869 -0.638
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.59804400 -0.45461690 0.01578365 0.50244658 3.68620792
Number of Observations: 108
Number of Groups: 27
> summary(m11)
=============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.2 **********************
=============================================================
logLik AIC BIC Method Converge
Value -0.008409636 8.016819 18.74534 REML TRUE
=============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:Subject.distance-distance 3.198 2.76014 1.159 Positive
units.distance-distance 1.914 0.09586 19.966 Positive
=============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 distance (Intercept) 16.3406 0.97646 16.7346
2 distance SexFemale 1.0321 1.52982 0.6747
3 distance age 0.7844 0.07734 10.1422
4 distance SexFemale:age -0.3048 0.12116 -2.5158
=============================================================
Groups and observations:
distance
u:Subject 27
=============================================================
Use the '$' sign to access results and parameters
- For the m3 and m33 I obtained satisfactory results
for the fixed part, however for the random part there are still differences, based on the information found in this link of the page 2: https://cran.r-project.org/package=sommer/vignettes/v5.sommer.vs.lme4.pdf :
> summary(m3)
Linear mixed-effects model fit by REML
Data: Orthodont
AIC BIC logLik
448.5817 469.7368 -216.2908
Random effects:
Formula: ~age | Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 2.4055009 (Intr)
age 0.1803455 -0.668
Residual 1.3100396
Fixed effects: distance ~ Sex * age
Value Std.Error DF t-value p-value
(Intercept) 16.340625 1.0185320 79 16.043311 0.0000
SexFemale 1.032102 1.5957329 25 0.646789 0.5237
age 0.784375 0.0859995 79 9.120691 0.0000
SexFemale:age -0.304830 0.1347353 79 -2.262432 0.0264
Correlation:
(Intr) SexFml age
SexFemale -0.638
age -0.880 0.562
SexFemale:age 0.562 -0.880 -0.638
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.168077732 -0.385939009 0.007104087 0.445154545 3.849463576
Number of Observations: 108
Number of Groups: 27
> summary(m33)
===============================================================
Multivariate Linear Mixed Model fit by REML
************************ sommer 4.2 ************************
===============================================================
logLik AIC BIC Method Converge
Value -0.1589647 8.317929 19.04645 REML TRUE
===============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Subject.distance-distance 2.434714 1.394632 1.746 Positive
age:Subject.distance-distance 0.009851 0.009058 1.088 Positive
units.distance-distance 2.007921 0.118293 16.974 Positive
===============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 distance (Intercept) 16.3406 0.97097 16.8291
2 distance SexFemale 1.0321 1.52122 0.6785
3 distance age 0.7844 0.08301 9.4493
4 distance SexFemale:age -0.3048 0.13005 -2.3439
===============================================================
Groups and observations:
distance
Subject 27
age:Subject 27
===============================================================
Use the '$' sign to access results and parameters
Checking random effects m3 and m33:
> VarCorr(m3)
Subject = pdLogChol(age)
Variance StdDev Corr
(Intercept) 5.78643477 2.4055009 (Intr)
age 0.03252449 0.1803455 -0.668
Residual 1.71620366 1.3100396
> sqrt(summary(m33)$varcomp[,1])
[1] 1.55463405 0.08802067 1.36549133
However, for model 2, I have been facing problem in structuring this model, compared to the model adjusted via lme() function. Note that the models do not exhibit the same coefficients, especially the random effect:
> summary(m2)
Linear mixed-effects model fit by REML
Data: Orthodont
AIC BIC logLik
447.1509 465.6617 -216.5755
Random effects:
Formula: ~age | Subject
Structure: Diagonal
(Intercept) age Residual
StdDev: 1.554607 0.08801657 1.365502
Fixed effects: distance ~ Sex * age
Value Std.Error DF t-value p-value
(Intercept) 16.340625 0.9408690 79 17.367587 0.0000
SexFemale 1.032102 1.4740585 25 0.700177 0.4903
age 0.784375 0.0794421 79 9.873548 0.0000
SexFemale:age -0.304830 0.1244618 79 -2.449182 0.0165
Correlation:
(Intr) SexFml age
SexFemale -0.638
age -0.858 0.547
SexFemale:age 0.547 -0.858 -0.638
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.623218206 -0.458236586 0.004989037 0.517505851 3.716965299
Number of Observations: 108
Number of Groups: 27
> summary(m22)
===============================================================
Multivariate Linear Mixed Model fit by REML
************************ sommer 4.2 ************************
===============================================================
logLik AIC BIC Method Converge
Value -4.161598 16.3232 27.05172 REML TRUE
===============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
age:Subject.distance-distance 0.02684 0.01981 1.355 Positive
age:units.distance-distance 0.19972 0.01017 19.633 Positive
===============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 distance (Intercept) 16.5180 0.87982 18.7744
2 distance SexFemale 0.8419 1.37841 0.6108
3 distance age 0.7682 0.09143 8.4030
4 distance SexFemale:age -0.2875 0.14324 -2.0075
===============================================================
Groups and observations:
distance
age:Subject 27
===============================================================
Use the '$' sign to access results and parameters
>
The sommer package version I am using is:
> packageVersion("sommer")
[1] ‘4.3.1’
Conclusion: Unfortunately I don't know what could be happening.
I think your model 2 is fitted like this
m22 <- mmer(distance ~ Sex * age,
random = ~ vsr(Subject) + vsr(dsr(age), Subject) ,
data = Orthodont)
summary(m22)
sqrt(summary(m22)$varcomp[,1])
That returns the same variance components and estimates than the model 2 using lme()
# > summary(m22)
# ===============================================================
# Multivariate Linear Mixed Model fit by REML
# ************************ sommer 4.2 ************************
# ===============================================================
# logLik AIC BIC Method Converge
# Value 0.2963668 7.407266 18.13579 NR TRUE
# ===============================================================
# Variance-Covariance components:
# VarComp VarCompSE Zratio Constraint
# u:Subject.distance-distance 2.416887 1.56526 1.5441 Positive
# age:Subject.distance-distance 0.007748 0.01122 0.6902 Positive
# units.distance-distance 1.864567 0.30408 6.1319 Positive
# ===============================================================
# Fixed effects:
# Trait Effect Estimate Std.Error t.value
# 1 distance (Intercept) 16.3406 0.94087 17.3676
# 2 distance SexFemale 1.0321 1.47405 0.7002
# 3 distance age 0.7844 0.07944 9.8736
# 4 distance SexFemale:age -0.3048 0.12446 -2.4492
# ===============================================================
# Groups and observations:
# distance
# u:Subject 27
# age:Subject 27
# ===============================================================
# Use the '$' sign to access results and parameters
# > sqrt(summary(m22)$varcomp[,1])
# [1] 1.55463405 0.08802067 1.36549133
That is using the mmer() function which uses direct inversion. If you want to use the mmec() function which uses Henderson's approach the model would be:
m22 <- mmec(distance ~ Sex * age,
random = ~ vsc(isc(Subject)) + vsc(dsc(age), isc(Subject)) ,
# emWeight = rep(1,30),
data = Orthodont)
summary(m22)$varcomp
sqrt(summary(m22)$varcomp[,1])