---
title: "Binary regression ELM-2 Ch.2"
author: ""
date: "18/10/2022"
output: html_document
---

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

## Binary data

```{r data}
data(wcgs, package = "faraway")
dim(wcgs)
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: The first fit is silly, luckily `R` gives a warning.
```{r modelling}
glm(chd ~ ., family = binomial, data = wcgs)
```

It turns out that `typechd` and `timechd` are confounded with `chd`; and that `behave` and `dibepB` are the same variable recorded differently. So we'll try:
```{r modelling again}
heartmod <- glm(chd ~ age + height + weight + sdp + dbp + chol + behave + cigs, family = binomial, data = wcgs)
sumary(heartmod)
heartmod2 <- update(heartmod, . ~ . - behave)
sumary(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 )
sumary(heartmod3)
```
In Sections 2.1-2.4 of ELM-2, just two covariates (i.e. potential explanatory variables) are used: `height` and `cigs`.
```{r follow ELM-2}
lmod <- glm(chd ~ height + cigs, family = binomial, data = wcgs)
beta <- coef(lmod)
plot(jitter(as.numeric(chd)-1, amount = 0.1) ~ jitter(cigs), wcgs, xlab="cigs", ylab = "Heart Disease", pch=".")
curve(ilogit(beta[1]+beta[2]*60 + beta[3]*x), add=TRUE)
curve(ilogit(beta[1]+beta[2]*78 + beta[3]*x), add=TRUE, lty=2, col="blue")
```

Inference for $\beta$ follows as in the challenger shuttle example; the only difference is that our response is binary instead of binomial. This makes the raw residuals $y - \hat y$ look weird, because we're always subtracting from either $1$ or $0$. In Section 2.4, Faraway groups the residuals and plots the average within each group
```{r residuals}
plot.glm.diag(heartmod3)
par(mfrow=c(2,2))
plot(heartmod3)
```
For *model selection*, $AIC = 2 \text{maximum log-likelihood} + 2p$  is very popular. The maximized value of the log-likelihood function plays the role of the *minimized* residual sum of squares.
```{r model selection}
heartmod4 <- MASS::stepAIC(heartmod3)
sumary(heartmod4)
```
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 it says in the help file. The original paper is available [on my web page](https://utstat.toronto.edu/reid/sta2101f/2021/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`: (not shown: `View(nodal2)`)
```{r nodal2}
library(tidyverse)
nodal2 <- nodal %>% 
  group_by(aged,stage,grade,xray,acid) %>% 
  summarise(total=sum(m), rtot = sum(r)) #Table 10.8
#View(nodal2)
```
First I'll fit with the binary data
```{r nodal3 }
par(mfrow=c(2,2))
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)
par(mfrow=c(2,2))
plot(Ex1018binom.glm$linear.predictors,residuals(Ex1018binom.glm, type="deviance"), ylab = "deviance residuals", xlab = "linear predictors", main = "Binomial Fit")

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

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

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."  

Model selection using AIC via the `step` function, leads to a simpler model with 3 covariates, `acid`, `stage`, `xray`. 
```{r model choice}
step(Ex1018binom.glm)
```
We can also compare the fits of the full model and the simpler one using the deviance.
```{r compare}
Ex1018binom2.glm <- update(Ex1018binom.glm, .~. - aged - grade)
pchisq(deviance(Ex1018binom2.glm) - deviance(Ex1018binom.glm), df = 2, lower = F)
```