emmeans: Is bias adjustment appropriate?

Question

I have the following model. The outcome, prop_correct, is a proportion of words correctly identified by participants in 5-word sentences (0 - 1) which is logit-transformed. Fixed effects include speaker group (disordered vs. non-disordered), Condition 1 (levels = on, off), and Condition 2 (levels = A, B, C; contrasts: A vs B+C, B vs C), and all possible interactions (all speakers participated in both conditions). Random effects include 1) by-participant intercept, 2) by-speaker intercepts and slopes for the two Condition 2 contrasts, and 3) by-sentence intercepts.

mod <- lmerTest::lmer(
  car::logit(prop_correct) ~
                    group*cond1*cond2 +
                    (1 | participant)+
                    (1 + (cond1AvBC + cond1BvC) | speaker)+
                    (1|sentence),
                  data=df)

This produces a three-way interaction, and I am using emmeans to investigate the difference between the two levels of Condition 1 for each group and Condition 2 level.

Given the discussion of bias adjustment in https://cran.r-project.org/web/packages/emmeans/vignettes/transformations.html#bias-adj, I assumed this was the most appropriate approach. However, the empirical means are much closer to the non-bias-adjusted emmeans estimates. Does this imply that I should use the non-bias-adjusted means? And, if so, why is this, given this is a (relatively complex) mixed model?

emmeans with bias adjustment:

sigma <- sqrt(sum(as.data.frame(lme4::VarCorr(mod))$vcov)) # 2.24
emm_biasadj <- emmeans::emmeans(mod, pairwise ~ cond1|cond2|group,
                                 tran = "logit",
                                 type = "response",
                                 bias.adj = TRUE,
                                 sigma = sigma,
                               adjust = "bonferroni"
                               )

Output:

cond2 = A, group = non-dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.707 0.06651 Inf     0.822     0.575
 on           0.644 0.06881 Inf     0.772     0.524

cond2 = B, group = non-dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.829 0.03511 Inf     0.888     0.750
 on           0.697 0.04953 Inf     0.787     0.599

cond2 = C, group = non-dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.771 0.06212 Inf     0.871     0.635
 on           0.659 0.07305 Inf     0.791     0.528

cond2 = A, group = dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.490 0.02545 Inf     0.583     0.488
 on           0.503 0.02127 Inf     0.479     0.519

cond2 = B, group = dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.688 0.04996 Inf     0.779     0.590
 on           0.480 0.00684 Inf     0.520     0.489

cond2 = C, group = dis:
 cond1 response      SE  df asymp.LCL asymp.UCL
 off           0.576 0.06658 Inf     0.716     0.485
 on           0.493 0.03030 Inf     0.599     0.488

Degrees-of-freedom method: asymptotic 
Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 
Bias adjustment applied based on sigma = 2.2404 

emmeans without bias adjustment:

emm_nobiasadj <- emmeans::emmeans(mod, pairwise ~ cond1|cond2|group,
                                 tran = "logit",
                                 type = "response",
                               adjust = "bonferroni"
                               )

Output:

cond2 = A, group = non-dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.894 0.0319 Inf     0.813     0.942
 on           0.861 0.0403 Inf     0.762     0.923

cond2 = B, group = non-dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.945 0.0129 Inf     0.914     0.966
 on           0.889 0.0244 Inf     0.832     0.929

cond2 = C, group = non-dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.923 0.0255 Inf     0.855     0.960
 on           0.870 0.0406 Inf     0.768     0.931

cond2 = A, group = dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.702 0.0704 Inf     0.549     0.820
 on           0.489 0.0841 Inf     0.331     0.649

cond2 = B, group = dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.885 0.0253 Inf     0.825     0.926
 on           0.658 0.0556 Inf     0.543     0.758

cond2 = C, group = dis:
 cond1 response     SE  df asymp.LCL asymp.UCL
 off           0.814 0.0540 Inf     0.685     0.898
 on           0.711 0.0733 Inf     0.550     0.832

Degrees-of-freedom method: asymptotic 
Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

Comparison of outputs with empirical means

|cond1 |cond2   |group   | empirical   | emm_adj   | emm_nonadj|
|:-----|:-------|:-------|------------:|----------:|----------:|
|on    |B       |non-dis |    0.8236285|  0.6970182|  0.8893415|
|on    |A       |non-dis |    0.7406944|  0.6443800|  0.8611386|
|on    |C       |non-dis |    0.7741667|  0.6590652|  0.8695067|
|off   |B       |non-dis |    0.9149006|  0.8293237|  0.9451524|
|off   |A       |non-dis |    0.8235185|  0.7066464|  0.8940260|
|off   |C       |non-dis |    0.8639120|  0.7709109|  0.9225107|
|on    |B       |dis     |    0.6233873|  0.4795816|  0.6583957|
|on    |A       |dis     |    0.4763141|  0.5028886|  0.4886361|
|on    |C       |dis     |    0.6390031|  0.4932362|  0.7106518|
|off   |B       |dis     |    0.7937821|  0.6876523|  0.8846596|
|off   |A       |dis     |    0.6522436|  0.4899969|  0.7022852|
|off   |C       |dis     |    0.7312788|  0.5757685|  0.8143192|

The mean absolute difference between the empirical vs. bias-adjusted means is 0.114, vs 0.065 for the empirical vs. non-bias-adjusted.

Based on this, I think it is appropriate to use the non-bias-adjusted emmeans, but I don't know for certain and I don't understand why.

In an answer to a related question here, Russ Lenth did advise to use bias-adjustment, throwing into question whether the non-bias adjusted means are actually more appropriate in my example.

Answer 1

I think that if a bias adjustment is appropriate, then it is appropriate; and if it is not, then it is not. That choice should not be based on how the estimates compare.

There are a number of issues to consider here. One biggie is that not all of the random effects in this example are random intercepts, as would be necessary for how sigma is computed. If you had only random intercepts, then they can be combined into a total SD via a Pythagorean sum as shown. But the term (1 + (cond1AvBC + cond1BvC) | speaker) includes additional random slopes that are either not accounted for, or not correctly accounted for in that computation; but those random effects provide different contributions to the variance depending on the level of cond1 (though from the description, it is cond2 that has levels A, B, and C; so I am confused as to what those predictors mean).

I will also mention that the linked article uses a generalized linear mixed model with a logit link, while this one uses a linear mixed model with logit response transformation. That makes a difference too, because in this case, the residual variance should be included in sigma, whereas in the GLMM it should not.

I am not sure what to suggest. It is difficult to do the bias adjustments with these random effects. If it is really important to predict proportion correct, as opposed to logit(prop correct), perhaps it would be better to omit the logit response transformation and model those proportions directly -- perhaps devising appropriate weights for each response when fitting the model.