# R code for Crab example
# Poisson Regression Model for Count Data

library(rsq)
library(ggplot2)

# https://rdrr.io/cran/rsq/man/hcrabs.html
data(hcrabs)

head(hcrabs)
colnames(hcrabs) = c("C","S","W","Sa","Wt")

write.table(hcrabs, file='crabr.txt', col.names=FALSE)

# plot of satellite count versus width
ggplot(hcrabs, aes(x=W, y=Sa)) + geom_point()

# width as only predictor
model = glm(Sa ~ W, family=poisson(link=log), data=hcrabs)
summary(model)
model$fitted.values
rstandard(model)
anova(model)

# the linear predictor values
model$linear.predictors
exp(model$linear.predictors)

# to get the predicted count for each observation: 
print = data.frame(crab,pred=model$fitted)
print

#### Based on the residual deviance the model does NOT fit well
#### e.g., 567.88/171 = 3.3209

1-pchisq(model$deviance, model$df.residual)

#### Diagnostics measures (like in logistic regression)
#### But these work for ungrouped data too, 
#### as long as there is a variable with counts

influence(model)
plot(influence(model)$pear.res)
plot(model$linear.predictors, residuals(model, type="pearson"))

#### To predict a new value
newdt=data.frame(W=26.3)
predict.glm(model, type="response", newdata=newdt)

#### Let's assume for now that we do not have other covariates
#### and we will adjust for overdispersion 
#### first look at the sample mean and variances 
## e.g., tapply(crab$Sa, crab$W,function(x)c(mean=mean(x),variance=var(x)))

#### To estimate dispersion parameter
model.disp = glm(Sa ~ W, family=quasipoisson(link=log), data=hcrabs)
summary.glm(model.disp)
summary.glm(model.disp)$dispersion

#### Fit a negative binomial model
#### the dispersion parameter is THETA
#### Here are two different ways, both must use library(MASS) 
#### In the first one, we fix theta to be 1.0
#### In the second one we let glm.nb() estimate it

library(MASS)
nb.fit = glm(Sa~W, data=hcrabs, family=negative.binomial(theta=1, link="identity"), start=model$coef)
summary(nb.fit)

nb.fit1=glm.nb(Sa~W, data=crab, init.theta=1, link=identity, start=model$coef)
summary(nb.fit1)

#### Adding a categorical predictor
#### make sure C is a factor

# setting 5th color as baseline
hcrabs$C = factor(hcrabs$C, levels=c('5','2','3','4'))
model = glm(Sa ~ W+C, family=quasipoisson(link=log), data=hcrabs)
summary(model)
anova(model, test='Chisq')

# color as a numeric predictor 
data(hcrabs)
colnames(hcrabs) = c("C","S","W","Sa","Wt")
hcrabs$C = as.numeric(hcrabs$C)+1
model = glm(Sa ~ W+C, family=quasipoisson(link=log), data=hcrabs)
summary(model)
anova(model, test='Chisq')

# converting to 8 predictor groups, based on width intervals

# creating the 8 width groups
width.long = cut(hcrabs$W, breaks=c(0,seq(23.25, 29.25),Inf))

# summarizing per group
Cases = aggregate(hcrabs$W, by=list(W=width.long),length)$x
lcases = log(Cases)
Width = aggregate(hcrabs$W, by=list(W=width.long),mean)$x
AvgSa = aggregate(hcrabs$Sa, by=list(W=width.long),mean)$x
SdSa = aggregate(hcrabs$Sa, by=list(W=width.long),sd)$x
VarSa = aggregate(hcrabs$Sa, by=list(W=width.long),var)$x
SaTotal = aggregate(hcrabs$Sa, by=list(W=width.long),sum)$x

plot(Width, AvgSa)
ggplot(hcrabs.grp, aes(x=Width, y=AvgSa)) + geom_point()
cbind(table(width.long), Width, AvgSa, SdSa, VarSa)

hcrabs.grp = data.frame(Width=Width, lcases=lcases, SaTotal=SaTotal)
model.grp = glm(SaTotal ~ Width, offset=lcases, family=poisson, data=hcrabs.grp)
summary(model.grp)
rstandard(model.grp)
predict.glm(model.grp, type="response", newdata=hcrabs.grp)

# create a plot of the results
plot(SaTotal, pch="o", col="blue", main="Plot of Observed and Predicted Sa vs. groups", 
     xlab="Width groups", ylab="Number of Satellites")
points(model.grp$fitted, pch="p", col="red")
legend(6,30,c("obs","pred"), pch=c("o","p"), col=c("blue","red"))