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(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(vcdExtra)
## Loading required package: vcd
## Loading required package: grid
## Loading required package: gnm
## 
## Attaching package: 'vcdExtra'
## The following object is masked from 'package:carData':
## 
##     Burt
data("Caesar", package = "vcdExtra")
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c('Type 1', 'Type 2'))
str(Caesar.df)
## 'data.frame':    24 obs. of  6 variables:
##  $ Infection  : Factor w/ 3 levels "Type 1","Type 2",..: 1 2 3 1 2 3 1 2 3 1 ...
##  $ Risk       : Factor w/ 2 levels "Yes","No": 1 1 1 2 2 2 1 1 1 2 ...
##  $ Antibiotics: Factor w/ 2 levels "Yes","No": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Planned    : Factor w/ 2 levels "Yes","No": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Freq       : num  0 1 17 0 1 1 11 17 30 4 ...
##  $ Infect     : num  1 1 0 1 1 0 1 1 0 1 ...
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)])
str(Caesar.df)
## 'data.frame':    24 obs. of  6 variables:
##  $ Infection  : Factor w/ 3 levels "Type 1","Type 2",..: 1 2 3 1 2 3 1 2 3 1 ...
##  $ Risk       : Factor w/ 2 levels "No","Yes": 2 2 2 1 1 1 2 2 2 1 ...
##  $ Antibiotics: Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 1 1 1 1 ...
##  $ Planned    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Freq       : num  0 1 17 0 1 1 11 17 30 4 ...
##  $ Infect     : num  1 1 0 1 1 0 1 1 0 1 ...
Caesar.logistic <- glm(Infect ~ Risk + Antibiotics + Planned, data = Caesar.df, family = binomial, weights=Freq)
summary(Caesar.logistic)
## 
## 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 2.84e-05 ***
## AntibioticsYes  -3.0011     0.4593  -6.535 6.37e-11 ***
## 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
effect <- exp(coef(Caesar.logistic)) - 1
effect.perc <- paste(round(100*effect, 2), "%", sep="")
effect.perc.abs <- paste(round(100*abs(effect), 2), "%", sep="")
effect.perc
## [1] "-54.77%" "521.52%" "-95.03%" "-59.6%"

when an antibiotic is present, the risk of infection decrease.the reduction goes in -95%

effect <- exp(coef(Caesar.logistic)) - 1
effect.perc <- paste(round(100*effect, 2), "%", sep="")
effect.perc.abs <- paste(round(100*abs(effect), 2), "%", sep="")
effect.perc
## [1] "-54.77%" "521.52%" "-95.03%" "-59.6%"
plot(Caesar.logistic)

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