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

Convert data into binary outcome of infection vs. no infection.

Caesar.df=as.data.frame(Caesar)
Caesar.df$Infect =as.numeric(Caesar.df$Infection %in% c("Type 1", "Type 2"))
Caesar.df
##    Infection Risk Antibiotics Planned Freq Infect
## 1     Type 1  Yes         Yes     Yes    0      1
## 2     Type 2  Yes         Yes     Yes    1      1
## 3       None  Yes         Yes     Yes   17      0
## 4     Type 1   No         Yes     Yes    0      1
## 5     Type 2   No         Yes     Yes    1      1
## 6       None   No         Yes     Yes    1      0
## 7     Type 1  Yes          No     Yes   11      1
## 8     Type 2  Yes          No     Yes   17      1
## 9       None  Yes          No     Yes   30      0
## 10    Type 1   No          No     Yes    4      1
## 11    Type 2   No          No     Yes    4      1
## 12      None   No          No     Yes   32      0
## 13    Type 1  Yes         Yes      No    4      1
## 14    Type 2  Yes         Yes      No    7      1
## 15      None  Yes         Yes      No   87      0
## 16    Type 1   No         Yes      No    0      1
## 17    Type 2   No         Yes      No    0      1
## 18      None   No         Yes      No    0      0
## 19    Type 1  Yes          No      No   10      1
## 20    Type 2  Yes          No      No   13      1
## 21      None  Yes          No      No    3      0
## 22    Type 1   No          No      No    0      1
## 23    Type 2   No          No      No    0      1
## 24      None   No          No      No    9      0
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 ...
  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.) ##Change the order of of factors so that “No” is the first level.
Caesar.df$Planned=relevel(x = Caesar.df$Planned,ref = "No")
Caesar.df$Antibiotics=relevel(x = Caesar.df$Antibiotics,ref = "No")
Caesar.df$Risk=relevel(x = Caesar.df$Risk,ref = "No")
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.model= glm(Infect ~ Risk + Antibiotics + Planned, weights = Freq, data= Caesar.df, family = "binomial" ))
## 
## Call:  glm(formula = Infect ~ Risk + Antibiotics + Planned, family = "binomial", 
##     data = Caesar.df, weights = Freq)
## 
## Coefficients:
##    (Intercept)         RiskYes  AntibioticsYes      PlannedYes  
##        -0.7935          1.8270         -3.0011         -0.9064  
## 
## Degrees of Freedom: 16 Total (i.e. Null);  13 Residual
## Null Deviance:       300.9 
## Residual Deviance: 236.4     AIC: 244.4
  1. Use summary () or car::Anova () to test the terms in this model
anova(Caesar.model)
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Infect
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev
## NULL                           16     300.85
## Risk         1    4.104        15     296.75
## Antibiotics  1   55.163        14     241.59
## Planned      1    5.230        13     236.36
  1. Interpret the coefficients in the fitted model in terms of their effect on the odds of infection.
(effect <- exp(coef(Caesar.model)) - 1)
##    (Intercept)        RiskYes AntibioticsYes     PlannedYes 
##     -0.5477367      5.2151576     -0.9502659     -0.5960215
(effect2<- paste(round(100*effect, 2), "%", sep=""))
## [1] "-54.77%" "521.52%" "-95.03%" "-59.6%"

The coefficients of the model indicate that with the odds of infection rish increases 521%, whereas the presence of antibodies decreases the likelihod by 95% and chance of infection by 59%.

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

#install.packages("carData")
#library(carData)
#install.packages("effects")
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.Separated=update(Caesar.model, . ~ . + Risk:Antibiotics)
plot(allEffects(Caesar.Separated))