### Method to identify disturbance events in palynological sequences for rangelpine (Figueroa Rangel et al. 2008) #in R: "obervations"-> rows, "variables"->columns rangelpine <- read.csv("FigueRangelPineF2008_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(rangelpine) # 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)) rangelpine2<-rangelpine[-c(1:2),-c(3:5)] #set all na values to zero values in the counts rangelpine2[is.na(rangelpine2)]<-0 str(rangelpine2) rangelpine2[,c(3:ncol(rangelpine2))]<-lapply(rangelpine2[,c(3:ncol(rangelpine2))],as.numeric) #for when columns are read in as integers str(rangelpine2) #to calculate percentages, I need the whole terrestrial pollen dataset #change group names: until here on 4.12. rangelpine_terr<-with(rangelpine2,rangelpine2[rangelpine2$group == "CTRSH" | rangelpine2$group == "CH" | rangelpine2$group == "I" | rangelpine2$group == "F" | rangelpine2$group == "PH" | rangelpine2$group == "POF" | rangelpine2$group == "GF" | rangelpine2$group == "GH" | rangelpine2$group == "SH" | rangelpine2$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) rangelpine_terr_num<-rangelpine_terr[,-c(1:2)] rangelpinepcts = sapply(names(rangelpine_terr_num), function(x) { (rangelpine_terr[x] / sum(rangelpine_terr_num[x]))*100 }) rangelpinepcts<-as.data.frame(rangelpinepcts) rangelpinepcts[is.na(rangelpinepcts)]<-0 #I need to reassign sample names, use column (variable-) names from rangelpine2 names(rangelpinepcts)<-names(rangelpine2[,c(3:ncol(rangelpine2))]) #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 rangelpine_pcts_names<-cbind(rangelpine_terr[,1:2],rangelpinepcts) rangelpine_arb_pcts<-with(rangelpine_pcts_names,rangelpine_pcts_names[rangelpine_pcts_names$group == "CTRSH" | rangelpine_pcts_names$group == "POF",]) #cloud forest trees and Pine-dominated forest trees rangelpine_arb_pcts_num<-rangelpine_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(rangelpinepcts) df<-as.data.frame(df) #add species names as columns names #get species names: rangelpinepcts_ch<-as.character(rangelpine_pcts_names$ï..Rangel_Pine) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(rangelpinepcts_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.8766 -> 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(rangelpine_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") rangelpine_depth<-rangelpine[1,c(6:ncol(rangelpine))] trangelpine_depth<-t(rangelpine_depth) df_depth<-cbind(trangelpine_depth,sumperc_ord,df) #windows(width=14, height=10) rangelpine_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangelpine_plot,clust,nZone = 3, col="red") rangelpine_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangelpine_plot_short,clust,nZone = 3, col="red") #__________________________________________________________________________ mean_rangelpine<-mean(sumperc_ord) sd_rangelpine<-sd(sumperc_ord) lower_th<-mean_rangelpine-sd_rangelpine df_lower_th<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th,]) df_lower_th [,1:2] #C4 8 53.08642,C12 24 52.03252,C14 28 53.73563,C20 40 49.05660,C36 72 37.85166,C40 80 52.27273,C42 84 45.39474 #___________________________________________________________________________ #P2 4 62.87263,P4 8 67.12565 #P23 46 73.26733,P24 48 70.17893,P27 54 71.08434,P35 70 74.13793 #P37 74 47.88136 df_depth[,1:2] rangelpine_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(rangelpine_plot_short,clust,nZone = 3, col="red") addZone(rangelpine_plot_short,upper=3, lower=5,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, start of disturbance not clear addZone(rangelpine_plot_short,upper=7, lower=9,col=rgb(1,0,0,alpha=0.3))#not as disturbance event,only one sample in "valley of disturbance" addZone(rangelpine_plot_short,upper=45, lower=49,col=rgb(1,0,0,alpha=0.3))#disturbance 1 addZone(rangelpine_plot_short,upper=53, lower=55,col=rgb(1,0,0,alpha=0.3))#disturbance 2 addZone(rangelpine_plot_short,upper=69, lower=71,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, start of disturbance not clear addZone(rangelpine_plot_short,upper=73, lower=74,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, start of disturbance not clear #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 "rangelpine" #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec #see excel fiel RecoveryRates_BiolLetters for calculations