---
title: "Lecture 13: Binomial Models"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

```{r setup, include=FALSE}
# April 18th
knitr::opts_chunk$set(echo = FALSE)
```

# Before we start...

Number of Exercises: ```ceiling(N/2)```

This week's exercises don't have to be handed in until after Easter (I will out them up tomorrow)


# The Binomial Distribution

We have seen the Binomial distribution a few times, now we'll take it seriously...

We have $N$ tests, of which $r$ are a 'success'. Problem is to find $Pr(r|N)$

- may depend on covariates


# Likelihood

The probability of $r$ is 

$$
Pr(r|N, p) = \frac{N!}{r! (N-r)!} p^r (1-p)^{N-r}
$$

So the likelihood is 
$$
\begin{aligned}
l(p|N,r) & = r \log{p} + (N-r) \log (1-p) + \log{N!} - \log r! - \log(N-r)! \\
         & = r \left( \log p  - \log(1-p) \right) + N \log(1-p) + f(N, r) \\
         & = r \log \left( \frac{p}{1-p} \right) + N \log(1-p) + f(N, r)
\end{aligned}
$$

This is a GLM with a logit link function

- we will see some other link functions soon

# The logit Link

$$
\mu = \log \frac{p}{1-p}
$$
The inverse is

$$
p = \frac{e^\mu}{1 + e^\mu} = \frac{1}{1 + e^{-\mu}}
$$

```{r logitlink, fig.height=2, fig.width=4}
x <- seq(-5,5,length=50)
p <- exp(x)/(1 + exp(x))

par(mar=c(4.1,4.1,1,1))
plot(x, p, type = "l", xlab=expression(mu), las=1)


```

# Interpreting the logit

The logit is the log-odds

Odds: p/(1-p)

Tiger Roll won the Grand National on Saturday with an odds of 10:1

- if I bet £1, I win £10
- if the odds are fair, this means that for every Grand National won by Tiger Roll, he would lose 10.

The probability is Success/(Success + Failure), i.e. 

$$1/(1 + 10) = 1/11 = 0.09$$

# More interpretation of the Logit

if we have a basline effect  $\alpha$ (so $p=e^{\alpha}/(1 + e^{\alpha})$) and we increase it by a small $\beta$ (i.e. $\beta = 0 + \varepsilon$, so $e^\beta \approx \varepsilon$) then 

$$
p = \frac{e^{\alpha + \beta}}{1+e^{\alpha + \beta}} =  \frac{e^{\alpha} e^{\beta}}{1+e^{\alpha}e^{\beta}} \approx \frac{e^{\alpha} }{1+e^{\alpha}} e^{\beta} \approx \frac{e^{\alpha} }{1+e^{\alpha}} (1+\beta)
$$

So a small effect (approximately) adds the effect to the probability 

 - especially if $p$ is small

# More interpretation of the Logit

e.g. if $\mu = -3$, 

$$p = \frac{e^{-3}}{(1+e^{-3})} = `r round(exp(-3)/(1+exp(-3)), 3)`$$

Now let $\beta = 0.2$, 

$$p = \frac{e^{-3+0.2}}{(1+e^{-3+0.2})} = `r round(exp(-3+0.2)/(1+exp(-3+0.2)), 3)`$$

and $0.047 \times e^{0.2} = 0.047 \times `r round(exp(0.2), 3)` = `r round(0.047*exp(0.2), 3)`$ 

# An example

```{r Data}
# 6882
# Dryocopus_pileatus
library(raster)
if(!file.exists("../Data/Dpileatus.csv")) {
  library(rBBS)
  library(cruts)
  library(maps)
  library(maptools)

  fityear <- 2010
  TimeRange <- c(paste0(fityear-30,"-06-01"), paste0(fityear,"-05-31"))
  proj <- CRS("+proj=longlat +ellps=WGS84") # Single projection

# Get sparrow data as spatial data frame
  # House sparrow: 6882
  # White-crowned sparrow: 5540
  load("../Data/Dryocopus_pileatusBBS.RData")
  Dpil <- Data$BBS
  

  # Get US polygon
  USmap <- map(database="usa", interior=FALSE, resolution=0, regions="main") # map
  USpolygon <- map2SpatialPolygons(USmap, ID="US", proj4string = proj)

# Get CRU precipitation data
  file.copy(from="~/git-repos/Map-of-Life-Models/Data/cru_ts3.24.01.1901.2015.pre.dat.nc.gz", to="CRUtmp.gz")
  CRUfile <- R.utils::gunzip("CRUtmp.gz")
  fitdat <- cruts2raster(CRUfile, timeRange = TimeRange, poly = USpolygon)

  Dpil[[paste0("prec.mean")]] <- raster::extract(mean(fitdat), y=Dpil)
  Dpil[[paste0("prec.seas")]] <- raster::extract(calc(fitdat, fun=sd)/mean(fitdat), y=Dpil)
  
  file.remove("CRUtmp")

# Get CRU temperature data
  file.copy(from="~/git-repos/Map-of-Life-Models/Data/cru_ts3.24.01.1901.2015.tmp.dat.nc.gz", to="CRUtmp.gz")
  CRUfile <- R.utils::gunzip("CRUtmp.gz")
  fitdat <- cruts2raster(CRUfile, timeRange = TimeRange, poly = USpolygon)

  Dpil[["temp.mean"]] <- raster::extract(mean(fitdat), y=Dpil)
  Dpil[["temp.seas"]] <- raster::extract(calc(fitdat, fun=sd)/mean(fitdat), y=Dpil)

  file.remove("CRUtmp")

# Tidy up data and save
  Dpil <- Dpil[!is.na(Dpil[["temp.mean"]]),]
  Dpil.df <- cbind(Dpil@coords, Dpil@data)
  
  write.csv(Dpil.df, file="../Data/Dpileatus.csv")
  save(Data, Dpil, VarsToKeep, file="../Data/Dpileatus.RData")
  rm(CRUfile, fitdat, fityear, Dpil, proj, TimeRange, USmap, USpolygon, VarsToKeep)
} else {
  load(file="../Data/Dpileatus.RData")
  DryoPil <- read.csv(file="../Data/Dpileatus.csv")
}

DryoPil$prec.mean.sc <- as.vector(scale(DryoPil$prec.mean))
DryoPil$temp.mean.sc <- as.vector(scale(DryoPil$temp.mean))
```

The data comne from the North American Breeding Bird survey. We have presence of the Pileated woodpecker (_Dryocopus pileatus_)

- a toal of `r nrow(DryoPil)` routes, at each one there are 50 stops (so 50 trials)
 
Can we explain its distribution with rain & temperature?


![By Joshlaymon (CC BY-SA 3.0), from Wikimedia Commons](PileatedWoodpeckerFeedingonTree.jpg)


# The data

Red: recorded in 2010, grey: not recorded in 2010

```{r PlotPileated}
plot(Dpil, col=c("tan", "red")[1 + (Dpil$NPres>0)], pch=16, cex=1)
```


# Fitting the model

First, we will ignore the 50 trials, and just look at whether the pileated woodpecker was seen at least one on each route

There are several ways to specify the response, depending on the data

- as a factor (first level is failure, rest is success)
- as 2 columns, with successes and failures

# Fitting the model: Method 1

As a logical vector

```{r LogicalGLM, echo=TRUE, comment=NA}
DryoPil$Present <- DryoPil$NPres>0
DryoPil$PresentF <- factor(DryoPil$Present)
DryoPil$PresentF[7:11]

mod.method1 <- glm(PresentF ~ temp.mean.sc, data=DryoPil, 
                family="binomial")
```


# Fitting the model: Method 2

As success and failure columns

```{r cbindGLM, echo=TRUE, comment=NA}
DryoPil$Absent <- 1-DryoPil$Present
cbind(DryoPil$Present, DryoPil$Absent)[7:11,]
mod.method2 <- glm(cbind(Present, Absent) ~ temp.mean.sc, 
                data=DryoPil, family="binomial")

```

# Fitting the model: Method 2

The advantage of this method is that it allows for a varying number of trials

```{r cbindGLMAll, echo=TRUE, comment=NA}
DryoPil$NAbs <- DryoPil$Ntrials-DryoPil$NPres
cbind(DryoPil[7:11,c("NPres", "NAbs")])
mod.mod.method2all <- glm(cbind(NPres, NAbs) ~ temp.mean.sc, 
                data=DryoPil, family="binomial")
```

# Model results

Methods 1 & 2 provide identical results

- just different ways of writing the same thing

```{r SommMethod1, echo=TRUE, comment=NA}
round(summary(mod.method1)$coefficients, 3)
```

```{r SommMethod2, echo=TRUE, comment=NA}
round(summary(mod.method2)$coefficients, 3)
```

# Model Interpretation: intercept

```{r Coefs, echo=FALSE}
odds <- coef(mod.method1)["(Intercept)"]
slope <- coef(mod.method1)["temp.mean.sc"]

```

```{r SommMethod1a, echo=TRUE, comment=NA}
round(summary(mod.method1)$coefficients, 3)

```

Mean temperature is scaled, so at the mean temperature ($`r round(mean(DryoPil$temp.mean), 1)`^\circ C$), the log odds of the species being observed is `r round(coef(mod.method1)["(Intercept)"], 2)`. This is the same as a probability of 

$$
\frac{e^{`r round(odds, 2)`}}{1 + e^{`r round(odds, 2)`}} = `r round( exp(odds)/(1+exp(odds)), 2)`
$$

# Model Interpretation: slope

if we change the mean temperature by 1 standard deviation ($`r round(sd(DryoPil$temp.mean), 1)`^\circ C$), the probability of observed presence is

$$
\frac{e^{`r round(odds, 2)` + `r round(slope, 2)` \times 1}}{1 + e^{`r round(odds, 2)` + `r round(slope, 2)` \times 1}} = `r round(exp(odds + slope)/(1 + exp(odds + slope)), 2)`
$$

# Model Interpretation: Response

We can plot the whole curve over the range of the data:

```{r WholeCurve, echo=FALSE, fig.height=4, fig.width=4}
PredData <- data.frame(temp.mean.sc = seq(min(DryoPil$temp.mean.sc), 
                                          max(DryoPil$temp.mean.sc), length=50))
Prob <- predict(mod.method1, newdata = PredData, 
                         type="response", se.fit=TRUE)
Prob$lower <- Prob$fit - Prob$se.fit
Prob$upper <- Prob$fit + Prob$se.fit
plot(PredData$temp.mean.sc, Prob$fit, type="n", las=1, ylim=c(0,1), 
     xlab=expression("Scaled Mean Temperature " ( degree~C)), ylab="Probability of Presence")
polygon(c(PredData$temp.mean.sc, rev(PredData$temp.mean.sc)), c(Prob$lower, rev(Prob$upper)),
        col="grey70", border=NA)
lines(PredData$temp.mean.sc, Prob$fit)

```

# Analysis of Deviance

Let's fit a bigger model

```{r GLMBigger, echo=TRUE, comment=NA}
mod.big <- glm(Present ~ prec.mean.sc + temp.mean.sc + I(prec.mean.sc^2) + I(temp.mean.sc^2), 
               data=DryoPil, family="binomial")
an.big <- paste(capture.output(print(anova(mod.big, test="Chisq"))), "\n")
cat(an.big[10:17])

```

Suggests a quadratic effect of precipitation (prec.mean.sc) and linear effect of temperature. 


# A better model

```{r SummGLMBigger, echo=TRUE, comment=NA}
mod.better <- glm(Present ~ prec.mean.sc + 
                    I(prec.mean.sc^2) + temp.mean.sc,
               data=DryoPil, family="binomial")
round(summary(mod.better)$coefficients, 3)

```

The quadratic term for precipitation is negative, so there is a maximum. 

But the temperature effect has reversed sign

# A better Prediction


```{r BetterWholeCurve, echo=FALSE, fig.height=4, fig.width=6}
PredData.prec <- data.frame(prec.mean.sc = seq(min(DryoPil$prec.mean.sc), 
                                          max(DryoPil$prec.mean.sc), length=50),
                            temp.mean.sc = 0)
ProbBetter.prec <- predict(mod.better, newdata = PredData.prec, 
                         type="response", se.fit=TRUE)
ProbBetter.prec$lower <- ProbBetter.prec$fit - ProbBetter.prec$se.fit
ProbBetter.prec$upper <- ProbBetter.prec$fit + ProbBetter.prec$se.fit

PredData.temp <- data.frame(temp.mean.sc = seq(min(DryoPil$temp.mean.sc), 
                                          max(DryoPil$temp.mean.sc), length=50),
                            prec.mean.sc = 0)
ProbBetter.temp <- predict(mod.better, newdata = PredData.temp, 
                         type="response", se.fit=TRUE)
ProbBetter.temp$lower <- ProbBetter.temp$fit - ProbBetter.temp$se.fit
ProbBetter.temp$upper <- ProbBetter.temp$fit + ProbBetter.temp$se.fit

par(mfrow=c(1,2), mar=c(4.1,4,1,0), oma=c(0,0,0,0))
plot(PredData.prec$prec.mean.sc, ProbBetter.prec$fit, type="n", las=1, ylim=c(0,1), 
     xlab="Scaled Mean Precipitation", ylab="Probability of Presence")
polygon(c(PredData.prec$prec.mean.sc, rev(PredData.prec$prec.mean.sc)), 
        c(ProbBetter.prec$lower, rev(ProbBetter.prec$upper)),
        col="grey70", border=NA)
lines(PredData.prec$prec.mean.sc, ProbBetter.prec$fit)

plot(PredData.temp$temp.mean.sc, ProbBetter.temp$fit, type="n", las=1, ylim=c(0,1), 
     xlab="Scaled Mean temperature", ylab="")
polygon(c(PredData.temp$temp.mean.sc, rev(PredData.temp$temp.mean.sc)), 
        c(ProbBetter.temp$lower, rev(ProbBetter.temp$upper)),
        col="grey70", border=NA)
lines(PredData.temp$temp.mean.sc, ProbBetter.temp$fit)

```

The maximum for precipitation is in the range of the data (this is not always the case!)

# Link Functions

Unlike the Poisson, the binomial distribution has 3 link functions that are used:

- logit
- probit
- cloglog

# logit

We've seen above that this is the natural link function 

- "canonical link"

# Probit: a threshold model

Imagine we have a normal distribution

```{r probitlink, fig.height=4}

xx <- seq(-5,5,length=100)
zz <- dnorm(xx, mean=0.2)
plot(xx, zz, type="l", yaxs="i", ylim=c(0, max(zz*1.04)), ann=FALSE)
x.poly <- c(0, max(xx), rev(xx[xx>0]))
y.poly <- c(0, 0, dnorm(rev(xx[xx>0]), mean=0.2))
polygon(x.poly, y.poly, border=NA, col="hotpink")
lines(xx, zz)
        
        
```

if the distribution is >0 we have a success, if it is < 0 we have a failure

- the large the mean the greater the probability of success
- but still a random chance

This is the same as a probit link

- the inverse normal distribution

# cloglog

Last week we looked at the Poisson distribution

Sometimes we have presence/absence for something that is really a count 

# Dilution assays

We take a sample that might be contaminated by a microorganism

We serially dilute the sample 

- concentrations $x, x/2, x/4, x/8, \dots x/2^n$

Streak out onto agar plates. See if anything grows

If a sample contains the microorganism, it will grow. Assume a Poisson distribution for the organism, then if we had a count of the organism we would model it with a Poisson distribution and log link. But we only have presence/absence. 


# c the log log
$$
p = Pr(growth) = Pr(n > 0) = 1-Pr(n = 0) = 1-\frac{\lambda_i^0 e^{-\lambda}}{0!} = 1 - e^{-\lambda}
$$

So, 

$$
\begin{aligned}
p &= 1- e^{-\lambda} \\
e^{-\lambda} &= 1 - p \\
\lambda &= - \log(1-p) \\
\log(\lambda) &= \log(-\log(1-p))
\end{aligned}
$$

and $\log(-\log(1-p))$ is the cloglog link function.

# The link functions

logit & probit almost the same

cloglog asymmetrical

```{r Linkfunctions}
N <- 1e3
x <- seq(-5,5, length=100)
p <- exp(x)/(1+exp(x))
r <- rbinom(length(x), N, p)
resp <- cbind(r, N-r)
lgt.mod <- glm(resp ~ x, family=binomial("logit"))
probit.mod <- glm(resp ~ x, family=binomial("logit"))
cloglog.mod <- glm(resp ~ x, family=binomial("cloglog"))

plot(x, predict(lgt.mod, type="response"), type="l", lwd=4, ylab="Probability", yaxs="i")
lines(x, predict(probit.mod, type="response"), col=2, lwd=4, lty=3)
lines(x, predict(cloglog.mod, type="response"), col=4, lwd=4)
legend(-4.2, 0.73,c("logit", "probit", "cloglog"), lty=c(1,3,1), col=c(1,2,4), lwd=4)
```

# Symmetry

For logit & probit

$$
\begin{aligned}
Pr(success | x=0) &= 0.5 \\
Pr(success | x-\mu) &= Pr(failure | -(x-\mu))
\end{aligned}
$$


```{r Symmetry, fig.height=4, fig.width=5}
x1 <- c(-10,-2, -2)
plot(x, predict(lgt.mod, type="response"), type="l", lwd=2, ylab="Probability", yaxs="i", yaxt="n")
lines(x, predict(probit.mod, type="response"), col=2, lwd=2, lty=3)
lines(x, predict(cloglog.mod, type="response"), col=4, lwd=2)
# lines(x1, c(predict(probit.mod, newdata=data.frame(x=rep(x1[2],2)), type="response"), 0), col=2, lwd=2, lty=1)
# lines(c(x1[1],-x1[-1]), c(predict(probit.mod, newdata=data.frame(x=rep(-x1[2],2)), type="response"), 0), col=2, lwd=2, lty=1)
# lines(c(-10,0,0), c(0.5,0.5,0), col=1, lwd=2, lty=1)
# # legend(-4.2, 0.73,c("logit", "probit", "cloglog"), lty=c(1,3,1), col=c(1,2,4), lwd=2, bg="white")
axis(2, at=c(0.25,0.5, 0.75), las=1)


```

# Symmetry

For logit & probit

$$
\begin{aligned}
Pr(success | x=0) &= 0.5 \\
Pr(success | x-\mu) &= Pr(failure | -(x-\mu))
\end{aligned}
$$


```{r Symmetry2, fig.height=4, fig.width=5}
plot(x, predict(lgt.mod, type="response"), type="l", lwd=2, ylab="Probability", yaxs="i", yaxt="n")
lines(x, predict(probit.mod, type="response"), col=2, lwd=2, lty=3)
lines(x, predict(cloglog.mod, type="response"), col=4, lwd=2)
lines(x1, c(predict(probit.mod, newdata=data.frame(x=rep(x1[2],2)), type="response"), 0), col=2, lwd=1.2, lty=1)
lines(c(x1[1],-x1[-1]), c(predict(probit.mod, newdata=data.frame(x=rep(-x1[2],2)), type="response"), 0), col=2, lwd=1.2, lty=1)
lines(c(-10,0,0), c(0.5,0.5,0), col=1, lwd=1.2, lty=1)
# legend(-4.2, 0.73,c("logit", "probit", "cloglog"), lty=c(1,3,1), col=c(1,2,4), lwd=2, bg="white")
axis(2, at=c(0.25,0.5, 0.75), las=1)


```

# When to use the different link functions

logit

- default. Usually makes sense

probit

- sometimes with mixed models it's easier to understand. Otherwise, use a logit

cloglog

- when you have Poisson-like counts 
- not all count data!
    - sometimes it is not the counting that dominates
    


# A better model, different links

```{r SummGLMLinks, echo=TRUE, comment=NA}
fm <- Present ~ prec.mean.sc+I(prec.mean.sc^2)+temp.mean.sc
mod.logit <- glm(fm, data=DryoPil, family=binomial("logit"))
round(mod.logit$coefficients, 3)
mod.probit <- glm(fm, data=DryoPil, family=binomial("probit"))
round(mod.probit$coefficients, 3)
mod.cloglog <- glm(fm, data=DryoPil, family=binomial("cloglog"))
round(mod.cloglog$coefficients, 3)

```

# A better model, different links

```{r Betterlinks, echo=FALSE, fig.height=4, fig.width=6}
GetPreds <- function(data, model) {
  pred <- predict(model, newdata = data, type="response", se.fit=TRUE)

  pred$lower <- pred$fit - pred$se.fit
  pred$upper <- pred$fit + pred$se.fit
  cbind(data, pred)
}
pred.prec.logit <- GetPreds(data=PredData.prec, model=mod.logit)
pred.temp.logit <- GetPreds(data=PredData.temp, model=mod.logit)

pred.prec.probit <- GetPreds(data=PredData.prec, model=mod.probit)
pred.temp.probit <- GetPreds(data=PredData.temp, model=mod.probit)

pred.prec.cloglog <- GetPreds(data=PredData.prec, model=mod.cloglog)
pred.temp.cloglog <- GetPreds(data=PredData.temp, model=mod.cloglog)


par(mfrow=c(1,2), mar=c(4.1,4,1,0), oma=c(0,0,0,0))
plot(pred.prec.logit$prec.mean.sc, pred.prec.logit$fit, type="n", las=1, ylim=c(0,1), 
     xlab="Scaled Mean Precipitation", ylab="Probability of Presence")
polygon(c(pred.prec.logit$prec.mean.sc, rev(pred.prec.logit$prec.mean.sc)), 
        c(pred.prec.logit$lower, rev(pred.prec.logit$upper)),
        col="grey90", border=NA)
lines(pred.prec.logit$prec.mean.sc, pred.prec.logit$fit)
lines(pred.prec.probit$prec.mean.sc, pred.prec.probit$fit, col=2)
lines(pred.prec.cloglog$prec.mean.sc, pred.prec.cloglog$fit, col=4)


plot(pred.temp.logit$temp.mean.sc, pred.temp.logit$fit, type="n", las=1, ylim=c(0,1), 
     xlab="Scaled Mean temperature", ylab="")
polygon(c(pred.temp.logit$temp.mean.sc, rev(pred.temp.logit$temp.mean.sc)), 
        c(pred.temp.logit$lower, rev(pred.temp.logit$upper)),
        col="grey70", border=NA)
lines(pred.temp.logit$temp.mean.sc, pred.temp.logit$fit)
lines(pred.temp.probit$temp.mean.sc, pred.temp.probit$fit, col=2)
lines(pred.temp.cloglog$temp.mean.sc, pred.temp.cloglog$fit, col=4)

legend(-1, 0.8, c("logit", "probit", "cloglog"), lty=1, col=c(1,2,4))
```


# Model Checking

With binary data residuals are useless

```{r BinaryResids}
plot(fitted(mod.better), resid(mod.better))
```

# Model Checking: grouping residuals

But we can group residuals and take the mean

```{r GroupedBinaryResids, echo=TRUE, fig.height=5}
# Split data into groups, and calc. mean of residuals
FitCut <- cut(fitted(mod.better), breaks=20, labels=FALSE)
ResGrp <- tapply(resid(mod.better), list(FitCut), mean)
FitGrp <- tapply(fitted(mod.better), list(FitCut), mean)
plot(FitGrp, ResGrp, xlab="Fitted (mean)", ylab="Mean residual", cex=1.4)

```



# Model Checking

Better when we have a larger number of trials

```{r BetterBinaryResids, echo=TRUE, comment=NA, fig.height=4}
mod.method3 <- glm(cbind(NPres, NAbs) ~ prec.mean.sc, 
                   data=DryoPil, family="binomial")
par(mfrow=c(1,2))
plot(fitted(mod.method3), resid(mod.method3))
plot(DryoPil$prec.mean.sc, resid(mod.method3))

```

Linear model; residuals negative at both extremes for precipitation



# Model Checking

Add a quadratic term...

```{r BetterBinaryResidsQuad, echo=TRUE, comment=NA, fig.height=4}
form <- cbind(NPres, NAbs)~prec.mean.sc+I(prec.mean.sc^2)
mod.method3 <- glm(form, data=DryoPil, family="binomial")
par(mfrow=c(1,2))
plot(fitted(mod.method3), resid(mod.method3))
plot(DryoPil$prec.mean.sc, resid(mod.method3))

```


# Overdispersion 

In the binomial, like the Poisson, the variance is controlled by the mean.

So if we have more variation than expected, we can have overdispesion

We can test for it in the same way.

# Overdispersion

Without, then with overdispersion

```{r Overdispersion, echo=TRUE} 
N <- 100
x.all <- seq(-5,5,length=50)
p1 <- exp(x.all)/(1+exp(x.all))
r1 <- rbinom(length(p1), size=N, p1)
mod1 <- glm(cbind(r1, N-r1) ~ x.all, family="binomial")

x2 <- rnorm(length(x.all), x.all, 1)
p2 <- exp(x2)/(1+exp(x2))
r2 <- rbinom(length(p2), size=N, p2)
mod2 <- glm(cbind(r2, N-r2) ~ x.all, family="binomial")
```

# Overdispersion

```{r PlotOverdispersion, echo=TRUE, fig.height=4} 
par(mfrow=c(1,2))
plot(x.all,r1/N, main="Without overdispersion")
lines(x.all, predict(mod1, type="response"))
plot(x.all,r2/N, main="With overdispersion")
lines(x.all, predict(mod2, type="response"))
```


# Overdispersion 

Without Overdispersion

```{r NoOverdispersionDev, echo=FALSE, comment=NA} 
cat(" Deviance = ", deviance(mod1), " with ", df.residual(mod1), "DF.\np is", pchisq(deviance(mod1), df.residual(mod1), lower.tail = FALSE),  " \n")
cat("Deviance ratio = ", round(deviance(mod1)/ df.residual(mod1), 2), "\n")
```

With Overdispersion

```{r OverdispersionDev, echo=FALSE, comment=NA} 
cat("Deviance = ", deviance(mod2), " with ", df.residual(mod2), "DF.  \np is", pchisq(deviance(mod2), df.residual(mod2), lower.tail = FALSE),  " \n")
cat("Deviance ratio = ", round(deviance(mod2)/ df.residual(mod2), 2), "\n")

```

# Sparseness

Overdisperson can't always be detected. If we only have 0 s and 1s, then the overdispersion affects the model

```{r Sparseness, echo=TRUE} 
x.all <- seq(-4,4,length=1e4)
p1 <- exp(x.all)/(1+exp(x.all))
r1 <- rbinom(length(p1), size=1, p1)
mod1 <- glm(r1 ~ x.all, family="binomial")

x2 <- rnorm(length(x.all), x.all, 5)
p2 <- exp(x2)/(1+exp(x2))
r2 <- rbinom(length(p2), size=1, p2)
mod2 <- glm(r2 ~ x.all, family="binomial")

```

# Plot Sparseness

```{r PlotSparseness, echo=TRUE, fig.height=4} 
par(mfrow=c(1,2))
plot(x.all,jitter(r1), cex=0.2, col="grey70")
lines(x.all, predict(mod1, type="response"), col=2, lwd=1.5)
plot(x.all,jitter(r2), cex=0.2, col="grey70")
lines(x.all, predict(mod2, type="response"), col=2, lwd=1.5)
```

Can't tell if we have a poor model or overdispersion!

# Next week

Summary of the course 

- this will give me time to write a curriculum etc.
