[R-sig-ME] How to obtain monotonic effect of ordered factor predictor in lme4 package
Henrik Singmann
@|ngm@nn @end|ng |rom gm@||@com
Fri Jul 30 14:14:01 CEST 2021
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 using 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 using 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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list