I was fitting mixed-effects models in a longitudinal dataset and noticed that specifying the same structure led to different outputs depending on if I used 'mixed' (from package afex) or lme4. I have provided reproducible code below (R version 4.2.2)
library(lme4)
library(afex)
library(tidyverse)
library(effectsize)
# Pull in toy data and convert vars to factors
dset_test1 <- mtcars %>%
mutate(car_id = row.names(mtcars)) %>%
mutate(vs = as.factor(vs)) %>%
mutate(am = as.factor(am))
# Create simulated longitudinal dataset
increase_longitudinal <- mtcars %>%
dplyr::select(am, wt) %>%
mutate(car_id = row.names(mtcars)) %>%
mutate(wt2 = (wt)*2)
# Merge and make long form
dset_test = merge(dset_test1, increase_longitudinal, on='car_id') %>%
pivot_longer(cols = c(wt, wt2), names_to = 'year', values_to = 'weight') %>%
mutate(year = as.factor(year))
# Set reference levels
dset_test$year <- relevel(dset_test$year, "wt")
# fit model with mixed
mix1 <- mixed(mpg~year * weight + hp +(1 | car_id), data=dset_test)
summary(mix1)
parameters::standardize_parameters(mix1)
# fit model with lme4
lmer1 <- lmer(mpg~year * weight + hp +(1 | car_id), data=dset_test)
summary(lmer1)
parameters::standardize_parameters(lmer1)
This produces the following standardized output:
# model with mixed:
Parameter | Std. Coef. | 95% CI
------------------------------------------
(Intercept) | -0.29 | [-0.40, -0.18]
year1 | -0.81 | [-0.98, -0.63]
weight | -1.12 | [-1.36, -0.88]
hp | -0.39 | [-0.53, -0.25]
year1:weight | -0.37 | [-0.45, -0.29]
# Model with lme4:
Parameter | Std. Coef. | 95% CI
--------------------------------------------
(Intercept) | -1.05 | [-1.29, -0.81]
yearwt2 | 1.54 | [ 1.21, 1.86]
weight | -1.42 | [-1.72, -1.12]
hp | -0.40 | [-0.54, -0.26]
yearwt2:weight | 0.71 | [ 0.56, 0.86]
As you can see, the model coefficients are quite different, particularly for the interaction term, where the sign of the coefficient even is different. Thanks for any insight into what might be going on here!
afex::mixed() automatically uses sum-to-zero contrasts by default, since these often make it easier to interpret main effects in models with interactions (but R's default is to use treatment contrasts).
If you set options(contrasts = c("contr.sum", "contr.poly")) before running the lmer() fit, you'll get identical coefficients. (There are other ways to set contrasts but this is probably the simplest.)
I hope by the way that this is not representative of your real data -- I got lots of warnings when I ran this code, and lots of red flags in the output (e.g. the residual standard error was very close to zero, as was the raw yearwt2 estimate in the original lmer fit). If you can it's generally best to generate more sensible test data so that people (like me) don't get distracted/think the problem might be some unidentifiability issue ...