R lme4 model: calculating effect size between continuous predictor's max-min value

user_11774257 0
Viewed 154

I'm struggling to calculate an effect size between a continuous predictor's max-min value while using an R lme4 multilevel model.

Simulated data: predictor "x" ranges from 1 to 3

library(tidyverse)
n = 100
a = tibble(y = rep(c("pos", "neg", "neg", "neg"), length.out = n), x = rep(3, length.out = n), group = rep(letters[1:7], length.out = n))
b = tibble(y = rep(c("pos", "pos", "neg", "neg"), length.out = n), x = rep(2, length.out = n), group = rep(letters[1:7], length.out = n))
c = tibble(y = rep(c("pos", "pos", "pos", "neg"), length.out = n), x = rep(1, length.out = n), group = rep(letters[1:7], length.out = n))
d = rbind(a, b)
df = rbind(d, c)
df = df %>% mutate(y = as.factor(y))
df

enter image description here

Model

library("lme4")
m = glmer(
  y ~ x + (x | group), 
  data = df, 
  family = binomial(link = "logit"))

Output

ggpredict(m, "x")

.

# Predicted probabilities of y

x | Predicted |       95% CI
----------------------------
1 |      0.75 | [0.67, 0.82]
2 |      0.50 | [0.44, 0.56]
3 |      0.25 | [0.18, 0.33]

Adjusted for:
* group = 0 (population-level)

I'm failing to calculate the effect size between the predictor's "x" max (3) and min (1) value

My best try

library("emmeans")
emmeans(m, "x", trans = "logit", type = "response", at = list(x = c(1, 3)))

 x response     SE  df asymp.LCL asymp.UCL
 1     0.75 0.0387 Inf     0.667     0.818
 3     0.25 0.0387 Inf     0.182     0.333

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

How can I calculate the effect size with CIs between the predictor's "x" max (3) and min (1) value? The effect size should be in probability scale.

lme4 logistic-regression regression r
1 Answers

I'll try to answer, though I'm still not sure what the question is. I am going to assume that what is wanted is the difference between the two probabilities.

There are a lot of moving parts in the emmeans call shown, so I will proceed in smaller steps. First, let's get estimates of the probabilities in question:

> library(emmeans)
> EMM = emmeans(m, "x", at = list(x = c(1, 3)), type = "response")
> EMM
 x prob     SE  df asymp.LCL asymp.UCL
 1 0.75 0.0387 Inf     0.667     0.818
 3 0.25 0.0387 Inf     0.182     0.333

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

The quickest way to obtain a pairwise comparison is via

> pairs(EMM)
 contrast odds.ratio   SE  df null z.ratio p.value
 1 / 3             9 2.94 Inf    1   6.728  <.0001

Tests are performed on the log odds ratio scale

As stated in the annotations (and also in the documentation, e.g. the vignette on comparisons, when a log or logit transformation is in place, the comparison is shown as a ratio. This happens because the tests are performed on the link (logit) scale, and the difference between logs is the log of a ratio.

If we want the difference between probabilities, it is necessary to create a new object where the primary quantities being estimated are the probabilities, rather than their logits. In emmeans, this may be done via the regrid() function:

> EMMP = regrid(EMM, transform = "response")
> EMMP
 x prob     SE  df asymp.LCL asymp.UCL
 1 0.75 0.0387 Inf     0.674     0.826
 3 0.25 0.0387 Inf     0.174     0.326

Confidence level used: 0.95

This output looks a lot like the summary of EMM; however, all memory of the logit transformation has been erased, thus the confidence intervals are different because they are calculated directly from the SEs of the prob estimates. For more information, see the vignette on transformations. So now if we compare these, we get the difference of the probabilities:

> confint(pairs(EMMP))
 contrast estimate     SE  df asymp.LCL asymp.UCL
 1 - 3         0.5 0.0612 Inf      0.38      0.62

Confidence level used: 0.95

(Note: I wrapped this in confint() so that we woul;d obtain a confidence interval, rather than the default summary of the t ratio and P value.)

This could be accomplished in one line of code as follows:

confint(pairs(emmeans(m, "x", transform = "response", at = list(x = c(1, 3)))))

The transform argument requests that the reference grid be immediately passed to regrid(). Note that the correct argument here is transform = "response", rather than transform = "logit" (that is, specify what you want to end with, not what you started with). The latter undoes, then redoes, the logit transformation, putting you back where you started.

The emmeans package provides a lot of options, and I really do recommend reading the vignettes.

user_39615660
Related