[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