## https://bookdown.org/drki_musa/dataanalysis/poisson-regression.html

rm(list=ls())

### Dunn & Smith residuals for Poi
Resid_Poisson = function(X,Y,beta){
     X = as.matrix(X)
     n = length(Y)
     one = rep(1,n)
     X = cbind(one,X)
     beta = as.matrix(beta)
     L = exp(X %*% beta)
     p = ppois(Y,L)     ##### prob(Poisson) lower tail
     b = p
     a = ppois(Y-1,L)
     U = runif(n,a,b)
     r = qnorm(U)
     return(r)
     }



D = read.csv("asthma.csv",head=TRUE)
attach(D)

## names(D)
## [1] "gender"  "res_inf" "ghq12"   "attack" 


### Poisson Estimate One Feature: Gen. Health  

out = glm(attack ~ ghq12,data=D,family="poisson")

## summary(out)
## Call:
## glm(formula = attack ~ ghq12, family = "poisson", data = D)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.230923   0.159128  -1.451    0.147    
## ghq12        0.059500   0.006919   8.599   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 229.56  on 119  degrees of freedom
## Residual deviance: 145.13  on 118  degrees of freedom
## AIC: 422.2
## 
## Number of Fisher Scoring iterations: 5

beta=out$coefficients

## print(beta)
## (Intercept)       ghq12 
## -0.23092329  0.05950017 

o = order(ghq12)
fun = exp(beta[1]+beta[2]*ghq12[o])

pdf("PoissonEstimate.pdf")
plot(ghq12[o],fun,type="l",col="purple",lwd=4,ylab="",xlab="")
points(out$y,pch=21,bg="pink",cex=1.4,col="purple")
title(ylab="Asthma attacks", line=2, cex.lab=2.5)
title(xlab="Score Gen. Health",line=3, cex.lab=2.3)
dev.off()



### Poisson Residuals -- Deviance

r = residuals(out,type="deviance")

pdf("PoissonDeviance.pdf")
plot(ghq12,r,xlab="",ylab="",ylim=c(-3.2,3.2),pch=21,bg="pink",cex=1.8)
abline(h=0,col="red",lwd=3)
title(ylab="Residuals", line=2.6, cex.lab=2)
title(xlab="Score Gen. Health",line=3,cex.lab=2)
dev.off()


### Randomized Quantile Residuals

Y = attack
n = length(Y)
x3 = as.numeric(D[,3])
X = x3

pdf("PoissonDunnSmith.pdf")
r = Resid_Poisson(X,Y,beta)
plot(x3,r,xlab="",ylab="",ylim=c(-3.2,3.2),pch=21,cex=1.8,bg="pink")
abline(h=0,col="red",lwd=3)
title(ylab="Residuals", line=2.6, cex.lab=2)
title(xlab="Score Gen. Health.",line=3, cex.lab=2.3)
dev.off()