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

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?

Here is a reproducible example:

require(multcomp)
#>
#> 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") -
predict(id.glm3, newdata = pred_detec3, type = "link")
#>          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") -
predict(id.glm3, newdata = pred_detec3, type = "link"))
#>         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") /
predict(id.glm3, newdata = pred_detec3, type = "response")
#>         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

• @stats_noob on your 1st question. Yes, I agree with the post you cited. I would use confint(glht4) and then exponentiate the endpoints of the interval. Oct 24, 2022 at 13:27