HW6

library("vcdExtra")

## Warning: package 'vcdExtra' was built under R version 3.4.4

## Loading required package: vcd

## Warning: package 'vcd' was built under R version 3.4.4

## Loading required package: grid

## Loading required package: gnm

## Warning: package 'gnm' was built under R version 3.4.4

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$Antibiotics2 <- factor(Caesar.df$Antibiotics, levels(Caesar.df$Antibiotics)[c(2,1)])
Caesar.df$Risk2 <- factor(Caesar.df$Risk, levels(Caesar.df$Risk)[c(2,1)])
Caesar.df$Planned2 <- factor(Caesar.df$Planned, levels(Caesar.df$Planned)[c(2,1)])
LogisticModel <- glm(Infect ~ Antibiotics2 + Risk2 + Planned2, data = Caesar.df, family = binomial, weights=Freq)

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

summary(LogisticModel)

## 
## Call:
## glm(formula = Infect ~ Antibiotics2 + Risk2 + Planned2, 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 .  
## Antibiotics2Yes  -3.0011     0.4593  -6.535 6.37e-11 ***
## Risk2Yes          1.8270     0.4364   4.186 2.84e-05 ***
## Planned2Yes      -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

3. Interpret the coefficients in the fitted model in terms of their effect on the odds of infection.

exp(coef(LogisticModel))

##     (Intercept) Antibiotics2Yes        Risk2Yes     Planned2Yes 
##      0.45226327      0.04973413      6.21515759      0.40397845

effect <- exp(coef(LogisticModel)) - 1
effect

##     (Intercept) Antibiotics2Yes        Risk2Yes     Planned2Yes 
##      -0.5477367      -0.9502659       5.2151576      -0.5960215

Model Coefficient interpretation: If all else holds equal, when antibiotics were used, the odds of getting infection decreases by 95%; when risk is involved (? no explaination on this variable), the odds for getting infection increases by 522%; When C-section is planned, the odds of infection decreases by 60%.

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

library(effects)

## Warning: package 'effects' was built under R version 3.4.3

## 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.

ModelEffects <- allEffects(LogisticModel)
plot(ModelEffects, rows=1, cols=3)

LogisticModel2<-update(LogisticModel, .~. + Antibiotics2:Planned2)
anova(LogisticModel2, 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              
## Antibiotics2           1   36.310        15     264.54 1.683e-09 ***
## Risk2                  1   22.956        14     241.59 1.657e-06 ***
## Planned2               1    5.230        13     236.36    0.0222 *  
## Antibiotics2:Planned2  1    1.179        12     235.18    0.2776    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(LogisticModel2), rows=1, cols=3)