##---------------------------------------------------------------------## ## Script for Generalized Linear Models ## ## John Fox ## ## York SPIDA 2010 ## ##---------------------------------------------------------------------## # Poisson regression for Ornstein's interlocking-directorate data library(car) ?Ornstein attach(Ornstein) (tab <- table(interlocks)) (x <- sort(unique(interlocks))) plot(x, tab, type='h', xlab='Number of Interlocks', ylab='Frequency') points(x, tab, pch=16) detach(Ornstein) mod.ornstein <- glm(interlocks ~ assets + nation + sector, family=poisson, data=Ornstein) summary(mod.ornstein) Anova(mod.ornstein) # quasi-Poisson model mod.ornstein.qp <- glm(interlocks ~ assets + nation + sector, family=quasipoisson, data=Ornstein) summary(mod.ornstein.qp) Anova(mod.ornstein.qp) Anova(mod.ornstein.qp, test="F") # effect plot, based on quasi-Poisson model library(effects) plot(allEffects(mod.ornstein.qp), ask=FALSE) # diagnostics # unusual data influencePlot(mod.ornstein.qp) outlier.test(mod.ornstein.qp) D <- dfbeta(mod.ornstein.qp) head(D) plot(D[,"assets"]) # index plot abline(h=c(-3.389e-06, 0, 3.389e-06), lty=2) # +/- SE(B) identify(D[,"assets"]) # nonlinearity cr.plot(mod.ornstein.qp, "assets", span=.9) mod.ornstein.qp.1 <- glm(interlocks ~ log10(assets) + nation + sector, family=quasipoisson, data=Ornstein) summary(mod.ornstein.qp.1) deviance(mod.ornstein.qp) cr.plot(mod.ornstein.qp.1, "log10(assets)") influencePlot(mod.ornstein.qp.1) # binomial logit model for American Voter data closeness <- factor(rep(c('one.sided', 'close'), c(3, 3)), levels=c('one.sided', 'close')) preference <- factor(rep(c('weak', 'medium', 'strong'), 2), levels=c('weak', 'medium', 'strong')) voted <- c(91, 121, 64, 214, 284, 201) did.not.vote <- c(39, 49, 24, 87, 76, 25) AmVoter <- data.frame(did.not.vote, voted, closeness, preference) AmVoter options("contrasts") options(contrasts=c("contr.sum", "contr.poly")) # "deviation" regressors amvoter <- glm(cbind(voted, did.not.vote) ~ closeness*preference, family=binomial, data=AmVoter) # "saturated" model summary(amvoter) amvoter0 <- update(amvoter, . ~ . - closeness:preference) anova(amvoter0, amvoter, test="Chisq") # LR test of interaction Anova(amvoter) Anova(amvoter, type=3) # loglinear model for American Voter data closeness <- rep(closeness, 2) preference <- rep(preference, 2) count <- c(voted, did.not.vote) turnout <- factor(rep(c("voted", "did.not.vote"), each=6)) AmVoter <- data.frame(closeness, preference, turnout, count) AmVoter amvoter <- glm(count ~ closeness*preference*turnout, family=poisson, data=AmVoter) # saturated model summary(amvoter) Anova(amvoter)