Poisson regression

library(rjags)
library(tidyverse)
library(rjags)
library(ggmcmc)
library(polspline)
library(propagate)
d<-read.csv("https://tinyurl.com/aqm-data/marineinverts.csv")
str(d)

## 'data.frame':    45 obs. of  4 variables:
##  $ richness: int  0 2 8 13 17 10 10 9 19 8 ...
##  $ grain   : num  450 370 192 194 197 ...
##  $ height  : num  2.255 0.865 1.19 -1.336 -1.334 ...
##  $ salinity: num  27.1 27.1 29.6 29.4 29.6 29.4 29.4 29.6 29.6 29.6 ...

GGplot

Its easy to add a glm to a ggplot scatterplot. However be careful to add in the methods.args.

library(ggplot2)
g0<-ggplot(d,aes(x=grain,y=richness))
glm1<-g0+geom_point()+stat_smooth(method="glm",method.args=list( family="poisson"), se=TRUE) +ggtitle("Poisson regression with log link function")
glm1

pois_mod<-"
model{

#Likelihood
for( i in 1:n)
{
  y[i] ~ dpois(lambda[i])
  log(lambda[i]) <- a + b*x[i]

}

a~dnorm(0, tau1)
b~dnorm(0, tau2)
tau1 ~dunif(0,100)
tau2 ~dunif(0,100)

}"

dat<-list(
  n=length(d$richness),
  x=d$grain,
  y=d$richness)

mod <- jags.model(textConnection(pois_mod), data = dat, n.chains = 2)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 45
##    Unobserved stochastic nodes: 4
##    Total graph size: 227
## 
## Initializing model

update(mod, 10000)
dic1<- dic.samples(model = mod, variable.names = c("a","b","pd","deviance"), n.iter = 10000)

mod_sim <- coda.samples(model = mod, variable.names = c("a","b"), n.iter = 10000)

ms <-ggs(mod_sim) 
ms %>% spread(Parameter, value) -> ms


x<-seq(min(d$grain),max(d$grain),length=100)

pred<-function(x=x,mt=ms){
      pred<-exp(mt$a + mt$b * x)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
}
    
pred_mod<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
names(pred_mod)<-c("x","lwr","mid","upr")

g1<-g0+ geom_point() + geom_point()+stat_smooth(method="glm",method.args=list( family="poisson"), se=TRUE)

g1+  geom_line(data=pred_mod,aes(x=x,y=lwr),col="blue") + geom_line(data=pred_mod,aes(x=x,y=upr),col="red") + geom_line(data=pred_mod,aes(x=x,y=mid),col="black")

Over dispersed Poisson

od_pois_mod<-"
model{

#Likelihood
for( i in 1:n)
{
  
  y[i] ~ dgamma( ((mu[i]^2)/sig), (mu[i]/sig)) ## Reparameterised gamma
  mu[i] ~ dpois(lambda[i])
  log(lambda[i]) <- a + b*x[i]

}

#priors
a~dnorm(0, 0.00001)
b~dnorm(0, 0.00001)
sig~dunif(0,100)
}"

################ Not working

nb_mod<-"
model{

#Likelihood
for( i in 1:n)
{
  
  y[i] ~ dnegbin(p[i],r)
  p[i] <- r/(r+lambda[i]) 
  log(lambda[i]) <- a + b*x[i]

}

#priors
a~dnorm(0,0.00001)
b~dnorm(0, 0.00001)
r ~ dunif(0,50)
}"

dat<-list(
  n=length(d$richness),
  x=d$grain,
  y=d$richness)

mod <- jags.model(textConnection(nb_mod), data = dat, n.chains = 2)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 45
##    Unobserved stochastic nodes: 3
##    Total graph size: 310
## 
## Initializing model

update(mod, 10000)
dic2<- dic.samples(model = mod, variable.names = c("a","b","pd","deviance"), n.iter = 10000)
mod_sim <- coda.samples(model = mod, variable.names = c("a","b"), n.iter = 10000)
plot(mod_sim)

ms <-ggs(mod_sim) 
ms %>% spread(Parameter, value) -> ms


x<-seq(min(d$grain),max(d$grain),length=100)

pred<-function(x=x,mt=ms){
      pred<-exp(mt$a + mt$b * x)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
}
    
pred_mod<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
names(pred_mod)<-c("x","lwr","mid","upr")

#g0+ geom_point() + geom_point() +geom_smooth(method="glm.nb", se=TRUE) +
g1+ geom_line(data=pred_mod,aes(x=x,y=lwr),col="blue") + geom_line(data=pred_mod,aes(x=x,y=upr),col="red") + geom_line(data=pred_mod,aes(x=x,y=mid),col="black")

diffdic(dic2,dic1)

## Difference: -15.19085
## Sample standard error: 9.056654

sum(dic1$deviance)

## [1] 249.3316