Exercise 7.5 The data set Caesar in vcdExtra gives a 3 × 2 frequency table classifying 251 women who gave birth by Caesarian section by Infection (three levels: none, Type 1, Type2) and Risk, whether Antibiotics were used, and whether the Caesarian section was Planned or not. Infection is a natural response variable. In this exercise, consider only the binary outcome of infection vs. no infection.

library(vcdExtra)
## Loading required package: vcd
## Loading required package: grid
## Loading required package: gnm
data("Caesar",package="vcdExtra")
Caesar
## , , Antibiotics = Yes, Planned = Yes
## 
##          Risk
## Infection Yes No
##    Type 1   0  0
##    Type 2   1  1
##    None    17  1
## 
## , , Antibiotics = No, Planned = Yes
## 
##          Risk
## Infection Yes No
##    Type 1  11  4
##    Type 2  17  4
##    None    30 32
## 
## , , Antibiotics = Yes, Planned = No
## 
##          Risk
## Infection Yes No
##    Type 1   4  0
##    Type 2   7  0
##    None    87  0
## 
## , , Antibiotics = No, Planned = No
## 
##          Risk
## Infection Yes No
##    Type 1  10  0
##    Type 2  13  0
##    None     3  9
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c("Type 1", "Type 2"))

(a) Fit the main-effects logit model for the binary response Infect. Note that with the data in the form of a frequency data frame you will need to use weights=Freq in the call to glm(). (It might also be convenient to reorder the levels of the factors so that “No” is the baseline level for each.)

Caesar.df$Anti <- factor(Caesar.df$Antibiotics, levels(Caesar.df$Antibiotics)[c(2,1)])
Caesar.df$Ri <- factor(Caesar.df$Risk, levels(Caesar.df$Risk)[c(2,1)])
Caesar.df$Plan <- factor(Caesar.df$Planned, levels(Caesar.df$Planned)[c(2,1)])
Caesar.logistic <- glm(Infect ~ Anti + Ri + Plan, data = Caesar.df, family = binomial, weights=Freq)

(b) Use summary () or car::Anova () to test the terms in this model.

summary(Caesar.logistic)
## 
## Call:
## glm(formula = Infect ~ Anti + Ri + Plan, family = binomial, data = Caesar.df, 
##     weights = Freq)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.7471  -0.4426   0.0000   3.2338   5.4201  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7935     0.4785  -1.658   0.0972 .  
## AntiYes      -3.0011     0.4593  -6.535 6.37e-11 ***
## RiYes         1.8270     0.4364   4.186 2.84e-05 ***
## PlanYes      -0.9064     0.4084  -2.219   0.0265 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 300.85  on 16  degrees of freedom
## Residual deviance: 236.36  on 13  degrees of freedom
## AIC: 244.36
## 
## Number of Fisher Scoring iterations: 6
anova(Caesar.logistic, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Infect
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                    16     300.85              
## Anti  1   36.310        15     264.54 1.683e-09 ***
## Ri    1   22.956        14     241.59 1.657e-06 ***
## Plan  1    5.230        13     236.36    0.0222 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(c) Interpret the coefficients in the fitted model in terms of their effect on the odds of infection.

exp(coef(Caesar.logistic))
## (Intercept)     AntiYes       RiYes     PlanYes 
##  0.45226327  0.04973413  6.21515759  0.40397845

Interpretation: Holding everything else equal - if Antibiotics were used then the odds of infection decreases by 95%; if risk is yes then the odds of infection increases by 522% and if the Caesarian section was planned then the odds of infection decreases by 60%.

(d) Make one or more effects plots for this model, showing separate terms, or their combinations.

library(effects)
## Loading required package: carData
## 
## Attaching package: 'carData'
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
Caesar.eff <- allEffects(Caesar.logistic)
plot(Caesar.eff, rows=1, cols=3)

Caesar.logistic2 <- update(Caesar.logistic, . ~ . + Anti:Ri)
plot(allEffects(Caesar.logistic2), rows=1, cols=2)

Caesar.logistic3 <- update(Caesar.logistic, . ~ . + Anti:Plan)
plot(allEffects(Caesar.logistic3), rows=1, cols=2)

Caesar.logistic4 <- update(Caesar.logistic, . ~ . + Ri:Plan)
plot(allEffects(Caesar.logistic4), rows=1, cols=2)