---
title: "Ants Seeds SIBER"
author: "Israel Del Toro"
date: "July, 2018"
output: 
  word_document :
   toc: true
---
##Prepare data and load necessary packages for plotting and GLMM
```{r}
#May the Code be with You#
setwd('~/Dropbox/Seed Dispersal/Data')

library(vegan)
library(ggplot2)
library (lme4)
library(lmerTest)
library(gridExtra)
source("summary_se.R")
library (MASS)
library(car)
library("FactoMineR")
library("devtools")
library("ggbiplot")

#Overdispersion test function- Modified from B. Bolker July 2018.
overdisp_fun <- function(model) {
    rdf <- df.residual(model)
    rp <- residuals(model,type="pearson")
    Pearson.chisq <- sum(rp^2)
    prat <- Pearson.chisq/rdf
    pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
    c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

#load and subset datasets----
data_part<- read.csv('seed_depot_data_10_04_16.csv', header=T)
data_part$Elevation<-as.factor(data_part$Elevation)
data_part$seeds.removed<-as.numeric (data_part$seeds.removed)
str(data_part)

#Partition the dataset by species
datura<-data_part[which(data_part$Seed=="Datura"),]
iris<-data_part[which(data_part$Seed=="Iris"),]
oat<-data_part[which(data_part$Seed=="Oat"),]
sumac<-data_part[which(data_part$Seed=="Sumac"),]
```

#Use lme4 to create GLMM and test for overdispersion 
```{r}
#GLMM function try with a poisson distribution for count data. Global model 
data_part$id <- with(data_part, paste(Site, Elevation, Depot.., sep = "_"))
m1<- glmer(seeds.removed ~ hours + Elevation*Seed+ (1|Site) + (1+hours |id), family= poisson, data_part)
summary(m1)
Anova(m1)

#Overdispersion test
overdisp_fun(m1)
#These results indicate the global model is not overdispersed. 

###Figure 1 Global model of total number of seeds removed over time.
M1<- summarySE(data_part, measurevar="seeds.removed", groupvars=c("hours","Elevation"))
total.seeds.graph<-ggplot(M1, aes(x=hours , y=seeds.removed, group = Elevation)) +
  geom_errorbar (aes(ymin=seeds.removed-se, ymax=seeds.removed+se), size = 1, width=.3, position=position_dodge(.01)) + theme_bw() +
  geom_line(aes(), size = 1) +
  geom_point(aes(shape = Elevation), size = 6) +
  theme(legend.position=("none")) + theme(text = element_text(size=20)) +    
  xlab ("Time") +ylab("Seeds removed") 
total.seeds.graph

tiff (filename = "Figure1.tiff", height= 4, width= 4, units= "in", res=300)
total.seeds.graph
dev.off()

###ANOVA Table
datura$id <- with(datura, paste(Site, Elevation, Depot.., sep = "_"))
Datura.lmer<- glmer (seeds.removed ~ hours + Elevation + Site + (1|Site) + (1+hours|id), family= poisson, datura)
summary(Datura.lmer)
Anova(Datura.lmer)
overdisp_fun(Datura.lmer)

iris$id <- with(iris, paste(Site, Elevation, Depot.., sep = "_"))
Iris.lmer<- glmer (seeds.removed ~ hours + Elevation + Site + (1|Site) + (1+hours|id), family= poisson, iris)
summary(Iris.lmer)
Anova(Iris.lmer)
#This results suggest an overdispersion effect on the Iris LMER We can still fit is assuming a pseudo Poisson distribution, 

overdisp_fun(Iris.lmer)
#This model is over dispersed so must be corrected to interpret the potential source of Type 1 error (below):
## extract summary table; you may also be able to do this via
##  broom::tidy or broom.mixed::tidy
cc <- coef(summary(Iris.lmer))
phi <- overdisp_fun(Iris.lmer)["ratio"]
cc <- within(as.data.frame(cc),
{   `Std. Error` <- `Std. Error`*sqrt(phi)
    `z value` <- Estimate/`Std. Error`
    `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE)
})
printCoefmat(cc,digits=3)
#The estimated pseudo-poisson fit suggest that the effect of Iris seen at High elevations is largely driven by the trends at two mountains LNF and MOG 

oat$id <- with(oat, paste(Site, Elevation, Depot.., sep = "_"))
Oat.lmer<- glmer (seeds.removed ~ hours + Elevation + Site + (1|Site) + (1+hours|id), family= poisson, oat)
summary(Oat.lmer)
Anova(Oat.lmer)
overdisp_fun(Oat.lmer)

sumac$id <- with(sumac, paste(Site, Elevation, Depot.., sep = "_"))
Sumac.lmer<- glmer (seeds.removed ~ hours + Elevation + Site + (1|Site) + (1+hours|id), family= poisson, sumac)
summary(Sumac.lmer)
Anova(Sumac.lmer)
overdisp_fun(Sumac.lmer)
#This model is over dispersed so must be corrected to interpret the potential source of Type 1 error (below):
## extract summary table; you may also be able to do this via
##  broom::tidy or broom.mixed::tidy
cc <- coef(summary(Sumac.lmer))
phi <- overdisp_fun(Sumac.lmer)["ratio"]
cc <- within(as.data.frame(cc),
{   `Std. Error` <- `Std. Error`*sqrt(phi)
    `z value` <- Estimate/`Std. Error`
    `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE)
})
printCoefmat(cc,digits=3)

####Figure 2 in the paper (4 panel with each species)
Datura<- summarySE(datura, measurevar="seeds.removed", groupvars=c("hours","Elevation"))
head (Datura)
datura.graph<-ggplot(Datura, aes(x=hours , y=seeds.removed, group = Elevation)) +
  geom_errorbar (aes(ymin=seeds.removed-se, ymax=seeds.removed+se), size = 1, width=.3, position=position_dodge(.01)) + theme_bw() +
  geom_line(aes(), size = 1) +
  geom_point(aes(shape = Elevation, colour = Elevation), size = 6) +
  theme(legend.position=("none")) + theme(text = element_text(size=20)) +    
  xlab ("Time") +ylab("Seeds removed") +
  annotate("text", x = 22, y = 22, label = "Datura")
datura.graph

Iris<- summarySE(iris, measurevar="seeds.removed", groupvars=c("hours","Elevation"))
Iris
iris.graph<-ggplot(Iris, aes(x=hours , y=seeds.removed, group = Elevation)) +
  geom_errorbar (aes(ymin=seeds.removed-se, ymax=seeds.removed+se), size = 1, width=.3, position=position_dodge(.01)) + theme_bw() +
  geom_line(aes(), size = 1) +
  geom_point(aes(shape = Elevation, colour = Elevation), size = 6) +
  theme(legend.position=("none")) + theme(text = element_text(size=20)) +    
  xlab ("Time") +ylab("Seeds removed")  +
  annotate("text", x = 22, y = 22, label = "Iris")
iris.graph

Oat<- summarySE(oat, measurevar="seeds.removed", groupvars=c("hours","Elevation"))
Oat
oat.graph<-ggplot(Oat, aes(x=hours , y=seeds.removed, group = Elevation)) +
  geom_errorbar (aes(ymin=seeds.removed-se, ymax=seeds.removed+se), size = 1, width=.3, position=position_dodge(.01)) + theme_bw() +
  geom_line(aes(), size = 1) +
  geom_point(aes(shape = Elevation, colour = Elevation), size = 6) +
  theme(legend.position=("none")) + theme(text = element_text(size=20)) +    
  xlab ("Time") +ylab("Seeds removed")  +
  annotate("text", x = 22, y = 22, label = "Oat")
oat.graph

Sumac<- summarySE(sumac, measurevar="seeds.removed", groupvars=c("hours","Elevation"))
Sumac
sumac.graph<-ggplot(Sumac, aes(x=hours , y=seeds.removed, group = Elevation)) +
  geom_errorbar (aes(ymin=seeds.removed-se, ymax=seeds.removed+se), size = 1, width=.3, position=position_dodge(.01)) + theme_bw() +
  geom_line(aes(), size = 1) +
  geom_point(aes(shape = Elevation), size = 6) +
  theme(legend.position=c(.8, 0.8)) + theme(text = element_text(size=20)) +  
  xlab ("Time") +ylab("Seeds removed")  +
  annotate("text", x = 22, y = 22, label = "Sumac")
sumac.graph

tiff
grid.arrange(datura.graph, iris.graph, oat.graph, sumac.graph, ncol=2)
dev.off()

df<-read.csv("~/Dropbox/Seed Dispersal/Data/Genera Abundance.csv", header=TRUE)
genera.env<-(df[, c(1:2)])
ant_abund<-df[, c(1,2,18)]
genera<-df[, c(3:17)]
str(genera)

```
# Plot therelationship between abundance and number of seeds removed 
```{r, echo=FALSE}
#relationship between abundance and number of seeds removed 
whisker1<-ggplot(ant_abund, aes(Elevation, Abundance)) + 
  geom_boxplot() + theme_bw() +
  scale_y_continuous("Abundance (# of of traps with ants)")
plot(whisker1)
```

#Principal component analyses of ant genera 
```{r}
cor(genera)
pca.data= genera
genera.pca = PCA(genera, scale.unit=TRUE, ncp=5, graph=T)

#PCA for coordinate plots 
pca <- prcomp(genera, scale. = TRUE)

#PCA with ellipses for Elevation and points for sites. 
pca.total.elev <- ggbiplot(pca, obs.scale = 1, var.scale = 1, group=df$Elevation, ellipse = T) +
  scale_colour_manual(values = c("dark green", "blue", "orange")) + scale_shape_manual(values = c(15, 16, 17, 18))+ 
  geom_point(size = 3, aes(colour=df$Elevation,shape=df$Elevation)) + theme_bw() + theme(panel.grid.minor=element_blank(), panel.grid.major=element_blank()) + 
  theme(legend.position="none") + theme(text=element_text(size = 18, colour="black"))
plot(pca.total.elev)
```