library(visreg)
## Warning: package 'visreg' was built under R version 3.5.1
library(vcd)
## Warning: package 'vcd' was built under R version 3.5.1
## Loading required package: grid
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 3.5.1
## Loading required package: gnm
library(effects)
## Warning: package 'effects' was built under R version 3.5.1
## 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.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.1
library(car)
## Warning: package 'car' was built under R version 3.5.1
library(ca)
## Warning: package 'ca' was built under R version 3.5.1
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.1
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following object is masked from 'package:vcdExtra':
## 
##     summarise
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lmtest)
## Warning: package 'lmtest' was built under R version 3.5.1
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(ISLR)
## Warning: package 'ISLR' was built under R version 3.5.1
## 
## Attaching package: 'ISLR'
## The following object is masked from 'package:vcd':
## 
##     Hitters
library(MASS)
## Warning: package 'MASS' was built under R version 3.5.1
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(ROCR)
## Warning: package 'ROCR' was built under R version 3.5.1
## Loading required package: gplots
## Warning: package 'gplots' was built under R version 3.5.1
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess

Exercise 7.5. The Caesar data set in vcdExtra:

data("Caesar", package = "vcdExtra")
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in%
                                 c("Type 1", "Type 2"))

To get a feel for this data set, let’s look at its structure:

str(Caesar)
##  'table' num [1:3, 1:2, 1:2, 1:2] 0 1 17 0 1 1 11 17 30 4 ...
##  - attr(*, "dimnames")=List of 4
##   ..$ Infection  : chr [1:3] "Type 1" "Type 2" "None"
##   ..$ Risk       : chr [1:2] "Yes" "No"
##   ..$ Antibiotics: chr [1:2] "Yes" "No"
##   ..$ Planned    : chr [1:2] "Yes" "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 "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 ...
str(Caesar.df$Infect)
##  num [1:24] 1 1 0 1 1 0 1 1 0 1 ...

To add a new level of “0” for “none,” we use the following:

Caesar.df$Infection <- factor(Caesar.df$Infection, levels = c("None", "Type 1", "Type 2"))

We can reorder the variables to have level “No” as “zero” so that “No” will not be assigned to “One.”

Caesar.df$Risk <- factor(Caesar.df$Risk, levels = c("No", "Yes"))

Let’s look at the relation between Infection and the other factors:

Infection vs no infection:

Caesar.df.Inf <- glm(Infection ~ Risk + Antibiotics + Planned, family = binomial, weights = Freq, data = Caesar)

Caesar.df.Inf
## 
## Call:  glm(formula = Infection ~ Risk + Antibiotics + Planned, family = binomial, 
##     data = Caesar, weights = Freq)
## 
## Coefficients:
##   (Intercept)         RiskNo  AntibioticsNo      PlannedNo  
##        3.9890         1.2414        -2.6335        -0.7402  
## 
## Degrees of Freedom: 16 Total (i.e. Null);  13 Residual
## Null Deviance:       179.7 
## Residual Deviance: 154.5     AIC: 162.5

(a) Fit the main-effects logit model for the binary response, Infect

Caesar.ef <-allEffects(Caesar.df.Inf)
Caesar.ef
##  model: Infection ~ Risk + Antibiotics + Planned
## 
##  Risk effect
## Risk
##       Yes        No 
## 0.9090564 0.9719023 
## 
##  Antibiotics effect
## Antibiotics
##       Yes        No 
## 0.9857920 0.8328685 
## 
##  Planned effect
## Planned
##       Yes        No 
## 0.9641865 0.9277597
plot(Caesar.ef)

Results: It seems there was no risk of infection in the Caesarian section–antibiotics were taken. It was difficult to determine whether the Caesarian section was planned.

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

summary(Caesar.ef, test = "Chisq")
##  model: Infection ~ Risk + Antibiotics + Planned
## 
##  Risk effect
## Risk
##       Yes        No 
## 0.9090564 0.9719023 
## 
##  Lower 95 Percent Confidence Limits
## Risk
##       Yes        No 
## 0.8490996 0.9149531 
## 
##  Upper 95 Percent Confidence Limits
## Risk
##       Yes        No 
## 0.9466865 0.9910885 
## 
##  Antibiotics effect
## Antibiotics
##       Yes        No 
## 0.9857920 0.8328685 
## 
##  Lower 95 Percent Confidence Limits
## Antibiotics
##       Yes        No 
## 0.9545677 0.7341776 
## 
##  Upper 95 Percent Confidence Limits
## Antibiotics
##       Yes        No 
## 0.9956544 0.8999138 
## 
##  Planned effect
## Planned
##       Yes        No 
## 0.9641865 0.9277597 
## 
##  Lower 95 Percent Confidence Limits
## Planned
##       Yes        No 
## 0.9164305 0.8521633 
## 
##  Upper 95 Percent Confidence Limits
## Planned
##       Yes        No 
## 0.9850960 0.9662316

Results: Based on the summary above, there seems to be no risk involved with infection by Caesarian section.

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

coef(Caesar.df.Inf)
##   (Intercept)        RiskNo AntibioticsNo     PlannedNo 
##     3.9890313     1.2413990    -2.6335441    -0.7401842

exp(coef(Caesar.df.Inf))
##   (Intercept)        RiskNo AntibioticsNo     PlannedNo 
##   54.00255375    3.46045136    0.07182346    0.47702603

exp(10*coef(Caesar.df.Inf))
##   (Intercept)        RiskNo AntibioticsNo     PlannedNo 
##  2.109322e+17  2.462224e+05  3.653116e-12  6.101278e-04

Results: When the individual variable is increased by three, the exponent of improvement is improved by 2.46 %. (I am not so sure here.)

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

plot(Caesar.ef)

Results: Again, it seems there was no risk of infection in the Caesarian section–antibiotics were taken. It was difficult to determine whether the Caesarian section was planned.

plot(Caesar)

Results: In this plot, there is strong evidence of no risk of infection.