### ordination analyses for leyden 1998 #in R: "obervations"-> rows, "variables"->columns leyden <- read.csv("dataset2396Leyden1996_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(leyden) # 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)) leyden2<-leyden[-c(1:2),-c(3:5)] #set all na values to zero values in the counts leyden2[is.na(leyden2)]<-0 str(leyden2) #to calculate percentages, I need the whole terrestrial pollen dataset leyden_terr<-with(leyden2,leyden2[leyden2$group == "TRSH" | leyden2$group == "UPHE" | leyden2$group == "TRSHA" | leyden2$group == "PALM" | leyden2$group == "MANG",]) # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) leyden_terr_num<-leyden_terr[,-c(1:2)] leydenpcts = sapply(names(leyden_terr_num), function(x) { (leyden_terr[x] / sum(leyden_terr_num[x]))*100 }) leydenpcts<-as.data.frame(leydenpcts) leydenpcts[is.na(leydenpcts)]<-0 #I need to reassign sample names, use column (variable-) names from leyden2 names(leydenpcts)<-names(leyden2[,c(3:ncol(leyden2))]) #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 leyden_pcts_names<-cbind(leyden_terr[,1:2],leydenpcts) leyden_arb_pcts<-with(leyden_pcts_names,leyden_pcts_names[leyden_pcts_names$group == "TRSH" | leyden_pcts_names$group == "PALM",]) leyden_arb_pcts_num<-leyden_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(leydenpcts) df<-as.data.frame(df) #add species names as columns names #get species names: leydenpcts_ch<-as.character(leyden_pcts_names$Leyden1996name) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(leydenpcts_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: 0.95869 -> 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(leyden_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) #2 zones leyden_depth<-leyden[1,c(6:ncol(leyden))] tleyden_depth<-t(leyden_depth) df_depth<-cbind(tleyden_depth,sumperc_ord,df) #windows(width=14, height=10) leyden_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(leyden_plot,clust,nZone = 2, col="red") leyden_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(leyden_plot_short,clust,nZone = 2, 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_leyden1<-mean(sumperc_ord[1:c(length(zone1))]) sd_leyden1<-sd(sumperc_ord[1:c(length(zone1))]) lower_th1<-mean_leyden1-sd_leyden1 df_lower_th1<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th1,]) df_lower_th1[,1:3] #S48985 45, S48986 49 zone1 #Zone2 rows mean_leyden2<-mean(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) sd_leyden2<-sd(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) lower_th2<-mean_leyden2-sd_leyden2 df_lower_th2<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th2,]) df_lower_th2[,1:3] #S48987 53 , S48993 77 zone2 #Zone3 rows mean_leyden3<-mean(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) sd_leyden3<-sd(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) lower_th3<-mean_leyden3-sd_leyden3 df_lower_th3<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th3,]) df_lower_th3[,1:3] # zone3 #Zone4 rows mean_leyden4<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) sd_leyden4<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) lower_th4<-mean_leyden4-sd_leyden4 df_lower_th4<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th4,]) df_lower_th4[,1:3] # zone4 #Zone5 rows mean_leyden5<-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_leyden5<-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_leyden5-sd_leyden5 df_lower_th5<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th5,]) df_lower_th5[,1:3] # zone5 #Zone6 rows mean_leyden6<-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_leyden6<-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_leyden6-sd_leyden6 df_lower_th6<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th6,]) df_lower_th6[,1:3] # zone6 #__________________________________________________________________________ #same but not divided by zones: #mean_leyden<-mean(sumperc_ord) #sd_leyden<-sd(sumperc_ord) #lower_th<-mean_leyden-sd_leyden #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 ####S48985 45, S48986 49, #S48987 53 , S48993 77 leyden_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(leyden_plot_short,clust,nZone = 2, col="red") addZone(leyden_plot_short,upper=43, lower=55,col=rgb(1,0,0,alpha=0.3)) addZone(leyden_plot_short,upper=75, lower=79,col=rgb(1,0,0,alpha=0.3)) #not as disturbance event, as just one point in disturbance valley, start of disturbance not clear 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 43 # 1:Fpre: S48988 57 (55.49133%); Fmax: S48982 33 (59.19118%); Fmin: S48986 49 (38.07339%); Trec=Tfmin-Tfmax: 445-300=145 years # 3:Fpre: leydenage1<-read.csv("dataset2396Leyden1996.csv",sep=",") #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec #(((85.71429-75.77640)/(83.33333-75.77640))*100)/247 #RR=24.62367