library(vcdExtra)
library(effects)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.
data("Caesar")
Caesar.df <- as.data.frame(Caesar)
Caesar.df$Infect <- as.numeric(Caesar.df$Infection %in% c("Type 1", "Type 2"))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$Antibiotics <- factor(Caesar.df$Antibiotics, levels(Caesar.df$Antibiotics)[c(2,1)])
Caesar.df$Risk <- factor(Caesar.df$Risk, levels(Caesar.df$Risk)[c(2,1)])
Caesar.df$Planned <- factor(Caesar.df$Planned, levels(Caesar.df$Planned)[c(2,1)])
Caesar.logistic <- glm(Infect ~ Antibiotics + Risk + Planned, data = Caesar.df, family = binomial, weights=Freq)
Caesar.logistic##
## Call: glm(formula = Infect ~ Antibiotics + Risk + Planned, family = binomial,
## data = Caesar.df, weights = Freq)
##
## Coefficients:
## (Intercept) AntibioticsYes RiskYes PlannedYes
## -0.7935 -3.0011 1.8270 -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.
summary(Caesar.logistic)##
## Call:
## glm(formula = Infect ~ Antibiotics + Risk + 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 .
## AntibioticsYes -3.0011 0.4593 -6.535 6.37e-11 ***
## RiskYes 1.8270 0.4364 4.186 2.84e-05 ***
## 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: 6
anova(Caesar.logistic, test = "Chisq")Interpret the coefficients in the fitted model in terms of their effect on the odds of infection.
odds <- exp(coef(Caesar.logistic))-1
odds.percent <- paste(round(100*odds, 2), "%", sep="")
names(odds.percent) <- names(odds)
odds.percent## (Intercept) AntibioticsYes RiskYes PlannedYes
## "-54.77%" "-95.03%" "521.52%" "-59.6%"
Make one or more effects plots for this model, showing separate terms, or their combinations.
plot(Caesar.logistic)plot(allEffects(Caesar.logistic), rows = 1, cols = 3)plot(allEffects(update(Caesar.logistic, . ~ . + Antibiotics:Risk)))plot(allEffects(update(Caesar.logistic, . ~ . + Antibiotics:Planned)))plot(allEffects(update(Caesar.logistic, . ~ . + Risk:Planned)))