### Looking at the fish data.


### First the daily sample sizes.

n <- c(1000,912,789,946,932,1045,997,789,956,1100)

### Now the number of dead fish.

dead <- c(948,841,737,906,897,1006,959,742,873,1045)

### Here is a point estimate of the proportion dead (assuming an
### a constant probability)

thetahat <- sum(dead)/sum(n)

### And a classical binomial confidence interval.

thetahat.low  <- 
         thetahat + qnorm(0.025) * sqrt(thetahat * (1 - thetahat) / sum(n))
thetahat.high <- 
         thetahat + qnorm(0.975) * sqrt(thetahat * (1 - thetahat) / sum(n))

### Find out individual confidence intervals for each day based on 
### this thetahat.


dayhat.low  <- 
        n * thetahat + qnorm(0.025) * sqrt(thetahat * (1 - thetahat) * n)
dayhat.high <- 
        n *  thetahat + qnorm(0.975) * sqrt(thetahat * (1 - thetahat) * n)


### Plot this information

motif()

plot(n, n * thetahat)
lines(n, dayhat.low,col=3)
lines(n, dayhat.high,col=3)
points(n,dead,pch = 2)


### Comments on the plot.



### What is the predicted variance from the model for an individula thetahat? 

predvar <- thetahat * (1 - thetahat) / n

### What is the observed variance? (This one is based on the data, 
### NOT the model)

obsvar <-  var(dead/n)


### Now find the ration of these two, the excess variance.

mean(sqrt(obsvar/predvar))



### Some models for this data:

glm.out <- glm(cbind(dead, n - dead) ~ 1,family = binomial)


### Summarize this one.

summary(glm.out)


### Try this: Residual deviance / Resid df
sqrt(49.16208/9)

### Comments.

### Note exp(2.861531)/(1 + exp(2.861531))

upper <-  2.861531 + qnorm(.025) * 0.1056719
exp(upper)/(1 + exp(upper))