#### #### #### 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  + depth2 + 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)

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,Ni5,Ni6,Ni7,Ni8,Ni9,Ni10,weights=TRUE,base=TRUE)
AICtab(Zi1,Zi3,Zi4,Ni5,Ni7,Ni8,Ni9,Ni10,weights=TRUE,base=TRUE)

summary(Ni5)
confint(Ni5)
summary(Ni7)

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(Ni5,list(depth=x[i],depth2=x[i]^2,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")