[R-sig-ME] How to obtain monotonic effect of ordered factor predictor in lme4 package
Cristobal Moya
cr|@tob@|moy@ @end|ng |rom gm@||@com
Tue Aug 10 10:49:21 CEST 2021
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 using 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 using 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 using 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 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. Benjamin Bolker
> Professor, Mathematics & Statistics and Biology, McMaster University
> Director, School of Computational Science and Engineering
> Graduate chair, Mathematics & Statistics
>
> _______________________________________________
> 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