ANLY 545- Homework #6
Juan Colunga
July 26, 2018
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.
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.4data("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 9Converting 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 0str(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 ...- 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- 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- 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 in the fitted model, in the terms of their effect on the odds of infection, indicate that 1) The existence of risk increases by 521% when infection, 2) The presence of antibiotics decreases the likelihood by -95% and decreases the probability of infection by 60%.
- Make one or more effects plots for this model, showing separate terms, or their combinations.
plot(Caesar.model)
library(effects)## Warning: package 'effects' was built under R version 3.4.4## Loading required package: carData## Warning: package 'carData' was built under R version 3.4.4##
## 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))