HW6
Raavi Anvesh
November 11, 2018
7.5
library(vcdExtra)## Loading required package: vcd## Loading required package: grid## Loading required package: gnmdata("Caesar",package = "vcdExtra")
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c("Type 1","Type 2"))
Caesar.df$Infect## [1] 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 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 ...a
Caesar.df$Antibiotics=relevel(x = Caesar.df$Antibiotics,ref = "No")
Caesar.df$Planned=relevel(x = Caesar.df$Planned,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.logitmodel= 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.4b
summary(Caesar.logitmodel)##
## 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: 6anova(Caesar.logitmodel,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
## Risk 1 4.104 15 296.75 0.04278 *
## Antibiotics 1 55.163 14 241.59 1.11e-13 ***
## Planned 1 5.230 13 236.36 0.02220 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1c
(effectoninfection <- exp(coef(Caesar.logitmodel)) - 1)## (Intercept) RiskYes AntibioticsYes PlannedYes
## -0.5477367 5.2151576 -0.9502659 -0.5960215(effectoninfection2<- paste(round(100*effectoninfection, 2), "%", sep=""))## [1] "-54.77%" "521.52%" "-95.03%" "-59.6%"d
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.plot(Caesar.logitmodel)
Caesar.Separated=update(Caesar.logitmodel, . ~ . + Risk:Antibiotics)
plot(allEffects(Caesar.Separated))
```