#### #### #### Part I: Preparing the workspace and exploring the data #### #### ####
rm(list=ls())

#getwd()
#setwd("pathway")
library(pscl)
library(bbmle)
library (MASS)

dataforstats <- read.table("dataforstats_final_BB121613.txt", header=T)
names (dataforstats)
dataforstats$depth2 <- dataforstats$depth^2
dataforstats$ranchn <- factor(dataforstats$ranchno)

par(mfrow = c(2,4))
m0 <- mean(dataforstats$fishct)
m1 <- mean(dataforstats$fishct[dataforstats$fishct >0])
Poiss_s0 <- rpois(length(dataforstats$fishct),m0)
Poiss_s1 <- rpois(length(dataforstats$fishct),m1)
Neg_bin_s0 <- rnbinom(length(dataforstats$fishct),mu=m0,size=0.58)
Neg_bin_s1 <- rnbinom(length(dataforstats$fishct),mu=m1,size=0.58)
Neg_bin_s0_100 <- rnbinom(length(dataforstats$fishct),mu=m0,size=100)
Neg_bin_s1_100 <- rnbinom(length(dataforstats$fishct),mu=m0,size=0.19)

b=seq(0,max(dataforstats$fishct)+30,1)
hist(dataforstats$fishct,breaks = b,ylim=c(0,700),main="all data", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(dataforstats$fishct[dataforstats$fishct>0],breaks = b,ylim=c(0,700),main="all data >0", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Poiss_s0,breaks = b,ylim=c(0,700),main="Poisson", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Poiss_s1,breaks = b,ylim=c(0,700),main="Poisson no zeros", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Neg_bin_s0,breaks = b,ylim=c(0,700),main="Negative Binomial 0.58", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Neg_bin_s1,breaks = b,ylim=c(0,700),main="Neg Bin no zeros 0.58", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Neg_bin_s0_100,breaks = b,ylim=c(0,700),main="Negative Binomial 100", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")
hist(Neg_bin_s1_100,breaks = b,ylim=c(0,700),main="Negative Binomial 0.19", xlim =c(0,20))
abline(v=m0,col="red")
abline(v=m1,col="blue")

#### #### #### Part II: Model building and selection  #### #### ####

f1 <- formula(fishct ~ depth * ranchn + depth2 * ranchn)
f2 <- formula(fishct ~ depth * ranchn + depth2 * ranchn | depth * ranchn + depth2 * ranchn)
f3 <- formula(fishct ~ depth * ranchn + depth2 * ranchn | depth)
f4 <- formula(fishct ~ depth * ranchn + depth2 * ranchn | ranch)
f5 <- formula(fishct ~ depth  + ranchn | depth)
f6 <- formula(fishct ~ depth  + depth2 | depth)

Zi1 <- zeroinfl(f1,dist="poisson", link="logit", data = dataforstats)
Zi2 <- zeroinfl(f2,dist="poisson", link="logit", data = dataforstats)
Zi3 <- zeroinfl(f3,dist="poisson", link="logit", data = dataforstats)
Zi4 <- zeroinfl(f4,dist="poisson", link="logit", data = dataforstats)
Zi5 <- zeroinfl(f5,dist="poisson", link="logit", data = dataforstats)

Ni5 <- zeroinfl(f1,dist="negbin", link="logit", data = dataforstats)
Ni6 <- zeroinfl(f2,dist="negbin", link="logit", data = dataforstats)
Ni7 <- zeroinfl(f3,dist="negbin", link="logit", data = dataforstats)
Ni8 <- zeroinfl(f4,dist="negbin", link="logit", data = dataforstats)
Ni9 <- zeroinfl(f5,dist="negbin", link="logit", data = dataforstats)
Ni10 <- zeroinfl(f6,dist="negbin", link="logit", data = dataforstats)

AICtab(Zi1,Zi2,Zi3,Zi4,Zi5,Ni5,Ni6,Ni7,Ni8,Ni9,Ni10,weights=TRUE,base=TRUE)
AICtab(Zi1,Zi3,Zi4,Zi5,Ni5,Ni7,Ni8,Ni9,Ni10,weights=TRUE,base=TRUE)

summary(Ni5)
confint(Ni5)
summary(Ni9)
confint(Ni9)
summary(Zi5)

model1 <- glm.nb(formula = f1, data = dataforstats, init.theta = 0.58, link = log)
summary(model1)

vuong(Ni5,model1)
AICtab(Ni5,model1,weights=TRUE,base = TRUE)

## Fig.4 ##
x<-seq(min(dataforstats$depth),max(dataforstats$depth),1)
name <-c(" Ald","Bir","Pal","Wil")
lr <- levels(dataforstats$ranchn)
par(mfrow=c (1,4))
for (k in 1:4) {
plot(dataforstats$depth,(dataforstats$fishct+1),type="n",log="y", ylim=c(1,100),xlim=c(1,60), main=name[k],xlab= "depth (cm)",ylab=("number of fish+1"),cex.lab=1.5,cex.main=1.7)
depth_dat <- dataforstats$depth[dataforstats$ranchn==lr[k]]
fish <- dataforstats$fishct[dataforstats$ranchn==lr[k]]
t <- table(fish,depth_dat)
dep <- unique(depth_dat)
pez <- unique(fish)
ord <- order(dep)
dep <- dep[ord]
orf <- order(pez)
pez <- pez[orf]
y11 <- y22 <- rep(0,length(x))
for (i in 1: length(dep)) {
  for (j in 1: length(pez)) {
  points(dep[i],(pez[j]+1),pch=1,cex=log(t[j,i])+1,col="blue")}}
for (i in 1:length(x)){
y11[i]=predict(Ni9,list(depth=x[i],ranchn=factor(k,levels=levels(dataforstats$ranchn))),type="response")}
lines(x,(y11+1),col="black",lwd=2)}

## Fig.5 ##
par(mfrow=c (1,3))
# Zero inflated model
Y_hats_Ni5 <- fitted(Ni5)
residuals_p <- residuals(Ni5, type="pearson")
plot(Y_hats_Ni5, residuals_p,pch=16,cex=0.5,main=" Pearson residuals of Model ")
boxplot(residuals_p~dataforstats$ranchn ,pch=16,cex=0.5,main=" Pearson residuals of Model ")
plot(dataforstats$depth, residuals_p,pch=16,cex=0.5,main=" Pearson residuals of Model ")

# Non-zero inflated model
Y_hats_model1 <- fitted(model1)
residuals_p_m <- residuals(model1, type="pearson")
plot(Y_hats_model1, residuals_p_m,pch=16,cex=0.5,main=" Pearson residuals of Model ")
boxplot(residuals_p_m~dataforstats$ranchn ,pch=16,cex=0.5,main=" Pearson residuals of Model ")
plot(dataforstats$depth, residuals_p_m,pch=16,cex=0.5,main=" Pearson residuals of Model ")

## Fig.6 ##
par(mfrow=c (1,2))
b <- seq(-1,8,1)
hist(residuals_p,breaks=b, main="zero inflated- negative binomial")
hist(residuals_p_m,breaks=b, main="negative binomial")

###Part II: with a Bayesian approach We do a simpler model only for ilustration
# We implement model f4

library(rjags)

n <- NROW(dataforstats)
y <- dataforstats$fishct
x <- dataforstats$depth
r1 <- rep(0,NROW(dataforstats))
r2 <- r3 <- r4 <- r1
r2[dataforstats$ranchn==2] <- 1
r3[dataforstats$ranchn==3] <- 1
r4[dataforstats$ranchn==4] <- 1

model =  "bayesian zero inflated poisson.R"

# Bundle data
# Bundle data
win.data=list(
  x=x,
  r2=r2,
  r3=r3,
  r4=r4,
  y=y,
  n=n
)

# Inits 
inits=list(
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
   beta4 = runif(1,-1,1)),
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
    beta4 = runif(1,-1,1)),
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
    beta4 = runif(1,-1,1)))

# Parameters to estimate
params <- c("alpha00","alpha11","beta00","beta1","beta2","beta3","beta4")

# MCMC settings
nc = 3 
ni=5000
nb=1000


# MCMC settings, start Gibbs sampler and plot results and diagnostics
jm=jags.model(model,data=win.data,inits=inits,n.chains=nc,n.adapt=nb)
update(jm, n.iter=ni)
zc=coda.samples(jm,variable.names=params,n.iter=1000)

# Display results
gelman.diag(zc)
plot(zc)
summary(zc)
Zi5

##################################################################

model =  "bayesian zero inflated neg binom.R"

# MCMC settings
nc = 3 
ni=100000
nb=50000


# Inits 
inits=list(
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
    beta4 = runif(1,-1,1), r = runif(1,1,3)),
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
    beta4 = runif(1,-1,1), r = runif(1,1,3)),
  list(
    alpha00 = runif(1,-1,1),alpha11 = runif(1,-1,1),  beta00 = runif(1,-1,1),
    beta1 = runif(1,-1,1), beta2 = runif(1,-1,1),beta3 = runif(1,-1,1),
    beta4 = runif(1,-1,1), r = runif(1,1,3)))

# Parameters to estimate
params <- c("alpha00","alpha11","beta00","beta1","beta2","beta3","beta4","r")

# MCMC settings, start Gibbs sampler and plot results and diagnostics
jm=jags.model(model,data=win.data,inits=inits,n.chains=nc,n.adapt=nb)
update(jm, n.iter=ni)
zcnb=coda.samples(jm,variable.names=params,n.iter=1000)

# Display results
gelman.diag(zcnb)
plot(zcnb)
summary(zcnb)
Ni9
summary(Ni9)