--- title: "Ontogencit shift in isotopic niche space" author: "Volker Rudolf" date: "5/17/2019" output: html_document: df_print: paged reference_docx: word_style_draft1.docx word_document: null --- ```{r global_options, include=FALSE} knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE) ``` # Goal This document includes all the final analyses and figures for the stable isotope paper to examine ontogenetic niche shifts in natural communities. Briefly, main goal of this project is to obtain and compare statistics to desribe ontogenetic niche shifts based on stable isotopes ```{r, messsage=FALSE} #Load R Packages #list of packages that are needed pack<-c("tidyverse", "magrittr", "SIBER", "rjags","lubridate", "knitr","broom", "kableExtra", "cowplot", "RColorBrewer", "flextable","ggsci","wesanderson", "viridis") #check whether packages are installed and load them, if not get and then load package.check <- lapply(pack, FUN = function(x){ if (!require(x, character.only = TRUE)) { install.packages(x, dependencies = TRUE) library(x, character.only = TRUE) } }) #search() #verify that packages are loaded ``` ```{r, Rawdata} #get data sidat<-read_csv("StableIsotope_clean.csv") #note this will depend on where you store the file sidat<-sidat%>%filter(Date.Collected!="NA")%>%mutate(date=mdy(Date.Collected)) sidat<-sidat%>%mutate(Year=as.factor(year(date)), Month=month(date, label=TRUE)) sidat$Date<-as.Date(sidat$date) sidat$Species<-as.factor(sidat$Species) sidat$Location<-as.factor(sidat$Location) sidat$Species<-recode(sidat$Species,"Hyla.versicolor"="Hyla") sidat$Species<-recode(sidat$Species,"Rana.sphenocephala"="Rana") sidat$Species<-recode(sidat$Species,"Tiger.Beetle"="Pelocoris") ``` ```{r PlotTheme} singlefigtheme<-theme(axis.title.x=element_text(size=18, face="bold"), axis.title.y=element_text(size=18, face="bold"), axis.text.x=element_text(size=14), axis.text.y=element_text(size=14)) ``` ## Data overview ```{r SummaryStats} library(lubridate) datSum<-sidat%>%filter(!is.na(DryMass))%>% mutate(dateC=mdy(Date.Collected))%>% mutate(year=year(dateC), month=month(dateC))%>% group_by(Species,year, month,Location)%>% summarize(count=length(DryMass),meanM=mean(DryMass),minMass=min(DryMass),maxMass=max(DryMass)) autofit(flextable(datSum)) ``` Several analysis look at differences among stages and require distinct stages. Naturally, species differ in the number of developmental stages they have, and it's not always clear, i.e. some species have no distinct stages with continuous transition and some have big morphological jumps across very few stages (e.g. no stages in tadpoles, 11-13 stages in dragonflies, vs. 3 in beetles) which makes it tricky to make any comparisons between stages across species. Furthermore, sometimes the isotopes reflect the previous stage if they just did this transition. So need to make a comparison, we can simply use on common metric: mass. We assume that dry mass increases with age (which is usually true but there are exceptions), and that allows us to make create equal number of groups (=stages) across species to make comparisons. Not shown here, but overall, this works well for describing stages overall and typically captures some of the big jumps in some species. ### Data Preparation & Analysis ```{r clID} #Create function to create fixed number of clusters/groups for a given species clsize<-5 #define the cluster size #function to divide data in cluster and return assigned cluster ID for each individual clID<-function(x){ a<-kmeans(x,clsize) b<-a$cluster return(b) } ``` ```{r sidatG} #we apply the function above to bin species in "clusters" which represent stages based on the dry mass sidatG<-sidat%>%filter(!is.na(DryMass),Species!="Physid.Snail")%>% group_by(Species)%>%arrange(DryMass)%>% mutate(cluster=ntile(DryMass, clsize)) sidatG$Species<-as.factor(sidatG$Species) sidatG$stage<-as.factor(sidatG$cluster) sidatG1<-sidatG%>%mutate(TAX=fct_recode(Species,"Anura"="Rana","Anura"="Hyla","Anura"="Acris", "Nepidae"="Belostoma","Nepidae"= "Notonecta","Nepidae"="Pelocoris", "Nepidae"="Buenoa","Coleoptera"="Cybister", "Odonatea"="Tramea","Odonatea"="Plathemis","Odonatea"="Pachydiplax", "Odonatea"="Erythemis","Odonatea"="Libellula")) ``` This clustering works pretty well shown by this plot, where we sort individuals by dry mass see where groups "break" ```{r ClusterGroups} ggplot(subset(sidatG), aes(y=(DryMass),x=seq_along(DryMass),color=factor(cluster)))+ geom_point()+ ylab("Dry mass(mg)")+ facet_wrap(Location~Species, scale="free")+ theme_bw()+scale_x_discrete(name=NULL,labels=NULL)+ scale_color_viridis(discrete=T, name="stage") ``` Note, we can change the cluster size, increasing or decreasing it, too little means not enough resolution, too many and sample sizes get to small. Here 5 hits the sweetsop but results are robust to number of stages. # Stastical analyses ## TS: scaling coefficeint for shift in tropic height (d15N) First, we'll look at body size and d15N. Higher levels would indicate higher trophic position, and we are looking for change over ontogeny (here with size) ```{r, d15NAnalysis,message=FALSE} library(lme4) library(car) #look at data distribution #ggplot(sidatG,aes(x=d15N))+geom_density()+facet_wrap(Species~Location) #run glmm DMN<-lmer(d15N~DryMass*Species+(1|Location), data=sidatG) a<-(Anova(DMN)) Term<-rownames(a) a<-cbind(Term,a) autofit(flextable(a)) detach(package:car) ``` ```{r d15NFIgures, warning=FALSE} #let's extrac coefficients DMNs<-lmer(d15N~(DryMass):Species+(1|Location), data=subset(sidatG, cluster>0)) Ncoef<-tidy(DMNs, term="fixed") #note this requires broom package, might not work much longer if they change it NcoefS<-Ncoef%>%filter(str_detect(term,"DryMass")) NcoefS$term<-gsub("DryMass:Species","",NcoefS$term) NcoefS1<-ungroup(NcoefS)%>%mutate(Species=factor(term, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis","Libellula","Physid.Snail"))) NcoefS1$TAX<-c("Anura","Nepidae","Nepidae","Coleoptera","Odonata","Anura","Odonata","Nepidae","Odonata","Odonata","Anura","Nepidae","Odonata") #taxonomic group #lets plot the coefficients dNDryMSlopes<-ggplot(NcoefS1,aes(y=Species, x=estimate, color=TAX))+ geom_vline(aes(xintercept=0), color="black", size=1.5)+ geom_point(size=6)+ geom_errorbarh(aes(xmin=estimate-1.96*std.error, xmax=estimate+1.96*std.error),height=0, size=1)+ ylab("Species")+ xlab("Trophic scaling (TS)")+ xlim(xmin=-0.4,xmax=0.6)+ theme_classic()+ scale_color_manual(values=wes_palette(name="GrandBudapest1"), name=NULL) TSPlot<-dNDryMSlopes+ singlefigtheme+ theme(legend.position=c(0.8,0.85))+ geom_hline(yintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) TSPlot library(viridis) #and compare it to the empirical pattern dNDryMEmp<-ggplot(subset(sidatG, DryMass!="NA" ), aes(x=DryMass,y=d15N,color=factor(cluster), group=Location, shape=Location))+ labs(y=expression(delta^{15}*N), x=expression("Dry mass (mg)"))+ geom_point()+ scale_color_viridis(discrete=T, name="Stage")+ stat_smooth(se=T, method="lm", color="black")+ facet_wrap(~Species, scale="free")+theme_classic() dNDryMEmp ``` Overall, this shows that there is often a significant relationship between body size and trophic position, but it varies wildly by speciesand this relationship can be negative, neutral, and positive. Interestingly, several dragonfly predators show decline in trophic position with body size, which is bit surprising. ## DV: variation in diet composition across stages (d13C) Here we don't necessarily expact any gradual change of d13C with body size, because it's not simple enrichment with trophic position, and it's a mix of diets. To allow for more complicated patterns, instead of looking at body size, we examine differences in d13C across stages (clusters), treating them as factors (so no direction implied) and look at the variation explained by stages. ```{r Figd13CDryMas} #smoothed pattern carbon vs body size ggplot(subset(sidatG, DryMass!="NA" ), aes(x=DryMass,y=d13C,color=factor(cluster), group=Location, shape=Location))+ geom_point()+ stat_smooth(se=T, method="lm", color="black")+ scale_color_viridis(discrete=T, name="Stage")+ ylab(expression({delta}^13*C~'\u2030'))+ facet_wrap(~Species, scale="free")+theme_classic() ``` That's the raw data, so lets' do some stats and get predicted values ```{r d13CAnalyses, message=FALSE} library(car) ClustC<-lmer(d13C~stage*Species+(1|Location), data=sidatG) #plot(ClustC) CL<-(Anova(ClustC)) Term<-rownames(CL) CL<-cbind(Term,CL) autofit(flextable(CL)) detach(package:car) #we keep detaching it bc it interfers with some plyer functions #now lets extract explained variance issue: it's species specific, so we need sparate analysis for each species #that requires some double nesting (by Species and Location) and the glance function in broom package library(broom) DVnest<-sidatG%>%group_by(Species)%>%nest(.key=by_Species) DVnest2<-DVnest%>%mutate(by_Species=map(by_Species,~.x%>% group_by(Location)%>%nest(.key=by_Location))) DVM<-DVnest2%>%mutate(by_Species=map(by_Species,~.x %>% mutate(models=map(by_Location,~lm(d13C~stage,data=.x))))) DVAll<-DVM%>%mutate(by_Species=map(by_Species,~.x%>% mutate(tidied=map(models,glance)))) #put it all in table format, so we get species by location specific model output and r square for stage effect DV<-DVAll%>%unnest%>%unnest(tidied)%>%arrange(r.squared) ``` ```{r DVFIG} DV1<-ungroup(DV)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis","Libellula","Physid.Snail"))) DV1<-DV1%>%mutate(TAX=fct_recode(Species,"Anura"="Rana","Anura"="Hyla","Anura"="Acris", "Nepidae"="Belostoma","Nepidae"= "Notonecta","Nepidae"="Pelocoris", "Nepidae"="Buenoa","Coleoptera"="Cybister", "Odonatea"="Tramea","Odonatea"="Plathemis","Odonatea"="Pachydiplax", "Odonatea"="Erythemis","Odonatea"="Libellula")) DVFIG<-ggplot(DV1,aes(y=Species, x=r.squared, color=TAX, shape=Location))+ geom_point(size=6,alpha=0.8 )+ xlab("Diet variation across stages (DV)")+ ylab("Species")+ theme_classic()+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE) DVPlot<-DVFIG+singlefigtheme+theme(legend.position=c(0.85,0.85)) + geom_hline(yintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) DVPlot ``` And predicted plot ```{r d13CFigures} library(effects) e<-as.data.frame(effect("stage:Species",ClustC)) d13CFicModel<-ggplot(e,aes(x=stage, y=fit))+ geom_point(size=3)+ geom_errorbar(aes(ymin=fit-se,ymax=fit+se), width=0)+ theme_bw()+ #geom_line(linetype="dashed")+ #geom_point(data=sidatG,aes(x=cluster,y=d13C),size=0.5)+ ylab(expression("Predicted"~{delta}^13*C~'\u2030'~"(mean, ± 1se)"))+ facet_wrap(~Species) d13CFicRaw<-ggplot(sidatG,aes(x=stage, y=d13C, shape=Location))+ geom_point(size=0.5)+ stat_summary(fun.data=mean_se, size=.5)+ ylab(expression({delta}^13*C~'\u2030'))+ xlab("Stage")+ theme_bw()+ facet_wrap(~Species) d13CFicModel d13CFicRaw ``` Note that predicted values tend differ from observed estimates from species because predicted values account for the location effect. If we remove location effect, they are spot on. But since location is significant, we do want to account for it. So overall, some species show a significant change in diet indicated by differences in d13C. This shift varies across Species, hence significant interaction effect. Predicted effects include 95% CI, showing differences among stages. # Trophic niche space So, let's integrate both N and C to look at the niche space occupied by a stage. To do that, we need to work with stages/clusters from now on so we can define distinct stages to get their niche space. We can assume that the trophic niche is given by the d13C and d15N coordination space, and the trophic niche of a stage of a given species is given by the individuals within that stage. One of the preferred and unbiased metrics (e.g. compared to hull) is the standard ellipse area (SEA) that includes x% of the population. A common value is 40%. The advantagne of this metric is that it is not sensitive to outliers, and pretty robust to sample size differences (although we have very similar sample sizes within a given species, there is some variation across species). Furthemore, we can use the ellipse space to calucate key metrics to look at niche shifts. Lets start with a basic ML based plot to show the standard ellipse for all species, with 40% boundary ```{r EllipsePlot} col<-brewer.pal(n=(clsize+4),name="Blues") C_N_EllFig<-sidatG%>%select(ID,Species,Location,Year,Month,d13C,d15N,stage)%>% ggplot(aes(y=d15N,x=d13C,color=stage, shape=Location, linetype=Location))+ geom_point(alpha=0.5)+ stat_ellipse(type="norm",level=0.40,size=1.2)+ labs(x = expression({delta}^13*C~'\u2030'), y= expression({delta}^15*N~'\u2030'))+ facet_wrap(~Species, scale="free")+ theme_bw()+ scale_color_viridis(discrete=T, name="Stage") #scale_color_manual(values=c(col[1:clsize+4]))+ theme(legend.position=c(0.7,0.075),legend.direction="horizontal", legend.margin=margin(-0.5,0,0,0,"cm")) C_N_EllFig ``` Next, we need to decide of whether to go with ML or bayesian approach. Overall, it seems that ML approach makes comparisons simpler across species, but it's more sensitive to sample size. But sample size corrected SEAc is supposed to be basically the same as SAEbayesian, and recent simulation suggest correction isn't totally worht it (Sycaranta et al 2013).SEAc=SEA*(n-1)/(n-2), where n is the sample size minimum 3, better 10. NOte that my sample sizes here are usually >30 per stage/cluster, one of the perks of clustering and lots of sampling. ```{r SEAcstats} #function to get ellipse area for each stage/cluster SEAcstats<-function(x) { Metrics<-matrix(NA,ncol=clsize, nrow=1) colnames(Metrics)<-c(paste0("SEAc_stage", as.character(1:clsize))) #now, we need to drop the species and ID and make a SIBER object #select the right columns of the data specSI<-as.data.frame(x%>%ungroup%>%select(iso1,iso2,group,community)) specsiber<-createSiberObject(specSI) #now we are ready to caluclate some SIBER stats stageN<-max(specSI$group) #determine number of stages SIstats<-data.frame(Area1=integer(),Area2=integer(),overlap=integer(),stage=integer()) com<-specSI$community[1] MLE<-groupMetricsML(specsiber) Metrics[1:clsize]<-as.vector(MLE[3,]) #this pulls out the small sample size corrected metric: SEAc Metrics<-as.data.frame(Metrics) return(Metrics) } ``` ```{r Eoverlap} #function to get ellipse overlap for all sttage combinatino within a given species Eoverlap<-function (x) { Metrics<-matrix(NA,ncol=10, nrow=1) #get data in siber format specSI<-as.data.frame(x%>%ungroup%>%select(iso1,iso2,group,community)) #select the right columns of the data specsiber<-createSiberObject(specSI) #now we are ready to caluclate some SIBER stats #ellipse overlap stageN<-max(specSI$group) #determine number of stages #create data SIstats<-data.frame(Area1=integer(),Area2=integer(),overlap=integer(),stage=integer()) com<-specSI$community[1] count=0 for (k in 1:(stageN-1)) { for(j in 2:stageN) { count=count+1 el1<-as.character(com+k/10) el2<-as.character(com+j/10) overlapB<-maxLikOverlap(el1,el2,specsiber,p.interval=0.95,n=100)[1:3] #proportional overlap relative to non-overlapping area overlapP<-overlapB[3]/(overlapB[1]+overlapB[2]-overlapB[3]) SIstats[count,1:2]<-as.vector(overlapB)[1:2] SIstats[count,3]<-overlapP #SIstats[count,1:3]<-as.vector(maxLikOverlap(el1,el2,specsiber,p.interval=0.95,n=100))[1:3] SIstats[count,4]<-paste(k,j) } } SIstatsOv<-SIstats%>%select(overlap, stage)%>% subset(stage!="2 2" & stage!="3 3" & stage!="4 4" & stage!="5 5" & stage!="6 6" & stage!="7 7" & stage!="3 2" & stage!="4 2" & stage!="4 3"& stage!="5 2" & stage!="5 3" & stage!="5 4")%>% spread(stage,overlap) Metrics<-as.data.frame(SIstatsOv) #Metrics[[1]]<-mean(SIstats$overlap) #Metrics[[2]]<-max(SIstats$overlap) #Metrics[[3]]<-min(SIstats$overlap) #Metrics<-as.data.frame(Metrics) return(Metrics) } ``` ```{r sidatBSI, message=FALSE} #Now that we have those functions, let's get some statistics #First get things in SIBER format #rename and extract variables we want sidatBSI<-sidatG%>%arrange(Species,Location,cluster)%>% select(ID,Species,Location,Year,Month,d13C,d15N,cluster)%>% rename(iso1=d13C,iso2=d15N, group=cluster) SITableSum<-sidatBSI%>%group_by(Species,group)%>%summarize(N=length(iso1))%>%spread(group,N) #autofit(flextable(SITableSum)%>%add_header_row(value=c("","stage"),colwidths=c(2,5))) #Then apply function to each species in the data set for a given location #makes new variable indicating the community and returns them as numbers Community<-sidatG%>%arrange(Species,Location,cluster)%>% transmute(ID,community=as.numeric(recode(Location,"DC108_1"="1","Huntsville"="2","NicksPond"="3"))) #add community to main file sidatBSI<-right_join(sidatBSI,Community) #and get it in long format for plotting sumSidat<-sidatBSI%>%filter(!is.na(iso1))%>% group_by(Species,group,community)%>% summarize(observations=length(iso1))%>% arrange(Species,community, group) autofit(flextable(sumSidat)) ``` From this, we see a few things: (1) Tramea has bare minimum samples size for some small stages in location 1 (DC108-1 and we could probably drop it) (2) Tiger bettles plathemis libelulla & hyla some of lowest sample sizes, with <=11 per stage (3) sample sizes are nice and even across stages for most parts ## SEA (for stages) ```{r SEAData, warning=FALSE,message=FALSE} SEAstat<-sidatBSI%>%select(-c(Location))%>% group_by(Species,community)%>%do(SE=SEAcstats(.)) SEAF<-unnest(SEAstat) #flextable(SEAF) ``` ```{r SEAPlot} #convert in long format to plot range<-clsize+2 SEAFL<-SEAF%>%gather("stage", "SEA",3:range) #ggplot(subset(SEAFL, Species!="Libellula"), aes(x=stage,y=SEA,color=factor(community)))+ # geom_point()+ # facet_wrap(~Species, scale="free")+ # theme_bw() SEAFL1<-ungroup(SEAFL)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis","Libellula","Physid.Snail"))) SEAFL1<-SEAFL1%>%mutate(TAX=fct_recode(Species,"Anura"="Rana","Anura"="Hyla","Anura"="Acris", "Nepidae"="Belostoma","Nepidae"= "Notonecta","Nepidae"="Pelocoris", "Nepidae"="Buenoa","Coleoptera"="Cybister", "Odonatea"="Tramea","Odonatea"="Plathemis","Odonatea"="Pachydiplax", "Odonatea"="Erythemis","Odonatea"="Libellula")) meanSEAsp<-SEAFL1%>%group_by(Species, TAX, community)%>%summarize(mean.SEA=mean(SEA),sdSEA=sd(SEA)/length(SEA)^0.5) SEACP<-ggplot(meanSEAsp,aes(x=mean.SEA, y=Species, color=TAX, shape=factor(community)))+ geom_errorbarh(aes(xmin=mean.SEA-sdSEA,xmax=mean.SEA+sdSEA), height=0, size=2)+ geom_point(size=6)+ geom_point(data=SEAFL1,aes(x=SEA, y=Species, color=TAX,shape=factor(community)),size=2.5, alpha=0.5)+ scale_shape_manual(values=c(17,16,18), guide=FALSE)+ xlab("Stage specific trophic niche are (SEAc)")+ ylab("Species")+ theme_classic()+theme(legend.position="none")+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL) SEACPlot<-SEACP+singlefigtheme+theme(legend.position=c(0.85,0.85)) + geom_hline(yintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) SEACPlot ``` Test for differences in SEAc across stages &Species ```{r SEATests, message=FALSE} library(car) SEAFL$community<-as.factor(SEAFL$community) SEAModel<-lmer(log(SEA)~stage*Species+(1|community),data=subset(SEAFL, Species!="Libellula")) #plot(SEAModel) SEAT<-as.data.frame(Anova(SEAModel)) Terms<-rownames(SEAT) SEAT<-cbind(Terms,SEAT) autofit(flextable(SEAT)) detach(package:car) ``` ## SEA overlap (P statistic) ```{r EOV, message=FALSE, warning=FALSE} Ovstat<-sidatBSI%>%select(-c(Location))%>% group_by(Species,community)%>% do(Ov=Eoverlap(.)) EOV<-unnest(Ovstat) saveRDS(EOV,file="EOV.rds") ``` ```{r OverlapPlot} EOV<-readRDS("EOV.rds") #flextable(EOV) #to plot we need to get it in long-format Crange<-clsize+2 EOVL<-EOV%>%gather("stages","overlap",3:Crange) EOVL$stages<-gsub(" ","_",EOVL$stages) EOVL1<-ungroup(EOVL)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis","Libellula","Physid.Snail"))) EOVL1<-EOVL1%>%mutate(TAX=fct_recode(Species,"Anura"="Rana","Anura"="Hyla","Anura"="Acris", "Nepidae"="Belostoma","Nepidae"= "Notonecta","Nepidae"="Pelocoris", "Nepidae"="Buenoa","Coleoptera"="Cybister", "Odonatea"="Tramea","Odonatea"="Plathemis","Odonatea"="Pachydiplax", "Odonatea"="Erythemis","Odonatea"="Libellula")) EOVLM<-EOVL1%>%group_by(Species, TAX,community)%>% summarize(meanO=mean(overlap), sdO=sd(overlap)/length(overlap)^0.5) EOVP<-ggplot(EOVLM,aes(x=meanO, y=Species, color=TAX, shape=factor(community)))+ geom_errorbarh(aes(xmin=meanO-sdO,xmax=meanO+sdO), height=0, size=2)+ geom_point(size=6)+ geom_point(data=EOVL1,aes(x=overlap, y=Species, color=TAX,shape=factor(community)),alpha=0.5, size=2.5)+ scale_shape_manual(values=c(17,16,18), guide=FALSE)+ xlab("Niche overlap between stages (P)")+ ylab("Species")+ theme_classic()+theme(legend.position="none")+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL) EOVPlot<-EOVP+singlefigtheme+theme(legend.position=c(0.9,0.85)) + geom_hline(yintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) EOVPlot ``` Overall, this shows (1) there is never 100% overlap (2) there is considerable variability in overlap among stages (3) and species We can look at whether stages that are closer to gether are more similar ```{r OverlapByStage} EOVLMS<-EOVL%>%group_by(stages)%>%summarize(meanSO=mean(overlap)) ggplot(EOVL,aes(x=stages, y=overlap))+ geom_point()+ stat_summary(fun.data=mean_se, size=1.5)+ ylab("SEA overlap (%)")+ xlab("Stage comparison")+ theme_classic()+singlefigtheme ``` ```{r OverlapTrendAnalysis} EOVLD<-EOVL%>%separate(stages,c("stageA","stageB"))%>% mutate_if(is.character,as.numeric)%>% mutate(stageDif=stageB-stageA) ovtrend<-lmer(overlap~stageDif*Species+(1|community), data=EOVLD) #summary(ovtrend) #plot(ovtrend) library(car) EOVT<-as.data.frame(Anova(ovtrend)) Terms<-rownames(EOVT) EOVT<-cbind(Terms,EOVT) autofit(flextable(EOVT)) detach(package:car) ``` ```{r OVerlapTrendFig} ggplot(EOVLD,aes(x=stageDif, y=overlap))+ geom_point(colour="gray")+ stat_summary(fun.data=mean_se, size=1.5)+ ylab("SEA overlap (%)")+ xlab("Difference in developmental stage")+ theme_classic()+singlefigtheme ggplot(EOVLD,aes(x=stageDif, y=overlap))+ geom_point(colour="gray")+ stat_smooth(method=lm, color="black")+ facet_wrap(~Species)+ ylab("SEA overlap (%)")+ xlab("Difference in developmental stage")+ theme_classic()+singlefigtheme ``` It does look like overlap decreases with stage differences, but this doesn't have to be the case. Not surprising given what we know about dN and dC patterns. Some changed, others didn't much. ALso note that in some species early stages look most different, (e.g. Acris), in others, last stage (e.g. Buenoa, Notonecta) ## HA: ontogenetic niche width We need to get a general idea for the how much nice space changes across all stages and how this differs across species in an integrated way. This is done by looking at the hull area given by the centroids of each stage-specific ellipse. ```{r HullArea} #get centroid Cen<-sidatBSI%>%group_by(Species,Location,group)%>%summarise(iso1=mean(iso1),iso2=mean(iso2)) hull<-Cen%>%group_by(Species,Location)%>%slice(chull(iso1,iso2)) #create the hull from the data ggplot(subset(Cen),aes(x=iso1,y=iso2, color=Location, group=Location))+ geom_point()+ theme_classic()+ geom_polygon(data=hull, alpha=0.2)+ scale_color_viridis(discrete=T)+ facet_wrap(~Species, scale="fixed") #lets make a little function to get hull area from the lymanMetrics in siber package TotalHull<-function(x){ dat<-as.data.frame(x%>%ungroup%>%select(iso1,iso2)) hull<-laymanMetrics(dat$iso1,dat$iso2)$metrics["TA"] return(hull) } #use function to get hull area for each species by location combination HullA<-unnest(Cen%>%group_by(Species,Location)%>%do(hullA=TotalHull(.))) #flextable(HullA) ``` ```{r HullAPlot, warning=FALSE,message=FALSE} #and plot it species_av <- HullA %>% summarize(avg = mean(hullA, na.rm = T)) %>%pull(avg) HullA1<-ungroup(HullA)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis","Libellula","Physid.Snail"))) HullA1<-HullA1%>%mutate(TAX=fct_recode(Species,"Anura"="Rana","Anura"="Hyla","Anura"="Acris", "Nepidae"="Belostoma","Nepidae"= "Notonecta","Nepidae"="Pelocoris", "Nepidae"="Buenoa","Coleoptera"="Cybister", "Odonatea"="Tramea","Odonatea"="Plathemis","Odonatea"="Pachydiplax", "Odonatea"="Erythemis","Odonatea"="Libellula")) HullAP<-ggplot(HullA1,aes(x=hullA,y=Species,color=TAX, shape=Location))+ geom_point(size=6, alpha=0.9)+ ylab("Species")+ xlab("Ontogenetic diversity (Hull area)")+ theme_classic()+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE) legendtheme<-theme(legend.position=c(0.9,0.9), legend.title=element_text(size=14,face="bold"), legend.text=element_text(size=12)) HullAPlot<-HullAP+singlefigtheme+legendtheme+ geom_hline(yintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) HullAPlot ggsave("HullArea.png",HullAPlot,height=5,width =7, dpi=600) ggsave("HullArea.pdf",HullAPlot,height=5,width =7, dpi=600) ``` ## Link between traits and ontogentic niche shifts Now that we have all these statistics, lets' put some species level metrics together and see how position in this multidimensional niche space is linked to species traits First get species traits ```{r TraitT} #first, we need species level trait metrics TraitT<-sidatG%>%group_by(Species,Location)%>% summarise(DMav=mean(DryMass),DMmin=quantile(DryMass,probs=0.05),DMmax=quantile(DryMass,probs=0.95))%>%mutate(DMR=DMmax-DMmin) TraitT$FG<-c("Herb","Pred","Pred","Pred","Pred","Pred","Pred","Pred","Herb","Pred","Pred","Pred","Pred","Herb","Pred","Pred","Pred") #addotion of feeding type TraitT$FT<-c("Scraping","sucking","sucking","sucking","sucking","chewing","chewing","chewing","Scraping","chewing","sucking","chewing","chewing","Scraping","sucking","chewing","chewing") #feeding mode TraitT$TAX<-c("Anura","Nepidae","Nepidae","Nepidae","Coleoptera","Odonata","Odonata","Odonata","Anura","Odonata","Nepidae","Odonata","Odonata","Anura","Nepidae","Odonata","Odonata") #taxonomic group TraitT<-TraitT%>%select(TAX, everything()) autofit(flextable(TraitT)) ``` Plot size range ```{r sizeplot} TraitT1<-ungroup(TraitT)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis", "Libellula","Physid.Snail"))) TraitT1$TAX<-factor(TraitT1$TAX,levels=c("Anura","Nepidae","Coleoptera","Odonata")) P1<-ggplot(TraitT1, aes(Species, color=TAX, shape=Location))+ coord_flip()+ geom_linerange(aes(ymin=log(DMmin),ymax=log(DMmax), color=TAX), size=1.5, position=position_dodge(width=0.6))+ geom_point(aes(Species,log(DMav)),size=6,position=position_dodge(width=0.6))+ ylab("log(Dry Mass) [mg]")+xlab("Species")+ singlefigtheme+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE) SizeFig<-P1+theme(legend.position=c(0.75,0.9),legend.text=element_text(size=12))+ geom_vline(xintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) SizeFig ggsave("SizeFig.png",SizeFig,height=5,width =7, dpi=600) ``` Then we get ontogenetic niche shift metrics ```{r SPecM, message=FALSE, warning=FALSE} #now lets get ontogentic niche shift metrics for all species #average ellipse area SEAav<-SEAFL%>%mutate(Location=(recode(community,"1"="DC108_1","2"="Huntsville","3"="NicksPond")))%>% group_by(Species,Location)%>%summarise(meanSEAc=mean(SEA)) #add hull area SPecM<-left_join(HullA,SEAav) #add average niche overlap SEAovStat<-EOVLM%>% mutate(Location=(recode(community,"1"="DC108_1","2"="Huntsville","3"="NicksPond")))%>% select(Species,Location,meanO,sdO) SPecM<-left_join(SEAovStat,SPecM) #add diet variation explained DVS<-DV%>%select(Species,Location, r.squared)%>%rename(CVar=r.squared) SPecM<-left_join(SPecM,DVS) #add trophic scaling TSA<-NcoefS%>%rename(Species=term, TS=estimate)%>%select(Species,TS) SPecM<-merge(SPecM,TSA) #note, we don't have a location specific scaling, so a Species always gets the same model estimate #get trophic height range TR from d15N based on differences in centroids dNR1<-sidatG%>%filter(cluster==1 | cluster==clsize)%>% group_by(Species,Location,cluster)%>% summarise(d15NM=mean(d15N))%>% summarize(TR=max(d15NM)-min(d15NM))%>% select(Species,Location, TR) SPecM<-left_join(SPecM,dNR1) #get diet range DR from d13C centroids dCR<-sidatG%>%group_by(Species,Location,stage)%>%summarise(md13C=mean(d13C))%>% summarise(d13Cmin=min(md13C),d13Cmax=max(md13C))%>% mutate(DR=d13Cmax-d13Cmin)%>%select(Species,Location,DR) SPecM<-left_join(SPecM,dCR) #Make plots SPecM1<<-ungroup(SPecM)%>%mutate(Species=factor(Species, levels=c("Rana","Hyla","Acris","Belostoma", "Notonecta","Pelocoris","Buenoa","Cybister", "Tramea","Plathemis","Pachydiplax","Erythemis", "Libellula","Physid.Snail"))) TRP<-ggplot(SPecM1,aes(x=Species,y=TR, color=TAX, shape=Location))+ geom_point(size=6, alpha=0.8)+ xlab("Species")+ ylab("Trophic range (TR)")+ coord_flip()+ theme_classic()+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE) legendtheme<-theme(legend.position=c(0.9,0.9), legend.title=element_text(size=14,face="bold"), legend.text=element_text(size=12)) TRPlot<-TRP+singlefigtheme+legendtheme+ geom_vline(xintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) TRPlot DRP<-ggplot(SPecM1,aes(x=Species,y=DR, color=TAX, shape=Location))+ geom_point(size=6, alpha=0.8)+ xlab("Species")+ ylab("Diet range (DR)")+ coord_flip()+ theme_classic()+ scale_color_manual(values=wes_palette(name="GrandBudapest1"),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE) legendtheme<-theme(legend.position=c(0.9,0.9), legend.title=element_text(size=14,face="bold"), legend.text=element_text(size=12)) DRPlot<-DRP+singlefigtheme+legendtheme+ geom_vline(xintercept=c(3.5,8.5,7.5),linetype=2)+ guides(color = guide_legend(reverse = TRUE)) DRPlot ``` ```{r} #Print the table autofit(flextable(SPecM)) ``` #and lets look real quick at correlations among traits ```{r TraitCorr} library(corrplot) SPecM1<-SPecM%>%select(-sdO) SPecM1<-SPecM1%>%rename(P=meanO, HA=hullA, SEAc=meanSEAc,DV=CVar) cormat<-round(cor(SPecM1[,4:10]),3) corrplot(cormat, type="upper", order="hclust", col=brewer.pal(n=8,name="PuOr"), addCoef.col="black", diag=F) cor.mtest <- function(mat, ...) { mat <- as.matrix(mat) n <- ncol(mat) p.mat<- matrix(NA, n, n) diag(p.mat) <- 0 for (i in 1:(n - 1)) { for (j in (i + 1):n) { tmp <- cor.test(mat[, i], mat[, j], ...) p.mat[i, j] <- p.mat[j, i] <- tmp$p.value } } colnames(p.mat) <- rownames(p.mat) <- colnames(mat) p.mat } # matrix of the p-value of the correlation p.mat <- cor.mtest(SPecM1[,4:10]) head(p.mat[, 1:7]) corrplot(cormat,diag=F, type="upper", order="hclust", p.mat=p.mat,sig.level=0.05, method="ellipse", col=brewer.pal(n=8,name="PuOr")) ``` Now we can (1) graph species in multidimensional space based on traits and niche metrics ```{r OrdiSpecNichFig, message=FALSE,warning=FALSE} library(vegan) #first, need to make some ordination. Best if we standardize metrics, since all are on totally different scales SPecMO<-SPecM%>%select(Species,TAX,Location,everything()) NormSPecM<-decostand(SPecMO[,4:11],"range") #sets range from 0-1 Species<-SPecM$Species Location<-SPecM$Location DistSpecM<-metaMDS(NormSPecM, distance="euclidean", trace=FALSE) scors<-as.data.frame(scores(DistSpecM, display="sites")) scors<-cbind(scors,Species) scors<-left_join(scors,TraitT) vf<-envfit(DistSpecM,NormSPecM) #this fits the Trait names to data frame #look at how ordination is related to species traits (used to work, need to fix) # vT<-envfit(DistSpecM,TraitT[,3:9]) #plot species in ontogenetic traits niche space spp.scrs <- as.data.frame(scores(vf, display = "vectors")) spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs)) #get hull hull2<-subset(scors)%>%group_by(TAX)%>%slice(chull(NMDS1,NMDS2)) TraitOrdiP<-ggplot(scors)+ geom_point(mapping = aes(x = NMDS1, y = NMDS2, color=Species,shape=FT), size=5)+ #stat_ellipse(aes(x=NMDS1,y=NMDS2, group=TraitT$FT, linetype=FT),type="norm",level=0.5)+ coord_fixed()+ geom_segment(data = spp.scrs,aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),arrow = arrow(length = unit(0.25, "cm")), colour = "grey")+ geom_text(data = spp.scrs, aes(x = NMDS1, y = NMDS2, label = Species),size = 3) #TraitOrdiP #vector loading spp.scrs$Otrait<-c("P","sdP","HA", "SEAc","DV","TS","TR","DR") VEc<-ggplot(spp.scrs)+ coord_fixed()+ geom_segment(data = subset(spp.scrs,Otrait!="sdP"), aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "black", size=1)+ geom_text(data = subset(spp.scrs,subset=Otrait %in% c("DV","DR","TR")), aes(x = NMDS1, y = NMDS2, label = Otrait),size = 3.5,hjust=1,vjust=-0.2)+ geom_text(data = subset(spp.scrs,Otrait=="P"), aes(x = NMDS1, y = NMDS2, label = Otrait),size = 3.5,hjust=-0.6,vjust=0)+ geom_text(data = subset(spp.scrs,subset=Otrait %in% c("SEAc","HA","TS")), aes(x = NMDS1, y = NMDS2, label = Otrait),size = 3.5,hjust=1,vjust=1)+ ylab("NMDS2")+xlab("NMDS1")+ coord_cartesian(xlim=c(-1.4,0.8),ylim=c(-1.1,1.1)) Vec<-VEc+theme_classic(base_size=8)+theme(axis.title=element_blank(), axis.text=element_blank()) #cybister groups with dragnflies, not with hemiptera. library(ggrepel) hull2$TAX<-factor(hull2$TAX,levels=c("Anura","Nepidae","Coleoptera","Odonata")) scors$TAX<-factor(scors$TAX,levels=c("Odonata","Coleoptera","Nepidae","Anura")) hullP<-ggplot(spp.scrs,aes(x=NMDS1,y=NMDS2))+ geom_polygon(data=hull2, aes(fill=factor(TAX)), alpha=0.4)+ geom_point(data=scors,aes(x=NMDS1,NMDS2,color=TAX, shape=Location),size=2)+ geom_text(data=scors,aes(x=NMDS1,NMDS2, label=Species, color=TAX),vjust=1,size=3)+ coord_cartesian(xlim=c(-1.3,0.8), ylim=c(-1,1.1))+ scale_color_manual(values=rev(wes_palette(name="GrandBudapest1")),name=NULL)+ scale_fill_manual(values=rev(wes_palette(name="GrandBudapest1")),name=NULL)+ scale_shape_manual(values=c(17,16,18), guide=FALSE)+ guides(color=FALSE) hullF<-hullP+coord_fixed(xlim=c(-1.4,0.8),ylim=c(-1.1,1.1),ratio=1)+ theme(legend.position=c(0.02,0.9),legend.text=element_text(size=12))+ singlefigtheme+ draw_plot(Vec, x = .1, y = .5, width = .9, height = .7) hullF ggsave("HullFigure.png",hullF,height=5,width =5, dpi=600) save_plot("HullFigure.png",hullF,base_aspect_ratio=1,base_height=6) ``` ```{r, warning=FALSE} SPecM2<-as.tibble(SPecM) SPecM2$Species<-as.factor(SPecM2$Species) TraitT$Species<-as.factor(TraitT$Species) AllTable<-left_join(TraitT,SPecM2, by="Species") AllTable<-AllTable%>%select(-TAX.y,-Location.y,-FT,-FG) #flextable(AllTable) cormatT<-round(cor(AllTable[,4:15]),3) flextable(as.data.frame(cormatT)) #corrplot(cormatT) ```