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