# NOT RUN {
data(periphrasticDo)
# add midpoints of time periods
periphrasticDo$year = periphrasticDo$begin +
(periphrasticDo$end-periphrasticDo$begin)/2
# and ad an indicator variable distinguishing the first three time periods
# from the others
periphrasticDo$Indicator = rep(c(rep(0, 3), rep(1, 8)), 4)
# fit a logistic regression model
periphrasticDo.glm = glm(cbind(do, other) ~
(year + I(year^2) + I(year^3)) * type + Indicator * type +
Indicator * year, data = periphrasticDo, family = "binomial")
anova(periphrasticDo.glm, test = "F")
# visualization of data and model predictions
periphrasticDo$predict = predict(periphrasticDo.glm, type = "response")
par(mfrow=c(2, 2))
for (i in 1:nlevels(periphrasticDo$type)) {
subset = periphrasticDo[periphrasticDo$type ==
levels(periphrasticDo$type)[i], ]
plot(subset$year,
subset$do/(subset$do + subset$other),
type = "p", ylab = "proportion", xlab = "year",
ylim = c(0, 1), xlim = c(1400, 1700))
mtext(levels(periphrasticDo$type)[i], line = 2)
lines(subset$year, subset$predict, lty = 1)
}
par(mfrow=c(1, 1))
# }
Run the code above in your browser using DataLab