---
title: "Toxoplasmosis"
author: "Nancy"
date: "07/12/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen=10, digits=4)
```

```{r libraries, include=FALSE}
library(faraway)
library(ggplot2)
library(dplyr)
library(kableExtra)
library(SMPracticals)
library(mgcv)
data(toxo)
```


## Logistic regression
Figure 10.12 in SM plots the data, here's a try at reproducing the right figure, along with 
linear, quadratic and cubic regression fits: 
```{r parametric}
toxo.glm1 <- glm(cbind(r,m-r) ~ rain, family = binomial, data = toxo)
toxo.glm2 <- glm(cbind(r,m-r) ~ rain + I(rain^2), family = binomial, data = toxo)
toxo.glm3 <- glm(cbind(r,m-r) ~ rain + I(rain^2) + I(rain^3), family = binomial, data = toxo)
with(toxo, plot(rain, r/m, pch=16, xlab = "Rainfall (mm)", ylab = "proportion positive"))
with(toxo, lines(rain,toxo.glm1$fitted.values))
with(toxo, lines(rain,toxo.glm2$fitted.values, lty=2, col="red"))
with(toxo, lines(rain,toxo.glm3$fitted.values, lty=3, col="blue"))
```


The linear and quadratic fits are nearly identical: 
```{r summary}
summary(toxo.glm3)
```

These coefficients are hard to read -- one of Cox \& Donelly's rules of thumb is that estimates should be on a comfortable scale.

```{r re-scale}
succ <- toxo$r
failure <- toxo$m - toxo$r
rainfall <- toxo$rain/100
summary(glm(cbind(succ,failure) ~ rainfall + I(rainfall^2) + I(rainfall^3), family = binomial))
```

Very strange that all the $z$-values are the same! There is some evidence of over-dispersion, let's try:
```{r quasi}
toxo.quasi <- glm(cbind(succ,failure) ~ rainfall + I(rainfall^2) + I(rainfall^3), family = quasibinomial)
summary(toxo.quasi)
```

The point estimates are the same, but the significance has reduced. I think the reason the coefficient estimates are all the same is that the correlation between rainfall and rainfall$^2$ and rainfall$^3$ is really high:
```{r check}
cor(rainfall,rainfall^2); cor(rainfall,rainfall^3)
```
This often happens with observational data. The function `poly` will construct *orthogonal* polynomials so the estimates become uncorrelated. You can often do almost as well by centering the $x$ variable first. 
```{r polynomial}
toxo$rainc <- toxo$rain/100 - mean(toxo$rain/100)
cor(toxo$rainc,toxo$rainc^2); cor(toxo$rainc,toxo$rainc^3)
toxo.glm4 <- glm(cbind(r,m-r) ~ rainc + I(rainc^2) + I(rainc^3), family = quasibinomial, data = toxo)

toxo.glm5 <- glm(cbind(r,m-r) ~ poly(rain/100, 3), family = quasibinomial, data = toxo)
summary(toxo.glm4)
summary(toxo.glm5)
plot(toxo$rain, toxo.glm5$fitted.values)
```

Davison goes on to fit local linear and local quadratic regression, using a version of `locpoly` for generalized linear models. It is explained on in Example 10.32. I used the `gam` package in `mgcv`. Using this is an art in itself, and the documentation isn't terribly helpful. Possibly better at this point to turn to *An Introduction to Statistical Learning* and follow their examples. The library `gam` also has a `gam` function, and it will do local polynomial fits using `gam.lo`..  

```{r smoothing}
toxo.gam <- gam(cbind(r, m-r) ~ s(rain), family = binomial, data = toxo)
summary(toxo.gam)
plot(toxo.gam)
toxo.gam2 <- gam(cbind(r, m-r) ~ s(rain,k=4), family = binomial, data = toxo) # k-1 is the approximate degrees of freedom used
summary(toxo.gam2)
with(toxo, plot(rain, r/m, pch=16, xlab = "Rainfall (mm)", ylab = "proportion positive", xlim = c(1600,2400)))
with(toxo, lines(rain,toxo.gam2$fitted.values))
```