[R-sig-ME] How to obtain monotonic effect of ordered factor predictor in lme4 package

Ben Bolker bbo|ker @end|ng |rom gm@||@com
Fri Jul 30 22:14:43 CEST 2021


  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



More information about the R-sig-mixed-models mailing list