source("http://192.38.117.59/~linearpredictors/datafiles/readFever.R")

##  NOTE: The row corresponding to Binomial in Table 7.2.5 is actually the
##  results corresponding to a log-link. The traditional logit-link gives
##  very slightly different results. They are both calculated in this program

fever$parity <- relevel(fever$parity, ref = "1")

summary(pois <- glm(episodes ~ parity + I((mage - 30)/10), data = fever,
                    family = poisson))

## prediction + CI
p.pois <- predict(pois, newdata = data.frame(mage = 30, parity = factor(1)),
                  type = "response",se.fit = T)
p.pois$fit + c(Prediction = 0, LCL = -1.96, UCL = 1.96)*p.pois$se.fit


summary(binlog <- glm(episodes/14 ~ parity + I((mage - 30)/10), data = fever,
                    family = binomial(link = "log")))

## prediction + CI
p.binlog <- predict(binlog, newdata = data.frame(mage = 30, parity = factor(1)),
                  type = "response",se.fit = T)
14*(p.binlog$fit + c(Prediction = 0, LCL = -1.96, UCL = 1.96)*p.binlog$se.fit)

summary(bin <- glm(episodes/14 ~ + parity + I((mage - 30)/10), data = fever,
           family = binomial(link = "logit")))
## prediction + CI
p.bin <- predict(bin, newdata = data.frame(mage = 30, parity = factor(1)),
                  type = "response",se.fit = T)
14*(p.bin$fit + c(Prediction = 0, LCL = -1.96, UCL = 1.96)*p.bin$se.fit)

## the glm procedure cannot find start values itself in this particular case,
## so we provide the estimates from the poisson fit as start values.
summary(normal <- glm(episodes ~ parity + I((mage - 30)/10), data = fever,
              family = gaussian(link = "log"), start = pois$coef))
p.norm <- predict(normal, newdata = data.frame(mage = 30, parity = factor(1)),
                  type = "response",se.fit = T)
(p.norm$fit + c(Prediction = 0, LCL = -1.96, UCL = 1.96)*p.norm$se.fit)