I want to obtain a monotonic effect for an ordinal predictor with the function lmer() from the package lme4. My reference is the estimate that can be obtained with mo() in the brms package (a more developed review of monotonic effects for ordinal predictors can be found here https://psyarxiv.com/9qkhj/).
Below a reprex with the desired estimate (leaving aside the different statistical approach behind both packages) in m1, and what happens by default (m2) when an ordered factor is used lmer()
library(brms)
library(lme4)
sleepstudy$Days <- factor(sleepstudy$Days, ordered = T)
m1 <- brm(Reaction ~ mo(Days) + (1 | Subject), data = sleepstudy, chains = 2)
#> Compiling Stan program...
summary(m1)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: Reaction ~ mo(Days) + (1 | Subject)
#> Data: sleepstudy (Number of observations: 180)
#> Samples: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup samples = 2000
#>
#> Group-Level Effects:
#> ~Subject (Number of levels: 18)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 39.73 7.85 27.73 58.76 1.01 538 751
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 257.11 10.87 235.14 277.72 1.00 468 774
#> moDays 10.12 1.01 8.16 12.09 1.00 1603 1315
#>
#> Simplex Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> moDays1[1] 0.07 0.06 0.00 0.20 1.00 1343 747
#> moDays1[2] 0.07 0.05 0.00 0.20 1.00 1275 524
#> moDays1[3] 0.13 0.07 0.01 0.29 1.00 1337 591
#> moDays1[4] 0.10 0.07 0.00 0.28 1.00 1600 850
#> moDays1[5] 0.16 0.09 0.01 0.34 1.00 1285 658
#> moDays1[6] 0.09 0.06 0.00 0.24 1.00 1543 840
#> moDays1[7] 0.09 0.07 0.00 0.25 1.00 1534 992
#> moDays1[8] 0.16 0.08 0.02 0.32 1.00 1897 906
#> moDays1[9] 0.13 0.08 0.01 0.31 1.00 1839 936
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 31.25 1.81 27.93 35.19 1.00 1726 1341
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
m2 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
summary(m2)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (1 | Subject)
#> Data: sleepstudy
#>
#> REML criterion at convergence: 1731.8
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.3473 -0.5293 0.0317 0.5328 4.2570
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> Subject (Intercept) 1375.5 37.09
#> Residual 987.6 31.43
#> Number of obs: 180, groups: Subject, 18
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 298.508 9.050 32.985
#> Days.L 95.074 7.407 12.835
#> Days.Q 7.744 7.407 1.045
#> Days.C -0.705 7.407 -0.095
#> Days^4 5.889 7.407 0.795
#> Days^5 1.754 7.407 0.237
#> Days^6 -6.036 7.407 -0.815
#> Days^7 -1.695 7.407 -0.229
#> Days^8 -4.161 7.407 -0.562
#> Days^9 6.435 7.407 0.869
#>
#> Correlation of Fixed Effects:
#> (Intr) Days.L Days.Q Days.C Days^4 Days^5 Days^6 Days^7 Days^8
#> Days.L 0.000
#> Days.Q 0.000 0.000
#> Days.C 0.000 0.000 0.000
#> Days^4 0.000 0.000 0.000 0.000
#> Days^5 0.000 0.000 0.000 0.000 0.000
#> Days^6 0.000 0.000 0.000 0.000 0.000 0.000
#> Days^7 0.000 0.000 0.000 0.000 0.000 0.000 0.000
#> Days^8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
#> Days^9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
How could an equivalent monotonic effect of moDays in m1 can be obtained with lme4?
It is not possible to do this with lme4::lmer(), but an alternative can be implemented with the package glmmTMB
.
I got the following answer by Ben Bolker in the R mixed-models mailing list:
"You can't constrain fixed-effect parameters in lme4::lmer (they are profiled out of the likelihood function, so they're not accessible to be constrained), but in glmmTMB you could do this"
data("sleepstudy", package= "lme4")
ss <- within(sleepstudy, {
Days <- factor(Days)
contrasts(Days) <- MASS::contr.sdif
})
library(glmmTMB)
m1 <- glmmTMB(Reaction ~ Days + (1 | Subject), data = ss)
m2 <- update(m1,
control = glmmTMBControl(
optArgs = list(lower=c(-Inf, rep(0,9), rep(-Inf, 2)))))