### Method to identify disturbance events in palynological sequences for rangel transitional forest (Figueroa Rangel et al. 2012) #in R: "obervations"-> rows, "variables"->columns rangeltrans <- read.csv("FigueRangelTransitF2012_red.csv", sep=",") #read in dataset with taxa columns in numeric format #### (removed 4 rows from original file: analysisUnitName, Thickness, Sample name, radiocarbon dates, ) str(rangeltrans) # 1. step: I need to subsample the dataset in a way that I can process with ordinations (remove all non-count parts of data # and remove all taxa that are not within the pollen sum (for this I need to keep the "group" column and taxa column)) rangeltrans2<-rangeltrans[-c(1:2),-c(3:5)] #set all na values to zero values in the counts rangeltrans2[is.na(rangeltrans2)]<-0 str(rangeltrans2) rangeltrans2[,c(3:ncol(rangeltrans2))]<-lapply(rangeltrans2[,c(3:ncol(rangeltrans2))],as.numeric) #for when columns are read in as integers str(rangeltrans2) #to calculate percentages, I need the whole terrestrial pollen dataset #change group names: until here on 4.12. rangeltrans_terr<-with(rangeltrans2,rangeltrans2[rangeltrans2$group == "CTRSH" | rangeltrans2$group == "CH" | rangeltrans2$group == "I" | rangeltrans2$group == "F" | rangeltrans2$group == "TH" | rangeltrans2$group == "TF" | rangeltrans2$group == "GF" | rangeltrans2$group == "GH" | rangeltrans2$group == "SH" | rangeltrans2$group == "SHT",]) #Cloud trees and shrubs, cloud herbs, plus ferns because they were included in sum, Indeterminate, Pine forest herbs, Pine dominated forest trees, #Gallery forest, Gallery forest herbs, Secondary habitat herbs, Secondary habitat trees # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) rangeltrans_terr_num<-rangeltrans_terr[,-c(1:2)] rangeltranspcts = sapply(names(rangeltrans_terr_num), function(x) { (rangeltrans_terr[x] / sum(rangeltrans_terr_num[x]))*100 }) rangeltranspcts<-as.data.frame(rangeltranspcts) rangeltranspcts[is.na(rangeltranspcts)]<-0 #I need to reassign sample names, use column (variable-) names from rangeltrans2 names(rangeltranspcts)<-names(rangeltrans2[,c(3:ncol(rangeltrans2))]) #For the clustering, I will use all the taxa in the terrestrial pollen sum, but I will also need to sum the percentages of the taxa found within #the vegetation type of interest #Here I only want the arboreal taxa from the Andean forest (and mangroves) #(this step comes later after having calculated percentages on whole terrestrial pollen sum) #I need to subsample according to groups again, so I need to reattach group column to dataset rangeltrans_pcts_names<-cbind(rangeltrans_terr[,1:2],rangeltranspcts) rangeltrans_arb_pcts<-with(rangeltrans_pcts_names,rangeltrans_pcts_names[rangeltrans_pcts_names$group == "CTRSH" | rangeltrans_pcts_names$group == "TF",]) #cloud forest trees and Pine-dominated forest trees rangeltrans_arb_pcts_num<-rangeltrans_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(rangeltranspcts) df<-as.data.frame(df) #add species names as columns names #get species names: rangeltranspcts_ch<-as.character(rangeltrans_pcts_names$ï..Rangel_Transition) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(rangeltranspcts_ch) # 3. Step: Transform variables (sqrt most probably, check) df_sqr<-sqrt(df) # 4. Step: Check if linear or unimodal ordination methods should be used: library(vegan) ord<-decorana(df_sqr) #does not work if there are no values different from zero in one row plot(ord) summary(ord) #length DCA axis 1: 1.8832 -> linear with PCA #Plot the diagram to identify clusters on diagram #For plotting also a summary curve of arobreal, etc. percentages that were included in ordination sumperc_ord<-apply(rangeltrans_arb_pcts_num,2,sum) ### CLUSTERING PER ZONE TO CALCULATE MEAN AND SD PER ZONE library(rioja) diss<-dist(df_sqr) clust<-chclust(diss, method="coniss") rangeltrans_depth<-rangeltrans[1,c(6:ncol(rangeltrans))] trangeltrans_depth<-t(rangeltrans_depth) df_depth<-cbind(trangeltrans_depth,sumperc_ord,df) #windows(width=14, height=10) rangeltrans_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangeltrans_plot,clust,nZone = 3, col="red") rangeltrans_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangeltrans_plot_short,clust,nZone = 3, col="red") #__________________________________________________________________________ #same but not divided by zones: I run this here because zones were non-significant and compare it to the result with the zones mean_rangeltrans<-mean(sumperc_ord) sd_rangeltrans<-sd(sumperc_ord) lower_th<-mean_rangeltrans-sd_rangeltrans df_lower_th<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th,]) df_lower_th [,1:2] #T11 88 25.32468,T12 96 27.66798,T13 104 24.84277,T14 112 29.11392 #___________________________________________________________________________ df_depth[,1:2] rangeltrans_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangeltrans_plot_short,clust,nZone = 3, col="red") addZone(rangeltrans_plot_short,upper=87, lower=113,col=rgb(1,0,0,alpha=0.3))#disturbance 1 #not as disturbance event, start of disturbance not clear/only one sample in "valley of disturbance" df_depth[,1:2] #Choose disturbance events where more than one sample is below the local threshold within the same "valley of percentages" (is this justifiable?) # -> # disturbance events: 1 #age can be seen in dataset "rangeltrans" #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec #see excel fiel RecoveryRates_BiolLetters for calculations