1
$\begingroup$

I am giving a talk on logistic regression and I was going to mention log-binomial models to estimate risk-ratios. I understand the difference between odds and probabilities and that they only converge when the outcome is rare. I also understand that the unit change for the odds ratio is constant across all values of a continuous predictor, whereas probabilities will depend on the value of the predictor.

I obviously hadn't fully reconciled all of that in my mind though, in that I was expecting that when I predict probabilities from either a logistic or log-binomial model, they should still be the same. Furthermore, I had anticipated that (in a logistic model) as the probabilities change, so to will the risk ratios across the values of the predictor. But of course, with a log-binomial model, it outputs a coefficient that can be converted to a risk ratio which also remains constant across the values of the predictor.

Do I just need to accept that the predicted probabilities from a logistic and log-binomial model will be different, with the former giving constant odds-ratios and the latter giving constant risk-ratios?

I have tried to illustrate this with the code below.

library(tidyverse)
library(emmeans)

# Simulate data
n <- 100                    
set.seed(1234)
x  <-  rnorm(n)    
z  <-  -x     
pr  <-  1/(1 + exp(-z))      
y  <-  rbinom(n, 1, pr)      
df <-  data.frame(y = y, x = x, z = z, pr = pr)
df
table(df$y)

# Logistic model
summary(m1 <- glm(y ~ x, data = df, 
                  family = binomial(link = "logit")))
# Probabilities at set values of x
emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response")
# Pairwise odds ratios 
emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    pairs(rev = TRUE)
# Probabilities at set values of x using log transform
# Does this emulate predicted risks from log-binomial model? 
   # No - as these are the same as above.
emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    regrid("log")
# Pairwise risk ratios?
emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    regrid("log") |> pairs(rev = TRUE)

# Log-binomial model
summary(m2 <- glm(y ~ x, data = df, family = binomial(link = "log"), 
    start = c(-1, 0)))
# Probabilities at set values of x
emmeans(m2, ~ x, at = list(x = c(0, 1, 2)), type = "response")
# Pairwise risk ratios 
emmeans(m2, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    pairs(rev = TRUE)

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Output from running above

> table(df$y)

 0  1 
44 56 
> 
> # Logistic model
> summary(m1 <- glm(y ~ x, data = df, 
    family = binomial(link = "logit")))

Call:
glm(formula = y ~ x, family = binomial(link = "logit"), data = df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.07137    0.23304   0.306    0.759    
x           -1.16432    0.28234  -4.124 3.73e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.19  on 99  degrees of freedom
Residual deviance: 113.45  on 98  degrees of freedom
AIC: 117.45

Number of Fisher Scoring iterations: 4

> # Probabilities at set values of x
> emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response")
 x   prob     SE  df asymp.LCL asymp.UCL
 0 0.5178 0.0582 Inf    0.4048     0.629
 1 0.2511 0.0755 Inf    0.1323     0.424
 2 0.0947 0.0561 Inf    0.0282     0.274

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 
> # Pairwise odds ratios 
> emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    pairs(rev = TRUE)
 contrast odds.ratio     SE  df null z.ratio p.value
 x1 / x0      0.3121 0.0881 Inf    1  -4.124  0.0001
 x2 / x0      0.0974 0.0550 Inf    1  -4.124  0.0001
 x2 / x1      0.3121 0.0881 Inf    1  -4.124  0.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 
> # Probabilities at set values of x using log transform
> # Does this emulate predicted risks from log-binomial model? No - as these are the same as above.
> emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    regrid("log")
 x response     SE  df asymp.LCL asymp.UCL
 0   0.5178 0.0582 Inf    0.4155     0.645
 1   0.2511 0.0755 Inf    0.1392     0.453
 2   0.0947 0.0561 Inf    0.0297     0.302

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 
> # Pairwise risk ratios?
> emmeans(m1, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    regrid("log") |> pairs(rev = TRUE)
 contrast ratio     SE  df null z.ratio p.value
 x1 / x0  0.485 0.1127 Inf    1  -3.114  0.0052
 x2 / x0  0.183 0.0989 Inf    1  -3.143  0.0048
 x2 / x1  0.377 0.1167 Inf    1  -3.151  0.0046

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log scale 
> 
> # Log-binomial model
> summary(m2 <- glm(y ~ x, data = df, 
    family = binomial(link = "log"), start = c(-1, 0)))

Call:
glm(formula = y ~ x, family = binomial(link = "log"), data = df, 
    start = c(-1, 0))

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.70543    0.10700  -6.593 4.32e-11 ***
x           -0.30073    0.04562  -6.593 4.32e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.19  on 99  degrees of freedom
Residual deviance: 117.87  on 98  degrees of freedom
AIC: 121.87

Number of Fisher Scoring iterations: 25

There were 27 warnings (use warnings() to see them)
> # Probabilities at set values of x
> emmeans(m2, ~ x, at = list(x = c(0, 1, 2)), type = "response")
 x  prob     SE  df asymp.LCL asymp.UCL
 0 0.494 0.0528 Inf     0.400     0.609
 1 0.366 0.0558 Inf     0.271     0.493
 2 0.271 0.0537 Inf     0.184     0.399

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 
> # Pairwise risk ratios 
> emmeans(m2, ~ x, at = list(x = c(0, 1, 2)), type = "response") |> 
    pairs(rev = TRUE)
 contrast ratio     SE  df null z.ratio p.value
 x1 / x0  0.740 0.0338 Inf    1  -6.593  <.0001
 x2 / x0  0.548 0.0500 Inf    1  -6.593  <.0001
 x2 / x1  0.740 0.0338 Inf    1  -6.593  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log scale 
$\endgroup$
3
  • $\begingroup$ I think it should work in the case where you are not relying on linearity (e.g. where you have only one categorical predictor, or where you have only categorical predictors with all interactions among them included). $\endgroup$
    – Ben Bolker
    Dec 14, 2023 at 0:51
  • $\begingroup$ Thanks Ben - I think that might be the answer. I changed x to categorical (if I've simulated it right) x <- sample(1:3, 100, replace = TRUE, prob=c(0.2, 0.5, 0.3)), and now the predicted probs seem to align for both logistic and log-binomial output. The ratios remain different, with OR's further from 1 than RR's, as expected. $\endgroup$
    – LucaS
    Dec 14, 2023 at 1:22
  • $\begingroup$ This sounds like an answer. Does either of you want to post it as such? Comments are fleeting. @BenBolker $\endgroup$ Dec 14, 2023 at 7:56

1 Answer 1

1
$\begingroup$

As per Ben Bolker's comment, changing x to categorical - e.g. x <- sample(1:3, 100, replace = TRUE, prob=c(0.2, 0.5, 0.3)), results in predicted probs that align for both logistic and log-binomial output.

$\endgroup$

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.