How to obtain monotonic effect of ordered factor predictor in lme4 package
Dear Henrik and Ben, Indeed, I hadn't seen Ben's email. Thank you very much for showing this option with glmmTMB, it provides much similar estimates compared to the brms mo() parameterization. Kind regards, Crist?bal
On Tue, Aug 10, 2021 at 7:26 AM Henrik Singmann <singmann at gmail.com> wrote:
Hi Crist?bal,
It looks like you did not get Ben's email (which only went to the list)
which I attach again below. It clarifies that monotonic effects in the way
implemented in brms are not possible in lme4 but in glmmTMB.
Cheers,
Henrik
---------- Forwarded message ---------
Von: Ben Bolker <bbolker at gmail.com>
Date: Fr., 30. Juli 2021 um 22:15 Uhr
Subject: Re: [R-sig-ME] How to obtain monotonic effect of ordered factor
predictor in lme4 package
To: <r-sig-mixed-models at r-project.org>
A few more thoughts
codingMatrices::code_diff() and MASS::contr.sdif() are basically
identical (in case you don't feel like installing another package --
MASS is automatically installed with base R)
Not quite sure why you're using code_diff_forward rather than code_diff?
The latter would naively make more sense to me.
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)))))
Notes: (1) the constraints aren't binding in this example because all
the parameters are estimated as being positive anyway (I was too lazy to
find another example/perturb this example); (2) you do need to know the
parameter ordering for this to work: in this case we have the intercept
(unconstrained) followed by 9 successive difference values, followed by
parameters for the RE variance and the dispersion (estimated on the log
scale, unconstrained). The combination of fixef(m2) and m2$obj$par can
help you figure this out.
On 7/30/21 8:14 AM, Henrik Singmann wrote:
Hi Crist?bal, As far as I know, there is no correspondence in lme4 to the monotonic effects as there is no way to enforce certain fixed effect estimate are strictly positive or negative. However, if the data is strictly ordered, you can get something similar
by
using appropriate coding. So instead of using an ordered factor you use a
regular factor and code_diff_forward() from package codingMatrices.
library("lme4")
library("codingMatrices")
sleepstudy$Days <- factor(sleepstudy$Days)
m2 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy,
contrasts = list(Days = "code_diff_forward"))
summary(m2)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (1 | Subject)
#> Data: sleepstudy
#>
#> REML criterion at convergence: 1729.5
#>
#> 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.5079 9.0499 32.985
#> Days0-1 -7.8439 10.4753 -0.749
#> Days1-2 -0.8661 10.4753 -0.083
#> Days2-3 -17.6301 10.4753 -1.683
#> Days3-4 -5.6574 10.4753 -0.540
#> Days4-5 -19.8690 10.4753 -1.897
#> Days5-6 -3.6598 10.4753 -0.349
#> Days6-7 -6.5723 10.4753 -0.627
#> Days7-8 -17.8789 10.4753 -1.707
#> Days8-9 -14.2217 10.4753 -1.358
#> [...]
This give you the differences between each two days as shown next:
library("tidyverse")
sleepstudy %>%
group_by(Days) %>%
summarise(mean = mean(Reaction)) %>%
mutate(meanlag = lag(mean)) %>%
mutate(diff = meanlag-mean)
#> # A tibble: 10 x 4
#> Days mean meanlag diff
#> <fct> <dbl> <dbl> <dbl>
#> 1 0 257. NA NA
#> 2 1 264. 257. -7.84
#> 3 2 265. 264. -0.866
#> 4 3 283. 265. -17.6
#> 5 4 289. 283. -5.66
#> 6 5 309. 289. -19.9
#> 7 6 312. 309. -3.66
#> 8 7 319. 312. -6.57
#> 9 8 337. 319. -17.9
#> 10 9 351. 337. -14.2
There is also the possibility to get exactly the same results by just
using
emmeans on the model fitted with default contrasts:
library("lme4")
sleepstudy$Days <- factor(sleepstudy$Days)
m2 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
library("emmeans")
emmeans(m2, "Days", contr = "consec")
#> $emmeans
#> Days emmean SE df lower.CL upper.CL
#> 0 257 11.5 42 234 280
#> 1 264 11.5 42 241 288
#> 2 265 11.5 42 242 288
#> 3 283 11.5 42 260 306
#> 4 289 11.5 42 266 312
#> 5 309 11.5 42 285 332
#> 6 312 11.5 42 289 335
#> 7 319 11.5 42 296 342
#> 8 337 11.5 42 314 360
#> 9 351 11.5 42 328 374
#>
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> 1 - 0 7.844 10.5 153 0.749 0.9887
#> 2 - 1 0.866 10.5 153 0.083 1.0000
#> 3 - 2 17.630 10.5 153 1.683 0.5279
#> 4 - 3 5.657 10.5 153 0.540 0.9988
#> 5 - 4 19.869 10.5 153 1.897 0.3775
#> 6 - 5 3.660 10.5 153 0.349 1.0000
#> 7 - 6 6.572 10.5 153 0.627 0.9965
#> 8 - 7 17.879 10.5 153 1.707 0.5098
#> 9 - 8 14.222 10.5 153 1.358 0.7629
#>
#> Degrees-of-freedom method: kenward-roger
#> P value adjustment: mvt method for 9 tests
Of course, all of this only works if the means are ordered in this way
(otherwise you will have effects that differ in sign).
Hope that helps,
Henrik
Am Di., 27. Juli 2021 um 19:06 Uhr schrieb Cristobal Moya <
cristobalmoya at gmail.com>:
Dear list members, I have a question regarding monotonic effects with lme4. I also posted
this
in StackOverflow (https://stackoverflow.com/q/68546489/6832021), which
is
what I reproduce below. 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. 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 in `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
...
How could an equivalent monotonic effect of `moDays` in `m1` can be
obtained with lme4?
I'm grateful for anyone who can provide some orientation. Kind regards,
--
Crist?bal Moya
Research Fellow
Chair for Social Structure Analysis of Social Inequalities
Faculty of Sociology
Bielefeld University
[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering Graduate chair, Mathematics & Statistics
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models -- Dr. Henrik Singmann Lecturer, Experimental Psychology University College London (UCL), UK http://singmann.org