Interpreting odds ratio of multiple comparisons from a logistic regression model (using R)

Question

I recently ran a logistic regression on categorical data and ran a Tukey multiple comparisons post hoc analysis using the glht function in multcomp package.

id.glm3<-glm(detec~apptreat+marker+exp+inter_MarEx, family=binomial, data=indiv_detec3)

Where inter_MarEx is the interaction between marker type (marker) and time after exposure (exp)

I understand that I must convert the coefficients from log odds ratios to odds ratios via exp(coefficient) and generally understand how to interpret them with mixed continuous. However, I do not quite understand how to interpret the odds ratio resulting from the multiple comparisons.

summary(glht(id.glm4, linfct=mcp(inter_MarEx="Tukey")))

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = detec ~ apptreat + inter_MarEx, family = binomial, 
data = indiv_detec3)

Linear Hypotheses:
                       Estimate Std. Error z value Pr(>|z|)    
EW 24 - EW 0 == 0       -1.5989     0.6196  -2.581   0.0937 .  
EW 48 - EW 0 == 0        0.9306     0.6190   1.503   0.6429   
EW 48 - EW 24 == 0       2.5295     0.3575   7.076   <0.001 *** 
Milk 24 - Milk 0 == 0    0.3010     0.4442   0.678   0.9828    
Milk 48 - Milk 0 == 0    4.0374     0.5190   7.779   <0.001 ***
Milk 48 - Milk 24 == 0   3.7364     0.4681   7.983   <0.001 ***
Milk 0 - EW 0 == 0       0.8109     0.6180   1.312   0.7632    
Milk 24 - EW 24 == 0     2.7108     0.2716   9.980   <0.001 ***
Milk 48 - EW 48 == 0     3.9177     0.5161   7.592   <0.001 ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Now I convert the coefficients to get the odds ratio from log odds ratio.

Linear Hypotheses:

                   Estimate Std. Error z value Pr(>|z|)    
EW 24 - EW 0 == 0        0.20     0.6196  -2.581   0.0937 .  
EW 48 - EW 0 == 0        2.54     0.6190   1.503   0.6429    
EW 48 - EW 24 == 0       12.55     0.3575   7.076   <0.001 ***
Milk 24 - Milk 0 == 0    1.35     0.4442   0.678   0.9828    
Milk 48 - Milk 0 == 0    56.68     0.5190   7.779   <0.001 ***
Milk 48 - Milk 24 == 0   41.95     0.4681   7.983   <0.001 ***
Milk 0 - EW 0 == 0       2.25     0.6180   1.312   0.7632    
Milk 24 - EW 24 == 0     15.04     0.2716   9.980   <0.001 ***
Milk 48 - EW 48 == 0     50.28     0.5161   7.592   <0.001 ***

How would I go about interpreting the odds ratios for the multiple comparisons above? For example, when comparing Milk 48 and EW 48 (odds ratio = 50.28), would I interpret the output as saying the odds of detecting EW at 48 hours after exposure are 50:1 greater than the odds of detecting Milk at 48 hrs after exposure? If this is incorrect, how should this be interpreted?

Answer 1

Here is a reproducible example:

require(multcomp)
#> Loading required package: multcomp
#> Loading required package: mvtnorm
#> Loading required package: survival
#> Loading required package: TH.data
#> Loading required package: MASS
#> 
#> Attaching package: 'TH.data'
#> The following object is masked from 'package:MASS':
#> 
#>     geyser

set.seed(102393)
N <- 200
indiv_detec3 <- data.frame(
  marker = factor(rbinom(N, 1, prob = c(0.5)), labels = c("EW", "Milk")),
  exp = factor(sample(c(0, 24, 48), size = N, replace = TRUE)),
  apptreat = sample(c("A", "B", "C"), size = N, replace = TRUE),
  detec = rbinom(N, 1, prob = c(0.7))
)

I recommend using the exp*marker notation since it is more explicit than making your own variable. Of course, for the glht procedure, you have to make your own.

id.glm3 <- glm(detec ~ apptreat + exp*marker, 
               family = binomial(link = "logit"), data = indiv_detec3)

summary(id.glm3)
#> 
#> Call:
#> glm(formula = detec ~ apptreat + exp * marker, family = binomial(link = "logit"), 
#>     data = indiv_detec3)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.2971   0.3763   0.5964   0.8003   0.9038  
#> 
#> Coefficients:
#>                  Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)       2.27104    0.78019   2.911  0.00360 **
#> apptreatB         0.07351    0.43732   0.168  0.86652   
#> apptreatC         0.34134    0.42393   0.805  0.42072   
#> exp24            -1.51956    0.81360  -1.868  0.06180 . 
#> exp48            -0.97584    0.84935  -1.149  0.25059   
#> markerMilk       -1.54518    0.82925  -1.863  0.06241 . 
#> exp24:markerMilk  3.01646    1.16348   2.593  0.00952 **
#> exp48:markerMilk  0.93434    0.99308   0.941  0.34678   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 213.27  on 199  degrees of freedom
#> Residual deviance: 202.28  on 192  degrees of freedom
#> AIC: 218.28
#> 
#> Number of Fisher Scoring iterations: 5

Make the interaction for creating confidence intervals:

indiv_detec3$inter_MarEx <- with(indiv_detec3, interaction(exp, marker))

id.glm4 <- glm(detec ~ apptreat + inter_MarEx, family = binomial(link = "logit"), data = indiv_detec3)

summary(id.glm4)
#> 
#> Call:
#> glm(formula = detec ~ apptreat + inter_MarEx, family = binomial(link = "logit"), 
#>     data = indiv_detec3)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.2971   0.3763   0.5964   0.8003   0.9038  
#> 
#> Coefficients:
#>                    Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)         2.27104    0.78019   2.911   0.0036 **
#> apptreatB           0.07351    0.43732   0.168   0.8665   
#> apptreatC           0.34134    0.42393   0.805   0.4207   
#> inter_MarEx24.EW   -1.51956    0.81360  -1.868   0.0618 . 
#> inter_MarEx48.EW   -0.97584    0.84935  -1.149   0.2506   
#> inter_MarEx0.Milk  -1.54518    0.82925  -1.863   0.0624 . 
#> inter_MarEx24.Milk -0.04828    1.04667  -0.046   0.9632   
#> inter_MarEx48.Milk -1.58669    0.81982  -1.935   0.0529 . 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 213.27  on 199  degrees of freedom
#> Residual deviance: 202.28  on 192  degrees of freedom
#> AIC: 218.28
#> 
#> Number of Fisher Scoring iterations: 5

Finally, using the multiple comparisons procedure:

glht4 <- multcomp::glht(id.glm4, linfct = multcomp::mcp(inter_MarEx = "Tukey"))
summary(glht4)
#> 
#>   Simultaneous Tests for General Linear Hypotheses
#> 
#> Multiple Comparisons of Means: Tukey Contrasts
#> 
#> 
#> Fit: glm(formula = detec ~ apptreat + inter_MarEx, family = binomial(link = "logit"), 
#>     data = indiv_detec3)
#> 
#> Linear Hypotheses:
#>                        Estimate Std. Error z value Pr(>|z|)
#> 24.EW - 0.EW == 0      -1.51956    0.81360  -1.868    0.406
#> 48.EW - 0.EW == 0      -0.97584    0.84935  -1.149    0.852
#> 0.Milk - 0.EW == 0     -1.54518    0.82925  -1.863    0.409
#> 24.Milk - 0.EW == 0    -0.04828    1.04667  -0.046    1.000
#> 48.Milk - 0.EW == 0    -1.58669    0.81982  -1.935    0.365
#> 48.EW - 24.EW == 0      0.54372    0.54280   1.002    0.912
#> 0.Milk - 24.EW == 0    -0.02562    0.51123  -0.050    1.000
#> 24.Milk - 24.EW == 0    1.47128    0.81733   1.800    0.449
#> 48.Milk - 24.EW == 0   -0.06713    0.49713  -0.135    1.000
#> 0.Milk - 48.EW == 0    -0.56934    0.56523  -1.007    0.910
#> 24.Milk - 48.EW == 0    0.92756    0.85277   1.088    0.879
#> 48.Milk - 48.EW == 0   -0.61084    0.55028  -1.110    0.870
#> 24.Milk - 0.Milk == 0   1.49690    0.83243   1.798    0.450
#> 48.Milk - 0.Milk == 0  -0.04150    0.51570  -0.080    1.000
#> 48.Milk - 24.Milk == 0 -1.53840    0.82214  -1.871    0.403
#> (Adjusted p values reported -- single-step method)

Note: in the original poster's question, not all pairwise comparisons are shown.

pred_detec3 <- data.frame(
  marker = factor(c(1, 0), labels = c("EW", "Milk")),
  exp = factor(c(48, 48), levels = c(0, 24, 48)),
  apptreat = factor(c("A","A"), levels = c("A", "B", "C"))
)

pred_detec3$inter_MarEx <- with(pred_detec3, interaction(exp, marker))

$$P_{milk,48,A} = P(Detect = 1 | marker = Milk, exposure = 48, apptreat = A) = 0.665$$

$$P_{EW,48,A} = P(Detect = 1 | marker = EW, exposure = 48, apptreat = A) = 0.785$$

plogis(predict(id.glm3, newdata = pred_detec3, type = "link"))
#>         1         2 
#> 0.6647106 0.7850259
predict(id.glm3, newdata = pred_detec3, type = "response")
#>         1         2 
#> 0.6647106 0.7850259

Difference in log odds: (same as the glht procedure 48.Milk - 48.EW == 0)

predict(id.glm3, newdata = pred_detec3, type = "link")[1] -
  predict(id.glm3, newdata = pred_detec3, type = "link")[2]
#>          1 
#> -0.6108418

Odds ratio:

$$\large \frac{\frac{P_{milk,48,A}}{1-P_{milk,48,A}}}{\frac{P_{EW,48,A}}{1-P_{EW,48,A}}} = 0.543$$

exp(predict(id.glm3, newdata = pred_detec3, type = "link")[1] -
      predict(id.glm3, newdata = pred_detec3, type = "link")[2])
#>         1 
#> 0.5428937

The ratio of probabilities: $$\frac{P_{milk,48,A}}{P_{EW,48,A}} = 0.847$$

predict(id.glm3, newdata = pred_detec3, type = "response")[1] /
  predict(id.glm3, newdata = pred_detec3, type = "response")[2]
#>         1 
#> 0.8467371

Original Question

Yes, $50.28 =$ the odds ratio of a detection given Milk and 48 relative to a detection given EW and 48 with the same AppTreat

Another note on the model:

The original posting worded it as "odds of detecting EW at 48 hours". If that is the case, then the model might need to be changed. The model that was fit was:

$$Detect = AppTreat + Marker*Exposure$$

This type of model says that the Marker is an independent variable that is "set" in the experiment. It also implies that the type of detection is the same if it is a Milk Marker or EW marker. Therefore, it is incorrect to say "milk detection" or "EW detection". If, on the other hand, what was intended was two types of detection, Milk and EW, then the model might need to be changed.

$${No\ detection,\ Milk\ detection,\ EW\ detection,\ Both?} = AppTreat + Exposure$$

Update: comment response

The reason that the apptreat variable has to be the same in the odds ratio is that you are trying to isolate the effects in the contrast.

This notation is loose, but it makes the point:

$$ln\frac{P(Y=1)}{1-P(Y=1)} = \beta_0 + \beta_1 apptreat + \beta_2 interMarEx$$

For this contrast: 48.Milk - 48.EW

$$\beta_0 + \beta_{1,A} + \beta_{2,48,Milk} - (\beta_0 + \beta_{1,A} + \beta_{2,48,EW}) = \beta_{2,48,Milk} - \beta_{2,48,EW} = ln \left(\frac{P_{milk,48,A}}{1-P_{milk,48,A}}\right) - ln \left(\frac{P_{EW,48,A}}{1-P_{EW,48,A}}\right)$$

If you did not have a consistent $\beta_{1,A}$ in both options, then they wouldn't cancel and create the odds ratio of interest.

^{Created on 2022-10-23 with reprex v2.0.2}