### Method to identify disturbance events in palynological sequences for La Cocha (Gonzalez Carranza et al. 2012) #in R: "obervations"-> rows, "variables"->columns gonzcarranza <- read.csv("GonzCarranza2012_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(gonzcarranza) # 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)) gonzcarranza2<-gonzcarranza[-c(1:2),-c(3:5)] #set all na values to zero values in the counts gonzcarranza2[is.na(gonzcarranza2)]<-0 str(gonzcarranza2) #gonzcarranza2[,c(3:ncol(gonzcarranza2))]<-lapply(gonzcarranza2[,c(3:ncol(gonzcarranza2))],as.numeric) #for when columns are read in as integers #str(gonzcarranza2) #to calculate percentages, I need the whole terrestrial pollen dataset #change group names: until here on 4.12. gonzcarranza_terr<-with(gonzcarranza2,gonzcarranza2[gonzcarranza2$group == "A" | gonzcarranza2$group == "B" | gonzcarranza2$group == "C" | gonzcarranza2$group == "D" | gonzcarranza2$group == "F",]) #A: Subandean forest, B: Andean forest, C: subparamo, D: paramo, F: Dry ind. # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) gonzcarranza_terr_num<-gonzcarranza_terr[,-c(1:2)] gonzcarranzapcts = sapply(names(gonzcarranza_terr_num), function(x) { (gonzcarranza_terr[x] / sum(gonzcarranza_terr_num[x]))*100 }) gonzcarranzapcts<-as.data.frame(gonzcarranzapcts) gonzcarranzapcts[is.na(gonzcarranzapcts)]<-0 #I need to reassign sample names, use column (variable-) names from gonzcarranza2 names(gonzcarranzapcts)<-names(gonzcarranza2[,c(3:ncol(gonzcarranza2))]) #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 gonzcarranza_pcts_names<-cbind(gonzcarranza_terr[,1:2],gonzcarranzapcts) gonzcarranza_arb_pcts<-with(gonzcarranza_pcts_names,gonzcarranza_pcts_names[gonzcarranza_pcts_names$group == "A" | gonzcarranza2$group == "B",]) #Subandean and andean forest gonzcarranza_arb_pcts_num<-gonzcarranza_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(gonzcarranzapcts) df<-as.data.frame(df) #add species names as columns names #get species names: gonzcarranzapcts_ch<-as.character(gonzcarranza_pcts_names$ï..ZaireCocha) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(gonzcarranzapcts_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.6653 -> 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(gonzcarranza_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") bstick(clust, ng=10, plot = TRUE) #6 zones gonzcarranza_depth<-gonzcarranza[1,c(6:ncol(gonzcarranza))] tgonzcarranza_depth<-t(gonzcarranza_depth) df_depth<-cbind(tgonzcarranza_depth,sumperc_ord,df) #windows(width=14, height=10) gonzcarranza_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(gonzcarranza_plot,clust,nZone = 6, col="red") gonzcarranza_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(gonzcarranza_plot_short,clust,nZone = 6, col="red") plot(clust) (x <- identify(clust)) #click on leaves of dendogramm groups and when finished click esc, results will print in console zone1<-x[[1]] zone2<-x[[2]] zone3<-x[[3]] zone4<-x[[4]] zone5<-x[[5]] zone6<-x[[6]] #calculate means and sd per zone #Zone 1 rows mean_gonzcarranza1<-mean(sumperc_ord[1:c(length(zone1))]) sd_gonzcarranza1<-sd(sumperc_ord[1:c(length(zone1))]) lower_th1<-mean_gonzcarranza1-sd_gonzcarranza1 df_lower_th1<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th1,]) df_lower_th1[,1:3] #LC_1 169,LC_26 221,LC_27 223,LC_29 227,LC_31 231 zone1 #Zone2 rows mean_gonzcarranza2<-mean(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) sd_gonzcarranza2<-sd(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) lower_th2<-mean_gonzcarranza2-sd_gonzcarranza2 df_lower_th2<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th2,]) df_lower_th2[,1:3] # LC_42 257,LC_45 271,LC_68 407,LC_69 409,LC_71 979,LC_72 981 zone2 #Zone3 rows mean_gonzcarranza3<-mean(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) sd_gonzcarranza3<-sd(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) lower_th3<-mean_gonzcarranza3-sd_gonzcarranza3 df_lower_th3<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th3,]) df_lower_th3[,1:3] #LC_74 987 zone3 #Zone4 rows mean_gonzcarranza4<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) sd_gonzcarranza4<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) lower_th4<-mean_gonzcarranza4-sd_gonzcarranza4 df_lower_th4<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th4,]) df_lower_th4[,1:3] #LC_85 1013,LC_86 1017 zone4 #Zone5 rows mean_gonzcarranza5<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5))]) sd_gonzcarranza5<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5))]) lower_th5<-mean_gonzcarranza5-sd_gonzcarranza5 df_lower_th5<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th5,]) df_lower_th5[,1:3] #LC_93 1031,LC_95 1035,LC_104 1079 zone5 #Zone6 rows mean_gonzcarranza6<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5)+length(zone6))]) sd_gonzcarranza6<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4)+length(zone5)+length(zone6))]) lower_th6<-mean_gonzcarranza6-sd_gonzcarranza6 df_lower_th6<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th6,]) df_lower_th6[,1:3] #LC_112 1095,LC_114 1099,LC_127 1125,LC_128 1127,LC_130 1131 zone6 #__________________________________________________________________________ #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_gonzcarranza<-mean(sumperc_ord) sd_gonzcarranza<-sd(sumperc_ord) lower_th<-mean_gonzcarranza-sd_gonzcarranza df_lower_th<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th,]) df_lower_th [,1:3] #LC_1 169,LC_26 221,LC_31 231,LC_80 1001,LC_82 1007,LC_85 1013,LC_86 1017,LC_112 1095,LC_114 1099, #LC_127 1125,LC_128 1127,LC_130 1131 #___________________________________________________________________________ ##LC_1 169,LC_26 221,LC_27 223,LC_29 227,LC_31 231 # LC_42 257,LC_45 271,LC_68 407,LC_69 409,LC_71 979,LC_72 981 #LC_74 987 #LC_85 1013,LC_86 1017 #LC_93 1031,LC_95 1035,LC_104 1079 #LC_112 1095,LC_114 1099,LC_127 1125,LC_128 1127,LC_130 1131 df_depth[,1:3] gonzcarranza_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(gonzcarranza_plot_short,clust,nZone = 3, col="red") addZone(gonzcarranza_plot_short,upper=220, lower=224,col=rgb(1,0,0,alpha=0.3))#first disturbance event addZone(gonzcarranza_plot_short,upper=226, lower=228,col=rgb(1,0,0,alpha=0.3)) addZone(gonzcarranza_plot_short,upper=230, lower=232,col=rgb(1,0,0,alpha=0.3)) addZone(gonzcarranza_plot_short,upper=1000, lower=1002,col=rgb(1,0,0,alpha=0.3)) addZone(gonzcarranza_plot_short,upper=1006, lower=1018,col=rgb(1,0,0,alpha=0.3)) addZone(gonzcarranza_plot_short,upper=1093, lower=1100,col=rgb(1,0,0,alpha=0.3)) addZone(gonzcarranza_plot_short,upper=1124, lower=1132,col=rgb(1,0,0,alpha=0.3))#start of disturbance not clear #not as disturbance event, start of disturbance not clear/only one sample in "valley of disturbance" df_depth[,1:3] #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 "gonzcarranza" #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec #see excel fiel RecoveryRates_BiolLetters for calculations