---
title: "Binary and Count data"
author: "Nancy"
date: "2/11/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(SMPracticals)
library(faraway)
library(tidyverse)
library(dplyr)
options(scipen=10, digits=4) # gets rid of exponential notation
```

## Binary data

```{r data}
data(wcgs, package = "faraway")
head(wcgs) #not run: str(wcgs); plot(wcgs); help(wcgs)
sum(is.na(wcgs)) # there are a few missing values; much further sleuthing finds arcus and chol have some missing values
```
Fitting a binary model:
```{r modelling}
glm(chd ~ ., family = binomial, data = wcgs)
## typechd and timechd confounded with chd; behave and dibepB are the same variable recorded differently
heartmod <- glm(chd ~ age + height + weight + sdp + dbp + chol + behave + cigs, family = binomial, data = wcgs)
summary(heartmod)
heartmod2 <- update(heartmod, . ~ . - behave)
summary(heartmod2)
anova(heartmod2, heartmod, test = "LRT") #default for test is NULL
```
There is a significant effect of behaviour; a check using `levels(wcgs$behave)` shows that it has levels A1, A2, B1, B2. The variable `dibep` collapses these to A and B. 
```{r typea}
heartmod3 <- update(heartmod, .~ . - behave + dibep )
summary(heartmod3)
```
Quick check -- how do we interpret the coefficent estimate for `dibepB`?  

## Forensics 

The coefficient above for `dibepB` is positive, and significantly so, suggesting an *increase* in risk of `chd` for "type B" relative to "type A". This is not consistent with the coefficients for `behave` above, where both B1 and B2 showed a decrease risk relative to A1. What is going on? 
```{r crosstabs}
xtabs(~ chd + behave, data = wcgs)
xtabs(~ chd + dibep, data = wcgs)
1155 + 331; 1177+234
```
It looks like `dibepB` is actually `dibepA`, if `dibepA` is meant to be the aggregation of A1 and A2, which is what is says in the help file. The original paper is available [on my web page](https://utstat.toronto.edu/reid/sta2101f/jama_233_8_016.html). In their Table 3 ("Intake all ages") we see 178 CHD in the Type A group and 79 CHD in the Type B group, as above once we switch the labels. 

## Example 10.18 from SM
Now I'll have a look at the nodal involvment data from SM. The data as presented in the book looks like
![ image](/Users/nancy/teaching/2201/2011-2/Table108.jpg)
But the dataset is recorded as binary:
```{r nodal}
data(nodal)
head(nodal); dim(nodal) # sum(nodal$m) = 53 not shown
```
We can group them using the `tidyverse`:
```{r nodal2}
nodal2 <- nodal %>% 
  group_by(aged,stage,grade,xray,acid) %>% 
  summarise(total=sum(m), rtot = sum(r)) #Table 10.8
head(nodal2)
```
First I'll fit with the binary data
```{r nodal3 }
Ex1018bin.glm <- glm(r ~ aged+stage+grade+xray+acid, data = nodal, family = binomial)
plot(Ex1018bin.glm$linear.predictors, residuals(Ex1018bin.glm, type="deviance"), main = "Binary Fit")
plot(Ex1018bin.glm$fitted.values, residuals(Ex1018bin.glm, type="deviance"), main = "Binary Fit")
summary(Ex1018bin.glm)
```
And now with the binomial version (grouped in covariate classes as above)
```{r nodal4}
Ex1018binom.glm <- glm(cbind(rtot,total-rtot) ~ 
                    aged + stage + grade + xray + acid, data = nodal2, 
                  family = binomial)

summary(Ex1018binom.glm)

attributes(Ex1018binom.glm)

plot(Ex1018binom.glm$linear.predictors,residuals(Ex1018binom.glm, type="deviance"), main = "Binomial Fit")

plot(Ex1018binom.glm$fitted.values,residuals(Ex1018binom.glm, type="deviance"), main = "Binomial Fit")

qqnorm(residuals(Ex1018binom.glm, type="deviance"))
qqnorm(residuals(Ex1018binom.glm, type="pearson"))

plot(Ex1018binom.glm, main = "Binomial Fit, R default")

plot.glm.diag(Ex1018binom.glm, main = "Binomial Fit, Davison's SM")

?plot.glm.diag 
```
From the help file: "Makes plot of jackknife deviance residuals against linear predictor, normal scores plots of standardized deviance residuals, plot of approximate Cook statistics against leverage/(1-leverage), and case plot of Cook statistic."  

## Poisson regression
This example in Chapter 3 of FELM is used to introduce Poisson regression. Faraway first considers ordinary linear regression, and notes that the residuals don't seem to have constant variance. Back in the day we'd take logs or square root of the response and then fit linear regression again, but a more modern approach is to fit a Poisson glm.

The residual deviance in a Poisson model is a test of fit of the model, as it is with a Binomial.  The second fit below allows an inflation factor on the variance, while still using the simplicity of the Poisson for fitting. 
```{r poisson gala}
data(gala)
head(gala)
help(gala)

gal.glm <- glm(Species ~ . - Endemics,
             family = poisson,
             data = gala)

summary(gal.glm)

gal.quasi <- glm(Species ~ . - Endemics,
                 family = quasipoisson,
                 data = gala)
summary(gal.quasi)
plot(gal.quasi$linear.predictors, gal.quasi$residuals)
```

The canonical link for the Poisson is the log link. Identity links are sometimes used as well; see the `cloth` data analysis in SM Example 10.28.
