Skip to content

How to obtain monotonic effect of ordered factor predictor in lme4 package

3 messages · Cristobal Moya, Henrik Singmann, Ben Bolker

Original
Cristobal Moya
# 
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
2 days later
Henrik Singmann
# 
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>:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
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: