Second part on logistic regression (first one here).

We used in the previous post a likelihood ratio test to compare a full and null model. The same can be done to compare a full and nested model to test the contribution of any subset of parameters:

mod1.nest <- glm(casecont ~ dcpct, family = binomial("logit"), + data = nocardia) lr.mod1.nest <- -(deviance(mod1.nest) / 2) (lr <- 2 * (lr.mod1 - lr.mod1.nest)) [1] 30.16203 1 - pchisq(lr, 2) [1] 2.820974e-07 ### or, more straightforward, using anova to compare the 2 models: anova(mod1, mod1.nest, test = "Chisq") Analysis of Deviance Table Model 1: casecont ~ dcpct + dneo + dclox Model 2: casecont ~ dcpct Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 104 107.99 2 106 138.15 -2 -30.162 2.821e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

**Interpretation of coefficients**

### example 16.2 nocardia$dbarn <- as.factor(nocardia$dbarn) mod2 <- glm(casecont ~ dcpct + dneo + dclox + dbarn, + family = binomial("logit"), data = nocardia) (mod2.sum <- summary(mod2)) Call: glm(formula = casecont ~ dcpct + dneo + dclox + dbarn, family = binomial("logit"), data = nocardia) Deviance Residuals: Min 1Q Median 3Q Max -1.6949 -0.7853 0.1021 0.7692 2.6801 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.445696 0.854328 -2.863 0.00420 ** dcpct 0.021604 0.007657 2.821 0.00478 ** dneoYes 2.685280 0.677273 3.965 7.34e-05 *** dcloxYes -1.235266 0.580976 -2.126 0.03349 * dbarntiestall -1.333732 0.631925 -2.111 0.03481 * dbarnother -0.218350 1.154293 -0.189 0.84996 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 149.72 on 107 degrees of freedom Residual deviance: 102.32 on 102 degrees of freedom AIC: 114.32 Number of Fisher Scoring iterations: 5 cbind(exp(coef(mod2)), exp(confint(mod2))) Waiting for profiling to be done... 2.5 % 97.5 % (Intercept) 0.08666577 0.01410982 0.4105607 dcpct 1.02183934 1.00731552 1.0383941 dneoYes 14.66230075 4.33039752 64.5869271 dcloxYes 0.29075729 0.08934565 0.8889877 dbarntiestall 0.26349219 0.06729031 0.8468235 dbarnother 0.80384385 0.08168466 8.2851189

Note: Dohoo do not report the profile likelihood-based confidence interval on page 404. As said in the previous post, the profile likelihood-based CI is preferable due to the Hauck-Donner effect (overestimation of the SE).

Using neomycin in the herd increases the log odds of *Nocardia* mastitis by 2.7 units (or using neomycin increases the odds 14.7 times). Since *Nocardia* mastitis is a rare condition, we can interpret the odds ratio as a risk ratio (RR) and say neomycin increases the risk of *Nocardia* mastitis by 15 times. Changing the percentage of dry cows treated from 50 to 75% increases the log odds of disease by units (or ). An increase of 25% in the percentage of cows dry-treated increases the risk of disease by about 72% (i.e. 1.72 times). Tiestall barns and other barn types both have lower risks of *Nocardia* mastitis than freestall barns.

The significance of the main effects can be tested with:

drop1(mod2, test = "Chisq") Single term deletions Model: casecont ~ dcpct + dneo + dclox + as.factor(dbarn) Df Deviance AIC LRT Pr(>Chi) <none> 102.32 114.32 dcpct 1 111.36 121.36 9.0388 0.002643 ** dneo 1 123.91 133.91 21.5939 3.369e-06 *** dclox 1 107.02 117.02 4.7021 0.030126 * as.factor(dbarn) 2 107.99 115.99 5.6707 0.058697 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The drop1 function tests each predictor relative to the full model.

**Presenting effects of factors on the probability scale**

Usually we think about the probability of disease rather than the odds, and the probability of disease is not linearly related to the factor of interest. A unit increase in the factor does not increase the probability of disease by a fixed amount. It depends on the level of the factor and the levels of other factors in the model.

mod1 <- glm(casecont ~ dcpct + dneo + dclox, family = binomial("logit"), + data = nocardia) nocardia$neo.pred <- predict(mod1, type = "response", se.fit = FALSE) library(ggplot2) ggplot(nocardia, aes(x = dcpct, y = neo.pred, group = dneo, + colour = factor(dneo))) + + stat_smooth(method = "glm", family = "binomial", se = FALSE) + + labs(colour = "Neomycin", x = "Percent of cows dry treated", + y = "Probability of Nocardia")

Next: evaluating logistic regression models.