###Part I: Preparing the data
rm(list=ls())
getwd()
library(rethinking)
library(rstan)
library(bbmle)
library(MASS)
##Set up your working directory manually or with the following functions##
#getwd()
#setwd("C:/Users/Pedro/Documents/Classes/Methods in Ecology II/Demos/")

## read and attach the COUNT data 
cd <- read.table("mantis.txt", header=T)
names(cd)

## calculate means
mean(cd$total[cd$type=="Total_Mantid"])
mean(cd$total[cd$type=="Total_Flower"])
mean(cd$total[cd$type=="zTotal_Control"])

## plot data
b = seq(0,35,1)
par( mfrow=c(1,3))
hist(cd$total[cd$type=="Total_Mantid"],breaks =b,xlab="visits", main="Mantis")
abline(v=8.121212, col="red")
hist(cd$total[cd$type=="Total_Flower"],breaks =b,xlab="visits", main="Flower")
abline(v=6.060606, col="red")
hist(cd$total[cd$type=="zTotal_Control"],breaks =b,xlab="visits", main="Control")
abline(v=0.4545455, col="red")

## models with different families
cd$tid <- 1
cd$tid[cd$type=="Total_Flower"] <- 2
cd$tid[cd$type=="zTotal_Control"] <- 3
cd$flo <- 1
cd$man <- 0
cd$man[cd$type=="Total_Mantid"] <- 1
cd$con <- 0
cd$con[cd$type=="zTotal_Control"] <- 1

# factor ID variable
cd$ID <- as.factor(cd$ID)

model1 <- ulam(
  alist(
    total ~ dpois(lambda),
    log(lambda) <- t[tid],
    t[tid] ~ dnorm(0,1)
  ),
  data = cd,chains =3, log_lik = TRUE
)

precis(model1,digits=1,depth=2)
par( mfrow=c(1,1))
plot(model1,pars=c("t[1]","t[2]","t[3]"), depth=2)
pairs(model1,horInd=c(1:3),verInd=c(1:3))

# plot sequence of generating chains
post1 <- extract.samples(model1)
par(mfrow=c(1,2))
for(i in names(post1)){
  dfpost1 <- data.frame(post1[i])
  plot(as.numeric(rownames(dfpost1)), dfpost1[,1], type = "l",
       xlab = "Index",
       ylab = i)
}


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


model2 <- ulam(
  alist(
    total ~ dgampois(pbar,scale),
    log(pbar) <- a + b*man + c*con,
    a ~ dnorm(0,10),
    b ~ dnorm(0,2),
    c ~ dnorm(0,2),
    scale ~ dcauchy(0,2)
  ),
  data = cd, log_lik = TRUE,
  constraints =list(theta ="lower=0"),
  iter=4000,warmup=1000,chains =3
)

precis(model2,digits=1,depth=2)
par(mfrow=c(1,1))
plot(model2)
pairs(model2,horInd=c(1:4),verInd=c(1:4))

# plot sequence of generating chains
post2 <- extract.samples(model2)
par(mfrow=c(1,5))
for(i in names(post2)){
  dfpost2 <- data.frame(post2[i])
  plot(as.numeric(rownames(dfpost2)), dfpost2[,1], type = "l",
       xlab = "Index",
       ylab = i)
}

compare(model1,model2)

model1 <- glm(total ~ type, family = poisson ,data=cd)
summary(model1)
model3 <- glm.nb(formula = total ~ type, init.theta = 0.1, link = log,data=cd)
summary(model3)

N <- nrow(cd)
E1 <- resid(model1,type="pearson")
p.pois <- length(coef(model1))
Dispersion.pois <- sum(E1^2)/(N-p.pois)
Dispersion.pois

E2 <- resid(model3,type="pearson")
p.nbin <- length(coef(model3))+1
Dispersion.nbin <- sum(E2^2)/(N-p.nbin)
Dispersion.nbin