rm(list=ls())
getwd()
#setwd("C:/Users/Pedro/Documents/Classes/Methods in Ecology II/Demos/")
par(mfrow=c (1,1))
dataforstats <- read.table("dataforstats_final_BB121613.txt", header=T)
names(dataforstats)

# install glmmADMB
# library(devtools)
# install_github("bbolker/R2admb",subdir="R2admb")
# devtools::install_github("bbolker/glmmadmb")
library(glmmADMB)
library(bbmle)
library(mgcv)
library(plot3D)

sub1 <-subset(dataforstats,!is.na(upland_elev_m ))
summary(sub1)

par(mfrow=c (1,2))
hist(sub1$macroct,20,main="Macroinvertebrates",xlab="Number of organisms")
hist(sub1$fishct,20,main="Fish",xlab="Number of organisms")

y<-as.numeric(sub1$macroct)
depth <- sub1$depth
depth2 <- depth^2 
fish <- log(sub1$fishct+1)
fish2 <- fish^2
rancho <- rep(1,length(sub1$ranch))
rancho[sub1$ranch=="bir"] <- 2
rancho[sub1$ranch=="pal"] <- 3
rancho[sub1$ranch=="wil"] <- 4
rancho <- factor(rancho)

par(mfrow=c (2,2),mar=c(3,4,1,2))
scatter3D((fish[rancho==1]+1),log(depth[rancho==1]),y[rancho==1],theta = 40, phi = 30,
          xlab="fish",ylab="depth",zlab="macroinvertebrates",main="ald")
scatter3D((fish[rancho==2]+1),log(depth[rancho==2]),y[rancho==2],theta = 40, phi = 30,
          xlab="fish",ylab="depth",zlab="macroinvertebrates",main="bir")
scatter3D((fish[rancho==3]+1),log(depth[rancho==3]),y[rancho==3],theta = 40, phi = 30,
          xlab="fish",ylab="depth",zlab="macroinvertebrates",main="pal")
scatter3D((fish[rancho==4]+1),log(depth[rancho==4]),y[rancho==4],theta = 40, phi = 30,
          xlab="fish",ylab="depth",zlab="macroinvertebrates",main="wil")

M1 <- lm(y ~ depth + depth2 + rancho + fish)
summary(M1)

M2<-glmmadmb(y~ depth + depth2  + fish + fish2 + rancho, data = sub1, zeroInflation = FALSE,family="nbinom")
summary(M2)


#M3 <- gam(y ~ s(depth, fx=FALSE, k=-1, bs="cr")+rancho,family=negbin(c(.01,1)))
#summary(M3)
#gam.check(M3)

M4 <- gam(y ~ s(depth, fx=FALSE, k=-1, bs="cr")+ s(fish,fx=FALSE,k=-1,bs="cr")+rancho,family=negbin(c(0.669)))
summary(M4)


AIC(M1,M2,M4)

##########################################FISHCT#########################################
#macroct
title <- levels(sub1$ranch)
x<-seq(min(depth),max(depth),1)
#################################################
par(mfrow=c (1,4),mar=c(3,4,3,2))
for (k in 1:4) {
  plot(depth,y+1,type="n",ylim=c(.01,max(sub1$macroct)),xlim=c(-1,max(depth)), main=title[k],xlab= "depth (cm)", ylab=("number of moacroinvertebrates"))
  ############################################
  ord <- order(depth[rancho==k & sub1$fishct ==0])
  depth_dat <- depth[rancho==k & sub1$fishct ==0]
  depth_dat <-depth_dat[ord]
  macro1 <- y[rancho==k & sub1$fishct ==0]
  macro1 <- macro1[ord]
  
  t <- table(macro1,depth_dat)
  dep <- unique(depth_dat)
  mac <- unique(macro1)
  ord <- order(dep)
  dep <- dep[ord]
  orm <- order(mac)
  mac <- mac[orm]
  
  for (i in 1: length(dep)) {
    for (j in 1: length(mac)) {
      points(dep[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="blue")}}
  
  ############################################
  ord <- order(depth[rancho==k & sub1$fishct >0])
  depth_dat <- depth[rancho==k & sub1$fishct >0]
  depth_dat <-depth_dat[ord]
  mac1 <- y[rancho==k & sub1$fishct >0]
  mac1 <- mac1[ord]
  
  t <- table(mac1,depth_dat)
  dep <- unique(depth_dat)
  mac <- unique(mac1)
  ord <- order(dep)
  dep <- dep[ord]
  orm <- order(mac)
  mac <- mac[orm]
  
  for (i in 1: length(dep)) {
    for (j in 1: length(mac)) {
      points(dep[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="black")}}
  ############################################
  y2=predict(M2,list(depth=x,depth2=x^2,fish=rep(5,length(x)),fish2=rep(5^2,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y2),col="red",lwd=2)
  y3=predict(M1,list(depth=x,depth2=x^2,fish=rep(5,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y3),col="green",lwd=2)
  #y4=predict(M4,list(depth=x,fish=rep(0,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  #lines(x,(y4),col="blue",lwd=2)
  y5=predict(M4,list(depth=x,fish=rep(5,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y5),col="black",lwd=2)
}

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


par(mfrow=c(3,3))
res1 <-residuals(M1)
res2 <-residuals(M2)
#res3 <-residuals(M3)
res4 <-residuals(M4)

plot((sub1$macroct)+1,res1, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates")
plot((sub1$macroct)+1,res2, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates")
#plot((sub1$macroct)+1,res3, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates")
plot((sub1$macroct)+1,res4, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates")

plot(fish+1,res1, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish")
plot(fish+1,res2, ylim =c(-4,18), log="x",ylab="residual", xlab="fish")
#plot(fish+1,res3, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish")
plot(fish+1,res4, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish")

plot(rancho,res1,ylim =c(-3,18),ylab="residuals", xlab="ranches")
plot(rancho,res2,ylim =c(-3,18),ylab="residuals", xlab="ranches")
#plot(rancho,res3,ylim =c(-3,18),ylab="residuals", xlab="ranches")
plot(rancho,res4,ylim =c(-3,18),ylab="residuals", xlab="ranches")

par(mfrow=c(1,3))
y_predicted1 <- predict(M1,type="response")
y_predicted2 <- predict(M2,type="response")
#y_predicted3 <- predict(M3,type="response")
y_predicted4 <- predict(M4,type="response")
tb1 <-table(sub1$macroct,round(y_predicted1,0))
tb2 <-table(sub1$macroct,round(y_predicted2,0))
#tb3 <-table(sub1$macroct,round(y_predicted3,0))
tb4 <-table(sub1$macroct,round(y_predicted4,0))
fishy <- unique(sub1$macroct)
o1 <- order(fishy)
fishy <-fishy[o1]
pfishy1 <- unique(round(y_predicted1,0))
pfishy2 <- unique(round(y_predicted2,0))
#pfishy3 <- unique(round(y_predicted3,0))
pfishy4 <- unique(round(y_predicted4,0))
o2 <- order(pfishy1)
o3 <- order(pfishy2)
#o4 <- order(pfishy3)
o5 <- order(pfishy4)
pfishy1 <- pfishy1[o2]
pfishy2 <- pfishy2[o3]
#pfishy2 <- pfishy2[o4]
pfishy3 <- pfishy2[o5]
plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20))
for (i in 1:dim(tb1)[1]) {
  for (j in 1:dim(tb1)[2]) {
    if (i==1) v <- tb1[i,j]/30+1
    if (i>1) v <-log(tb1[i,j])+0.5
    points(fishy[i],pfishy1[j]-1,pch=1,cex=(v),col="black")
  }}
plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20))
for (i in 1:dim(tb2)[1]) {
  for (j in 1:dim(tb2)[2]) {
    if (i==1) v <- tb2[i,j]/30+1
    if (i>1) v <-log(tb2[i,j])+0.5
    points(fishy[i],pfishy2[j]-1,pch=1,cex=(v),col="black")
  }}
#plot(fishy,fishy,type="n",ylim=c(1,6),xlab="observed",ylab="predicted", xlim=c(-5,20))
#for (i in 1:dim(tb3)[1]) {
#  for (j in 1:dim(tb3)[2]) {
#  if (i==1) v <- tb3[i,j]/30+1
#  if (i>1) v <-log(tb3[i,j])+0.5
#     points(fishy[i],pfishy3[j],pch=20,cex=(v),col="black")
#                          }}
plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20))
for (i in 1:dim(tb4)[1]) {
  for (j in 1:dim(tb4)[2]) {
    if (i==1) v <- tb4[i,j]/30+1
    if (i>1) v <-log(tb4[i,j])+0.5
    points(fishy[i],pfishy4[j]-1,pch=1,cex=(v),col="black")
  }}

##########################################FISHCT#########################################
#macroct & fish
title <- levels(sub1$ranch)
x<-seq(min(fish),max(fish),0.1)
#################################################
par(mfrow=c (1,4))
for (k in 1:4) {
  plot(fish,y+1,type="n",ylim=c(-1,max(sub1$macroct)),xlim=c(-1,max(fish)), main=title[k],xlab= "log(fish+1)", ylab=("number of moacroinvertebrates"))
  ############################################
  ord <- order(fish[rancho==k])
  fish_dat <- fish[rancho==k]
  fish_dat <-fish_dat[ord]
  mac <- y[rancho==k]
  mac <- mac[ord]
  
  t <- table(mac,fish_dat)
  fisha <- unique(fish_dat)
  mac <- unique(mac)
  ord <- order(fisha)
  fisha <- fisha[ord]
  orc <- order(mac)
  mac <- mac[orc]
  
  for (i in 1: length(fisha)) {
    for (j in 1: length(mac)) {
      points(fisha[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="blue")}}
  
  
  ############################################
  y2=predict(M2,list(depth=rep(10,length(x)),depth2=rep(10^2,length(x)),fish=x,fish2 =x^2,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y2),col="red",lwd=2)
  y3=predict(M1,list(depth=rep(10,length(x)),depth2=rep(10^2,length(x)),fish=x,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y3),col="green",lwd=2)
  #y4=predict(M4,list(depth=x,fish=rep(0,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  #lines(x,(y4),col="blue")
  y5=predict(M4,list(depth=rep(10,length(x)),fish=x,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response")
  lines(x,(y5),col="black",lwd=2)
}

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

par(mfrow=c(2,2),mar=c(4,4,4,4))

q<-seq(min(depth),41,1)
y1<-y2<-y3<-y4 <- array(0,c(40,40))
for (i in 1:length(x)){
  for (j in 1:length(q)){
    y1[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("1",levels=levels(rancho))),type="response")
    y2[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("2",levels=levels(rancho))),type="response")
    y3[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("3",levels=levels(rancho))),type="response")
    y4[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("4",levels=levels(rancho))),type="response")
  }}

image(q,x,y1,xlab="depth",ylab="fish",main="ald",zlim=c(0,4), col=gray((0:10)/10))
image(q,x,y2,xlab="depth",ylab="fish",main="bir",zlim=c(0,4), col=gray((0:10)/10))
image(q,x,y3,xlab="depth",ylab="fish",main="pal",zlim=c(0,4), col=gray((0:10)/10))
image(q,x,y4,xlab="depth",ylab="fish",main="wil",zlim=c(0,4), col=gray((0:10)/10))  

image(q,x,y1,xlab="depth",ylab="fish",main="ald",zlim=c(0,4), col=hcl.colors(10, "Greens"))
image(q,x,y2,xlab="depth",ylab="fish",main="bir",zlim=c(0,4), col=hcl.colors(10, "Greens"))
image(q,x,y3,xlab="depth",ylab="fish",main="pal",zlim=c(0,4), col=hcl.colors(10, "Greens"))
image(q,x,y4,xlab="depth",ylab="fish",main="wil",zlim=c(0,4), col=hcl.colors(10, "Greens")) 

q2 <- c(rep(q,5))
x2 <- c(rep(x[1],40),rep(x[2],40),rep(x[3],40),rep(x[4],40),rep(x[5],40))
y1a <- y1[1:200]
length(x2)
length(q2)
length(y1a)
par(mfrow=c(1,1))
scatter3D(q2,x2,y1a)