### Method to identify disturbance events in palynological sequences #Neotoma dataset id 21665, from publication Behling et al. 1998 behling <- read.csv("dataset21665Behling1998_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(behling) # 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)) behling2<-behling[-c(1:2),-c(3:5)] #set all na values to zero values in the counts behling2[is.na(behling2)]<-0 str(behling2) #to calculate percentages, I need the whole terrestrial pollen dataset behling_terr<-with(behling2,behling2[behling2$group == "TRSH" | behling2$group == "UPHE" | behling2$group == "TRSHA" | behling2$group == "PALM" | behling2$group == "MANG",]) # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) behling_terr_num<-behling_terr[,-c(1:2)] behlingpcts = sapply(names(behling_terr_num), function(x) { (behling_terr[x] / sum(behling_terr_num[x]))*100 }) behlingpcts<-as.data.frame(behlingpcts) behlingpcts[is.na(behlingpcts)]<-0 #I need to reassign sample names, use column (variable-) names from behling2 names(behlingpcts)<-names(behling2[,c(3:ncol(behling2))]) #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 #I need to subsample according to groups again, so I need to reattach group column to dataset behling_pcts_names<-cbind(behling_terr[,1:2],behlingpcts) behling_arb_pcts<-with(behling_pcts_names,behling_pcts_names[behling_pcts_names$group == "TRSH" | behling_pcts_names$group == "PALM",]) behling_arb_pcts_num<-behling_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(behlingpcts) df<-as.data.frame(df) #add species names as columns names #get species names: behlingpcts_ch<-as.character(behling_pcts_names$Behling1998name) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(behlingpcts_ch) # 3. Step: Transform variables 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.56047 -> unimodal 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(behling_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) #no significant zones, but in publication 4 zones -> use 4 zones behling_depth<-behling[1,c(6:ncol(behling))] tbehling_depth<-t(behling_depth) df_depth<-cbind(tbehling_depth,sumperc_ord,df) #windows(width=14, height=10) behling_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(behling_plot,clust,nZone = 2, col="red") behling_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(behling_plot_short,clust,nZone = 4, 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]] #calculate means and sd per zone #Zone 1 rows mean_behling1<-mean(sumperc_ord[1:c(length(zone1))]) sd_behling1<-sd(sumperc_ord[1:c(length(zone1))]) lower_th1<-mean_behling1-sd_behling1 df_lower_th1<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th1,]) df_lower_th1[,1:3] #S199970 20, S199972 30, S199974 40, S199979 65, S199987 105 zone1 #Zone2 rows mean_behling2<-mean(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) sd_behling2<-sd(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) lower_th2<-mean_behling2-sd_behling2 df_lower_th2<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th2,]) df_lower_th2[,1:3] #S199996 150, S199997 155 zone2 #Zone3 rows mean_behling3<-mean(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) sd_behling3<-sd(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) lower_th3<-mean_behling3-sd_behling3 df_lower_th3<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th3,]) df_lower_th3[,1:3] #S200017 255 zone3 #Zone4 rows mean_behling4<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) sd_behling4<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) lower_th4<-mean_behling4-sd_behling4 df_lower_th4<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th4,]) df_lower_th4[,1:3] #0 zone4 #__________________________________________________________________________ #same but not divided by zones: mean_behling<-mean(sumperc_ord) sd_behling<-sd(sumperc_ord) lower_th<-mean_behling-sd_behling df_lower_th<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th,]) df_lower_th [,1:3] #___________________________________________________________________________ # -> local mean - 1 sd per zone ###S199970 20, S199972 30, S199974 40, S199979 65, S199987 105, ##S199996 150, S199997 155, #S200017 255 #without subdivision into zones: S199972 30, S199996 150, S199997 155, S200020 270, S200021 275, S200022 295 behling_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(behling_plot_short,clust,nZone = 4, col="red") addZone(behling_plot_short,upper=27.5, lower=32.5,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, as just one point in disturbance valley addZone(behling_plot_short,upper=147.5, lower=157.5,col=rgb(1,0,0,alpha=0.3)) addZone(behling_plot_short,upper=267.5,lower=295,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, as just one point in disturbance valley 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 # from depths 21 # 1:Fpre: S200001 175 (83.33333%); Fmax: S199994 140 (85.71429%); Fmin: S199997 155 (75.77640%); Trec=Tfmin-Tfmax: 2858-64=247 years 2989/2858/2683-2781/2611/2380 # 2:Fpre: # 3:Fpre: behlingage1<-read.csv("dataset21665Behling1998.csv",sep=",") #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec (((85.71429-75.77640)/(83.33333-75.77640))*100)/247 #RR=24.62367