---
title: "Poisson prediction interval"
author: "Will Townes"
output: html_document
---

Poisson prediction interval based on [Kim et al 2022](https://doi.org/10.1002/wics.1568)

```{r}
n<-100
x<-scale(1:n)
tot<-ceiling(rlnorm(n,3,.1))
offsets<-log(tot)
mu<-exp(offsets-x^2)
plot(x,mu)
y<-rpois(100,mu)
plot(x,y)
d<-data.frame(x=x,y=y,tot=tot)
fit<-glm(y~poly(x,2),family=poisson,data=d,offset=log(tot))
pred<-predict(fit,type="response",se.fit=TRUE)
plot(x,y,main="naive prediction interval")
lines(x,pred$fit)
lines(x,pred$fit-1.96*pred$se.fit,lty=2)
lines(x,pred$fit+1.96*pred$se.fit,lty=2)
```

```{r}
#gamma1 from Kim et al 2022
sf<-summary(fit)
Xo<-model.matrix(fit)
lam0<-pred$fit
#note: trace(V*xx') is equivalent to x'Vx
term1<-apply(Xo,1,function(x){x%*% crossprod(sf$cov.scaled, x)})
#term1<-diag(Xo%*% sf$cov.scaled %*% t(Xo)) #more expensive way to do same thing
Vo<- 1+(1)^2*lam0*term1
sqlamV<-sqrt(lam0*Vo)
lo<- lam0-1.96*sqlamV
hi<- lam0+1.96*sqlamV
plot(x,y,main="Gamma 1 prediction interval")
lines(x,lam0)
lines(x,lo,lty=2)
lines(x,hi,lty=2)
```