5
$\begingroup$

I have a dataset in which I measure some probability of success using a logistic regression approach. In one of the levels of the predictor variable, all observations were a success (or 1). When I look at the estimates and confidence intervals via the emmeans() function, I get a probability of 1 for this level as well as confidence intervals covering the entire range. My guess is that a more precise calculation of the confidence interval is somehow not possible when all observation for a given group are either zero or one since there is zero variation available and hence the confidence has to cover the entire range... As soon as I swap a 1 for a 0 in the aforementioned level of the predictor, I get a more precise confidence interval.

My questions are: What would be the proper statistical approach to analyze such a dataset, or more more importantly how would one interpret such a confidence interval?

Below is an example dataset to illustrate my question:

d <- structure(list(treat = c("Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "Control", "Control", "Control", "Control", 
"Control", "Control", "A", "A", "A", "A", "A", "A", "A", "A", 
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", 
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", 
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", 
"A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "C", "C", "C", 
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", 
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", 
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", 
"C", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D"), rep = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L), surv = c(1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 
0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 
1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L)), class = "data.frame", row.names = c(NA, 
-240L))


m <- glm(
  surv ~ treat ,
  family = binomial(link = "logit"),
  data = d
)

emmeans(m, specs = ~ treat, type = "response")
#treat       prob         SE  df asymp.LCL asymp.UCL
#A       0.916667 0.03989280 Inf  0.798080  0.968368
#B       0.937500 0.03493856 Inf  0.823371  0.979702
#C       0.979167 0.02061518 Inf  0.866394  0.997073
#Control 0.666667 0.06804138 Inf  0.523242  0.784699
#D       1.000000 0.00000494 Inf  0.000000  1.000000

#Confidence level used: 0.95 
#Intervals are back-transformed from the logit scale
$\endgroup$
5
  • 1
    $\begingroup$ Well, it is a fact of life. A probability of 1 comes out as + infinity on the logit scale, and it's pretty hard to construct a confidence interval around that. If you regrid it first, it will instead use the SE of zero to construct a CzmI of 1 to 2, which is equally ridiculous. $\endgroup$
    – Russ Lenth
    Aug 11, 2022 at 23:20
  • $\begingroup$ Thanks Russ! If you want to draft up a quick answer I will accept it. @RussLenth $\endgroup$
    – Stefan
    Aug 12, 2022 at 0:59
  • $\begingroup$ OK I did so. Sorry about the typos in the comment - I have no idea where "Czml" came from but that's a fact of fat fingers and a phone. And it should say 1 to 1. $\endgroup$
    – Russ Lenth
    Aug 13, 2022 at 14:34
  • $\begingroup$ You could perform a Bayesian regression in which your credible interval is not only influenced by the data but also by your prior. So that would open the door to a sensible credible interval but only if you can define a sensible prior. Depending on the application that may be helpful or not. $\endgroup$
    – Bernhard
    Aug 13, 2022 at 14:42
  • $\begingroup$ @Bernhard True. Bayesian methods are generally a lot cleaner when we get close a boundary of the parameter space. The Agresti-Caffo metho I mention is basically using a beta(2,2) prior. $\endgroup$
    – Russ Lenth
    Aug 13, 2022 at 18:56

2 Answers 2

7
$\begingroup$

Well, this is kind of a fact of life for logit models. A probability of 1 comes out as $+\infty$ on the logit scale, and it's pretty hard to construct a confidence interval around that.

If you regrid it first, e.g. regrid(emmeans(m, "treat")), this resets the whole structure on the response (probability) scale with no record that it was originally on the logit scale. You'll find then that the estimate will be 1.0 with an SE of 0.0. It then uses that zero SE to construct a confidence interval of 1.0 to 1.0 -- this of course is equally ridiculous.

There is a 1999 article by Agresti and Caffo that suggests, in a simple scenario, that tossing-in two extra successes and two extra failures leads to better results (though biased towards 0.5). In your example, that approach is possible:

dd = data.frame(
    treat = factor(c("control", "A", "B", "C", "D"),
                   levels = c("control", "A", "B", "C", "D")),
    succ = c(32, 44, 45, 47, 48),
    fail = c(16, 4, 3, 1, 0))

dd.glm = glm(cbind(succ + 2, fail + 2) ~ treat, data = dd, 
         family = binomial)

emmeans(dd.glm, "treat", type = "response")
##  treat    prob     SE  df asymp.LCL asymp.UCL
##  control 0.654 0.0660 Inf     0.516     0.770
##  A       0.885 0.0443 Inf     0.766     0.947
##  B       0.904 0.0409 Inf     0.789     0.959
##  C       0.942 0.0323 Inf     0.836     0.981
##  D       0.962 0.0267 Inf     0.859     0.990
##
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale
$\endgroup$
1
  • $\begingroup$ That's great! Thanks! $\endgroup$
    – Stefan
    Aug 30, 2022 at 22:19
3
$\begingroup$

This just in... Try using the logistf package instead of glm:

library(logistf)
> mf <- logistf(
+     surv ~ treat ,
+     family = binomial(link = "logit"),
+     data = d
+ )

> emmeans(mf, specs = ~ treat, type = "response")
 treat    prob     SE  df lower.CL upper.CL
 A       0.908 0.0413 235    0.789    0.963
 B       0.929 0.0368 235    0.813    0.975
 C       0.969 0.0246 235    0.861    0.994
 Control 0.663 0.0675 235    0.521    0.781
 D       0.990 0.0144 235    0.855    0.999

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

Note:

As this is written, emmeans support has just now been added to a branch of the logistf package, and it will take a while to reach CRAN. But in the meantime, you can use qdrg() instead as follows:

> mf.rg <- emmeans::qdrg(object = mf, data = d, link = "logit")
> emmeans(mf.rg, specs = ~ treat, type = "response")
 
    (same output)
$\endgroup$
1
  • $\begingroup$ Thanks Russ for following up on this! Much appreciated! $\endgroup$
    – Stefan
    Mar 16, 2023 at 20:47

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.