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