I have a question on statistical mixed linear model analysis with nlme in R. There are a couple of posts that are related to my question, I think (see here and here), but I'm not sure they give the solution to my problem.
I would like to fit a linear mixed model to data from an experiment with crossed random effects for subjects and stimuli. Here are simplified data simulated based on my data and with great help from this paper and the accompanying R code.
library(tidyverse)
library(lme4)
library(nlme)
set.seed(1234)
# parameters for group
beta_0 <- 55.145 # intercept
omega_0 <- 3.644 # by-item random intercept sd
tau_0 <- 7.671 # by-subject random intercept sd
sigma <- 26.683 # residual (error) sd
# number of subjects and items
n_subj <- 150
n_items <- 36
items <- data.frame(
item_id = seq_len(n_items),
O_0i = rnorm(n = n_items, mean = 0, sd = omega_0)
)
subjects <- data.frame(
subj_id = seq_len(n_subj),
T_0i = rnorm(n = n_subj, mean = 0, sd = tau_0)
)
# add subject IDs
subjects1$subj_id <- seq_len(n_subj)
trials <- crossing(subjects, items) %>%
mutate(e_si = rnorm(nrow(.), mean = 0, sd = sigma))
dat_sim <- trials %>%
mutate(y = beta_0 + T_0i + O_0i + e_si) %>%
select(subj_id, item_id, y)
I would typically use lme4::lmer() to fit a model to this data but I have a between subjects factor in my experiment and the data turned out to be heteroscedastic. As far as I know, I need to use nlme::lme() to fit different residual variances for each between-subjects group, so lme4::lmer() is not an option. Now, I would like to fit a model in nlme::lme() that matches my original model in lme4::lmer().
For demonstration, I fit intercept-only models. (In reality I would have fixed effects, too. But this is not relevant for the question I have.)
My lmer() model would look like this:
m.lmer <- lmer(y ~ 1 + (1|subj_id) + (1|item_id),
data = dat_sim,
REML = T,
control = lmerControl(optimizer="nloptwrap"))
With the help of this post, I pieced this model together for lme():
dat_sim$dummy <- factor(1)
m.lme <- lme(y ~ 1,
random = list(dummy =
pdBlocked(list(pdIdent(~subj_id-1),
pdIdent(~item_id-1)))),
data = dat_sim,
method = "REML")
And here is my problem: The two models give different estimates for fixed and especially random effects:
summary(m.lmer)$coefficients
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 54.08182 0.8576796 102.7496 63.05597 5.819007e-84
VarCorr(m.lmer)
# Groups Name Std.Dev.
# subj_id (Intercept) 7.6215
# item_id (Intercept) 2.8036
# Residual 26.4984
summary(m.lme)
# Fixed effects: y ~ 1
# Value Std.Error DF t-value p-value
# (Intercept) 53.66397 0.9505988 5399 56.45281 0
VarCorr(m.lme)
# dummy = pdIdent(subj_id - 1), pdIdent(item_id - 1)
# Variance StdDev
# subj_id 6.408743e-04 0.02531550
# item_id 6.743363e-03 0.08211798
# Residual 7.655888e+02 27.66927555
But if I remove one of the random effects, my results are identical across lmer() and lme():
m.lmer.s <- update(m.lmer, ~ 1 + (1|subj_id))
VarCorr(m.lmer.s)
# Groups Name Std.Dev.
# subj_id (Intercept) 7.6072
# Residual 26.6463
m.lme.s <- update(m.lme, random = ~1|subj_id)
VarCorr(m.lme.s)
# subj_id = pdLogChol(1)
# Variance StdDev
# (Intercept) 57.86926 7.607185
# Residual 710.02435 26.646282
m.lmer.i <- update(m.lmer, ~ 1 + (1|item_id))
VarCorr(m.lmer.i)
# Groups Name Std.Dev.
# item_id (Intercept) 2.7336
# Residual 27.5727
m.lme.i <- update(m.lme, random = ~1|item_id)
VarCorr(m.lme.i)
# item_id = pdLogChol(1)
# Variance StdDev
# (Intercept) 7.472733 2.73363
# Residual 760.252039 27.57267
Why do the models with both random effects give so such different results and why does this seem to depend on whether only one or both random intercepts are fitted? Is my model specification for lme() wrong? Or is there an issue with a local maximum when fitting two random effects with lmer() (see this post)? Could anyone explain to me what is happening here?
