I am using logistic regression to model the effect of two production variables on the occurrence of tuberculosis. My independent variables are A (farm size: 0 for small, 1 for large) and B (cattle co-occurrence: 0 for no cattle, 1 for cattle presence). Reference levels for the veriables are:
- Farm size = 0
- Cattle presence = 0
The final model includes the interaction term A x B. Before presenting my main question, please see model output below:
global <- glm(TB ~ Farm size + Cattle + Farm size * Cattle, data = datos, family = binomial) summary(global) Deviance Residuals: Min 1Q Median 3Q Max -2.2101 -0.9660 0.4265 0.7461 1.4046 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.5199 0.2692 -1.931 0.053481 . Farm: small 2.8713 0.7873 3.647 0.000265 *** Cattle: + 1.3953 0.5965 2.339 0.019324 * Farm size * Cattle -2.3253 1.0394 -2.237 0.025280 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 179.00 on 134 degrees of freedom Residual deviance: 147.59 on 131 degrees of freedom AIC: 155.59 Number of Fisher Scoring iterations: 4
Here, because of the presence of interaction, the main effects don't have a direct interpretation. Rather, we have to interpret them in the context of such phenomenon. At least in epidemiology, we work a lot using odds ratios by exponentiating variables' coefficients, and their 95% CI using the formula:
Upper CI = exp(coefficient + z * SE) (z being 1.96)
Lower CI = exp(coefficient - z * SE)
However, the calculation of odds ratios and their CI for the interaction term is not as straightforward. Following Hosmer and Lemeshow's book Applied Logistic Regression (2nd Ed., page 74, section 3.7), estimation of odds ratio in the presence of interaction should be derived from:
In our case, then, OR is calculated as follows:
Here, f represents Farm size, then f1 = 1 and f0 = 0. Also, x represents Cattle
OR <- round(exp(2.8713+(-2.3253)), digits = 2) print(paste("OR interaction =", OR)) "OR interaction = 1.73"
This OR represents the farm size effect (small - large) in farms with cattle, versus the farm size effect in farms without cattle. Or similarly, the odds of TB in small farms with cattle is 1.73 (73%) times the odds in small farms with no cattle.
The calculation of 95% CI for the interaction term is my main issue. Instead of using the term SE (1.0394, as calculated in the main model above), we need to use a variance-covariance matrix of the model and proceed as follows:
varcovar <- vcov(global) varcovar
Using this matrix, the 95% CI for the interaction term's OR is calculated as follows according to the above authors (page 76, same book) and others:
Replacing matix values in equation 3.15 above, we obtain the variance, from which the SE can be derived:
SE <- sqrt(0.619774148 + 1.080447013 * 1 + (2 * (-0.619774148))) SE  0.6787289
Lastly, we use the estimated SE to derive 95% CI for the interaction term:
OR_lower_CI <- exp(0.546 - 1.96 * SE) OR_lower_CI  0.4564352 OR_upper_CI <- exp(0.546 + 1.96 * SE) OR_upper_CI  6.529358
The interaction term of the model is highly significant (p = 0.025280), and the 95% CI for the OR should not include 1.
¿How can this be possible? One would expect the CI not to include 1 (e.g., OR: 1.73; 95% CI: 1.4-6.5).