###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

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

precis(model1,digits=1,depth=2)
plot(model1)
pairs(model1)

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


model2 <- map2stan(
  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,
  constraints =list(theta ="lower=0"),
  start =list(a=1,b=1,c=1,scale=3),
  iter=4000,warmup=1000,chains =3
  )

precis(model2,digits=1,depth=2)
plot(model2)
pairs(model2)

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