---
title: "Ants Seeds SIBER"
author: "Israel Del Toro"
date: "March 19, 2018"
output: word_document
---

```{r cars}
#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")

#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"),]

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

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

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)

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)

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)

####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("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)

#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)")
```

```{r whisker1, echo=FALSE}
plot(whisker1)
```

PCA
```{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"))
```


```{r PCA, echo=FALSE}
plot(pca.total.elev)
```

