```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  results = 'markdown'
  )
```

```{r load-packages}
library(mosaic)
library(Stat2Data)
library(arm) # for binnedplot()
```

```{r load-data}
data(ICU)
## Don't run in Markdown
?ICU
```

We want to estimate patients' chance of surviving in the ICU, using Systolic Blood Pressure.

```{r}
lin.BP.model <- glm(Survive ~ SysBP, data = ICU, family = "binomial")
plotModel(lin.BP.model)
```

Check the binned residuals

```{r}
binnedplot(fitted(lin.BP.model), residuals(lin.BP.model, type = "pearson"), 
           nclass = 25)
```

Looks like there might be some nonlinearity here. Maybe we should try a quadratic?

```{r}
quad.BP.model <- 
  glm(Survive ~ SysBP + I(SysBP^2), family = "binomial", data = ICU)
plotModel(quad.BP.model)
```

```{r}
binnedplot(fitted(quad.BP.model), residuals(quad.BP.model, type = "pearson"), nclass = 15)
```

Not a whole lot better, to be honest.

```{r}
summary(quad.BP.model)
```

Based on the z test, the quadratic term seems to be helpful, though...

Let's look at the residuals for the first-order model:

```{r}
## Deviance Residuals
dresids <- residuals(lin.BP.model, type = "deviance")
histogram(dresids, breaks = 30)
```

```{r}
## Pearson Residuals
presids <- residuals(lin.BP.model, type = "pearson")
histogram(presids, breaks = 30)
```

Oof. That point way over on the left looks bad.

```{r}
plot(presids)
```

Negative 6 Pearson residual. Double oof.

Let's try to figure out what point that is...
```{r}
plot(presids)
## Don't run when Knitting
identify(presids)
```

Let's see what that is.
```{r}
ICU[196,]
```

Their Systolic BP is pretty high.  Where is it in the distribution?

```{r}
favstats(~SysBP, data = ICU)
gf_boxplot(~SysBP, data = ICU)
```

Not terribly surprising; it's the highest BP in the whole dataset. That's giving it a ton of **leverage**,
on top of the influence it has from being an outlier for Survival relative to BP.

Let's check out Cook's distance (Note: doesn't measure outlier status quite right with 
logistic regression, but it's ok as a rough tool)

```{r}
plot(lin.BP.model, which = 4)
```

Let's try leaving it out of the fit, in case it changes the picture dramatically:

```{r}
lin.BP.model.no.outlier <- 
  glm(Survive ~ SysBP, data = filter(ICU, ID != 921), family = "binomial")
plotModel(lin.BP.model.no.outlier)

quad.BP.model.no.outlier <- 
  glm(Survive ~ SysBP + I(SysBP^2), family = "binomial", data = filter(ICU, ID != 921))
plotModel(quad.BP.model.no.outlier  )
```

How are our residuals now?

```{r}
binnedplot(
  fitted(quad.BP.model.no.outlier), 
  residuals(quad.BP.model.no.outlier, type = "pearson"), 
  nclass = 15)
```

```{r}
presids <- residuals(quad.BP.model.no.outlier, type = "pearson")
histogram(presids, breaks = 30)
```


```{r}
binnedplot(
  fitted(lin.BP.model.no.outlier), 
  residuals(lin.BP.model, type = "pearson"), 
  nclass = 25)
```

```{r}
presids <- residuals(lin.BP.model.no.outlier, type = "pearson")
histogram(presids, breaks = 30)
```

