2
$\begingroup$

Let's say I have this simple GLMM model in R: model = glmer(correct ~ treatment + (1|id), data = DATA, family = binomial(link=logit)), where correct is my dependent variable (0s or 1s) and treatment is my experimental variable including 3 distinct levels of my intervention.

I would like to get estimated marginal means using lsmeans(model, pairwise ~ treatment), all working well. However the means and SEs seem to be in logits. These are the marginal means:

             lsmean   SE
Placebo      2.51     0.151
Treatment1   2.50     0.151
Treatment2   2.67     0.152

When I transform the lsmeans using inv.logit(), I can easily get the mean probabilities on the original response scale 0 to 1 (I need this scale for my report). However, when I do the same thing for SEs, I get very weird values ~0.54, which makes no sense and thus suggesting I am doing something incorrectly.

Interestingly, when I compute the same model and marginal means in JASP, which automatically transposes the marginal means and SEs to the original response scale, the means JASP shows are the same as in R after using the inv.logit() transformation. But! And this is may main problem: JASP shows SEs = ~0.01, which seems to be correct.

Obviously, the 0.01 from JASP <<<< inv.logit(0.151) from R. Thus, my question would be, how should I transform the SEs from R (values: 0.151, 0.152) to have the SEs on the original response scale (from 0 to 1)?

Thank you, best, M.

EDIT: would it be correct to estimate SE on the original response scale this way: SE = inv.logit(2.51) - inv.logit(2.51-0.151)? This would be ~0.01, as estimated in JASP.

$\endgroup$
2
  • 2
    $\begingroup$ All you need to do is summarize with type = "response". See the documentation for. summary.emmGrid and the vignettes that come with the emmeans package. $\endgroup$
    – Russ Lenth
    Commented Dec 11, 2020 at 16:23
  • $\begingroup$ (+1) Nice comment, Russ - I used it as the basis for my response. Is it possible to compute pairs of treatments in terms of differences of probabilities (rather than ratios of odds) when using lsmeans? It would be great if it were possible! $\endgroup$ Commented Dec 12, 2020 at 4:44

1 Answer 1

4
$\begingroup$

To elaborate on the comment provided by Russ, here are some toy data:

correct <- c(0,1,1,
         1,0,1,
         1,1,0, 
         0,1,1, 
         1,1,0, 
         0,0,0)
treatment <- rep(c("A","B","A","B","A","B"), each = 3)
id <- rep(1:6, each = 3)

DATA <- data.frame(correct, treatment, id) 
DATA

DATA$id <- factor(DATA$id) 
DATA$treatment <- factor(DATA$treatment)

To keep things simple, assume the treatment variable has only 2 levels: A and B.

Using these data, you can fit your model:

library(lme4)

model = glmer(correct ~ treatment + (1|id),
           data = DATA, family = binomial(link=logit))             

For the purpose of this answer, ignore R's singular fit warning (but for your own data, you should not ignore that warning).

Now you can use the lsmeans() function:

library(emmeans)

emm.1 <- lsmeans(model, "treatment")

If you apply the summary() function as is to emm.1:

summary(emm.1)

you will get, for each treatment level, the log odds of your response variable taking the value 1. (The response variable is the one called correct.) R will warn you that the results provided by summary() are given on the logit (not the response) scale and will label the logits as lsmean. (Logit is the same thing as log odds.) Here is the R output for summary(emm.1):

 treatment lsmean    SE  df asymp.LCL asymp.UCL
 A          0.693 0.707 Inf    -0.693      2.08
 B         -0.223 0.671 Inf    -1.538      1.09

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 

If, however, you apply the summary() function to emm.1 with the option type = "response":

summary(emm.1, type = "response")

that will prompt R to report, for each treatment level, the probability that your response variable takes the value 1. This probability will be clearly labeled as prob in the R output:

 treatment  prob    SE  df asymp.LCL asymp.UCL
 A         0.667 0.157 Inf     0.333     0.889
 B         0.444 0.166 Inf     0.177     0.749

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

If you are interested in pairwise treatment comparisons, then use the command:

emm.2 <- lsmeans(model, pairwise ~ treatment) 

Applying the summary() command as is to emm.2 will allow you to compare pairs of treatment levels in terms of the log odds of your response variable taking the value 1:

summary(emm.2) 

The corresponding R output will look like this:

$lsmeans
 treatment lsmean    SE  df asymp.LCL asymp.UCL
 A          0.693 0.707 Inf    -0.693      2.08
 B         -0.223 0.671 Inf    -1.538      1.09

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast estimate    SE  df z.ratio p.value
 A - B       0.916 0.975 Inf 0.940   0.3472 

Results are given on the log odds ratio (not the response) scale. 

Using the option type="response" with the summary() command for emm.2 will allow you to compare pairs of treatment levels by computing the corresponding odds ratios and the associated confidence intervals:

summary(emm.2, type = "response") 

The R output for this last command will look like this:

$lsmeans
 treatment  prob    SE  df asymp.LCL asymp.UCL
 A         0.667 0.157 Inf     0.333     0.889
 B         0.444 0.166 Inf     0.177     0.749

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

$contrasts
 contrast odds.ratio   SE  df z.ratio p.value
 A / B           2.5 2.44 Inf 0.940   0.3472 

Tests are performed on the log odds ratio scale 

I do not know if it is possible with lsmeans() to compare pairs of treatment levels directly on the probability scale by computing differences in probabilities rather than ratios of odds.

Addendum 1

Russ kindly clarified that it is possible to use lsmeans() to compare treatments via differences of probabilities. Here is the R code that would facilitate that comparison for the toy data example above:

emm.3 <- lsmeans(model,  ~ treatment)

pairs(regrid(emm.3))

Here is the R output for the pairs() command:

>     pairs(regrid(emm.3))
 contrast estimate    SE  df z.ratio p.value
 A - B       0.222 0.228 Inf 0.973   0.3304 

Note that all lsmeans() commands in this answer can be replaced by emmeans() commands.

Addendum 2

Note that the computed log odds, odds ratios and probabilities are conditional quantities which refer to the "typical" subject in your study (where each subject is identified by their id), that is, the subject for whom the random intercept is equal to 0.

$\endgroup$
2
  • 1
    $\begingroup$ Use the regrid() function before doing the comparisons. There are lots of examples in the vignettes $\endgroup$
    – Russ Lenth
    Commented Dec 12, 2020 at 4:54
  • 1
    $\begingroup$ Cool! Regrid worked like a charm - I added an addendum to my response to illustrate how it works. Thank you so much for your help, Russ! I learned a lot from your comments on this thread! $\endgroup$ Commented Dec 12, 2020 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.