library(tidyverse)
library(magrittr)
library(vcd)
library(vcdExtra)
library(effects)

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.

data("Caesar", package = "vcdExtra")
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c("Type 1", "Type 2"))
  1. 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 <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c("Type 1", "Type 2"))
Caesar.df$Risk <- factor(Caesar.df$Risk, levels(Caesar.df$Risk)[c(2,1)])
Caesar.df$Antibiotics <- factor(Caesar.df$Antibiotics, levels(Caesar.df$Antibiotics)[c(2,1)])
Caesar.df$Planned <- factor(Caesar.df$Planned, levels(Caesar.df$Planned)[c(2,1)])
mod1 <- glm(Infect ~ Risk + Antibiotics + Planned, family = binomial, data = Caesar.df, weights = Freq)
  1. Use summary () or car::Anova () to test the terms in this model.
summary(mod1)
## 
## Call:
## glm(formula = Infect ~ Risk + Antibiotics + Planned, 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 .  
## RiskYes          1.8270     0.4364   4.186 0.0000283815346 ***
## AntibioticsYes  -3.0011     0.4593  -6.535 0.0000000000637 ***
## PlannedYes      -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(mod1)
  1. Interpret the coefficients in the fitted model in terms of their effect on the odds of infection.
ef <- exp(coef(mod1)) - 1
ef.perc <- paste(round(100 * ef, 2), "%", sep = "")
ef.perc.abs <- paste(round(100 * abs(ef), 2), "%", sep = "")
ef.perc
## [1] "-54.77%" "521.52%" "-95.03%" "-59.6%"

From the result we can conclude that the odds of infection decreases by 95% when an antibiotic is present, and decrease by 60% when it is planned.

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

mod1.inter <- update(mod1, . ~ . + Risk:Antibiotics)
plot(allEffects(mod1.inter))