#----------------------------------------------#

# R in Action (2nd ed): Chapter 13             #

# Generalized linear models                    #

# requires packages AER, robust, gcc           #

# install.packages(c("AER", "robust", "gcc"))  #

#----------------------------------------------#

## Logistic Regression

# get summary statistics

data(Affairs, package="AER")

summary(Affairs)

table(Affairs$affairs)

# create binary outcome variable

Affairs$ynaffair[Affairs$affairs > 0] <- 1

Affairs$ynaffair[Affairs$affairs == 0] <- 0

Affairs$ynaffair <- factor(Affairs$ynaffair,

                           levels=c(0,1),

                           labels=c("No","Yes"))

table(Affairs$ynaffair)

# fit full model

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +

                  religiousness + education + occupation +rating,

                data=Affairs,family=binomial())

summary(fit.full)

# fit reduced model

fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness +

                     rating, data=Affairs, family=binomial())

summary(fit.reduced)

# compare models

anova(fit.reduced, fit.full, test="Chisq")

# interpret coefficients

coef(fit.reduced)

exp(coef(fit.reduced))

# calculate probability of extramariatal affair by marital ratings

testdata <- data.frame(rating = c(1, 2, 3, 4, 5),

                       age = mean(Affairs$age),

                       yearsmarried = mean(Affairs$yearsmarried),

                       religiousness = mean(Affairs$religiousness))

testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")

testdata

# calculate probabilites of extramariatal affair by age

testdata <- data.frame(rating = mean(Affairs$rating),

                       age = seq(17, 57, 10),

                       yearsmarried = mean(Affairs$yearsmarried),

                       religiousness = mean(Affairs$religiousness))

testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")

testdata

# evaluate overdispersion

fit <- glm(ynaffair ~ age + yearsmarried + religiousness +

             rating, family = binomial(), data = Affairs)

fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +

                rating, family = quasibinomial(), data = Affairs)

pchisq(summary(fit.od)$dispersion * fit$df.residual,

       fit$df.residual, lower = F)

## Poisson Regression

# look at dataset

data(breslow.dat, package="robust")

names(breslow.dat)

summary(breslow.dat[c(6, 7, 8, 10)])

# plot distribution of post-treatment seizure counts

opar <- par(no.readonly=TRUE)

par(mfrow=c(1, 2))

attach(breslow.dat)

hist(sumY, breaks=20, xlab="Seizure Count",

     main="Distribution of Seizures")

boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons")

par(opar)

# fit regression

fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())

summary(fit)

# interpret model parameters

coef(fit)

exp(coef(fit))

# evaluate overdispersion

deviance(fit)/df.residual(fit)

library(qcc)

qcc.overdispersion.test(breslow.dat$sumY, type="poisson")

# fit model with quasipoisson

fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,

              family=quasipoisson())

summary(fit.od)