I am using a dataset from the sommer package in R, where Env is a factor of 3 levels and Name is also a factor.
library(lme4)
library(sommer)
data(DT_example)
DT <- DT_example
mod_a <- lmer(Yield ~ Env + (Env | Name), data = DT)
print(VarCorr(mod_a), comp = "Variance")
Groups Name Variance Cov
Name (Intercept) 15.9918
EnvCA.2012 8.9229 -9.820
EnvCA.2013 10.9479 -9.626 3.829
Residual 4.3863
mod_b <- lmer(Yield ~ Env + (0 + Env | Name), data = DT)
print(VarCorr(mod_b), comp = "Variance")
Groups Name Variance Cov
Name EnvCA.2011 15.9930
EnvCA.2012 5.2743 6.172
EnvCA.2013 7.6897 6.366 0.375
Residual 4.3862
mod_c <- lmer(Yield ~ Env + (1 | Name/Env), data = DT)
print(VarCorr(mod_c), comp = "Variance")
Groups Name Variance
Env:Name (Intercept) 5.1732
Name (Intercept) 3.6819
Residual 4.3662
In mod_a, there is not variance parameter for the first level of Env. In mod_b, there is no random-intercept variance. If I understand correctly, mod_c is reasonable only when Env is nested in Name. There is only one variance parameter returned for Env.
In sommer, I can obtain variance parameter of random-intercept for Name and varance for each level of Env through:
mod_d <- mmer(Yield ~ Env, random = ~ Name + vsr(dsr(Env), Name),
rcov = ~ units, data = DT, verbose = FALSE)
summary(mod_d)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.3 **********************
============================================================
logLik AIC BIC Method Converge
Value -18.16164 42.32327 51.98434 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Name.Yield-Yield 2.965 1.5055 1.969 Positive
CA.2011:Name.Yield-Yield 10.424 4.4544 2.340 Positive
CA.2012:Name.Yield-Yield 2.658 1.8032 1.474 Positive
CA.2013:Name.Yield-Yield 5.702 2.5113 2.271 Positive
units.Yield-Yield 4.398 0.6517 6.748 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.511 0.8269 19.967
2 Yield EnvCA.2012 -5.809 0.8593 -6.760
3 Yield EnvCA.2013 -6.423 0.9358 -6.864
============================================================
Groups and observations:
Yield
Name 41
CA.2011:Name 41
CA.2012:Name 41
CA.2013:Name 41
============================================================
Use the '$' sign to access results and parameters
I want to ask is it possible to do the same in nlme or lme4?
To be honest I'm surprised that all the variance components in this model are identifiable, but they seem to be (I'll have to think about it).
This is slightly easier in glmmTMB (which has a built-in capability for fitting multivariate RE with diagonal covariance), but can also be hacked without too much trouble in lme4 using the framework described in ?modular:
library(lme4)
library(sommer)
library(glmmTMB)
library(broom.mixed)
data(DT_example)
fit1 <- glmmTMB(Yield ~ Env + (1|Name) + diag(0 + Env | Name), data = DT_example,
REML = TRUE)
## variances
tidy(fit1, effects = "ran_pars", scales = "vcov")
lmod <- lFormula(Yield ~ Env + (1|Name) + (0 + Env | Name), data = DT_example, REML = TRUE)
df <- do.call(mkLmerDevfun, lmod)
## wrapper function: take parameters for a diagonal cov matrix,
## pass `df()` the values it expects plus zeros for the covariances
df2 <- function(theta) {
d <- diag(theta[2:4])
th <- c(theta[1], d[lower.tri(d, diag = TRUE)])
df(th)
}
df2(c(1,1,1,1)) ## test deviance function
opt <-nloptwrap(c(1,1,1,1), fn = df2, lower = rep(0, 4), upper = rep(Inf, 4))
fit2 <- mkMerMod(environment(df), opt, lmod$reTrms, fr = lmod$fr)
tidy(fit2, effects = "ran_pars", scales = "vcov")