### ordination analyses for larmon dataset larmon <- read.csv("dataset15099Larmon2013_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(larmon) # 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)) larmon2<-larmon[-c(1:2),-c(3:5)] #set all na values to zero values in the counts larmon2[is.na(larmon2)]<-0 str(larmon2) #to calculate percentages, I need the whole terrestrial pollen dataset larmon_terr<-with(larmon2,larmon2[larmon2$group == "TRSH" | larmon2$group == "UPHE" | larmon2$group == "PALM" | larmon2$group == "MANG" ,]) # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) larmon_terr_num<-larmon_terr[,-c(1:2)] larmonpcts = sapply(names(larmon_terr_num), function(x) { (larmon_terr[x] / sum(larmon_terr_num[x]))*100 }) larmonpcts<-as.data.frame(larmonpcts) larmonpcts[is.na(larmonpcts)]<-0 #I need to reassign sample names, use column (variable-) names from larmon2 names(larmonpcts)<-names(larmon2[,c(3:ncol(larmon2))]) #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 larmon_pcts_names<-cbind(larmon_terr[,1:2],larmonpcts) larmon_arb_pcts<-with(larmon_pcts_names,larmon_pcts_names[larmon_pcts_names$group == "TRSH" | larmon_pcts_names$group == "PALM",]) larmon_arb_pcts_num<-larmon_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(larmonpcts) df<-as.data.frame(df) #add species names as columns names #get species names: larmonpcts_ch<-as.character(larmon_pcts_names$Larmon2013name) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(larmonpcts_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.8277 -> 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(larmon_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) #4 zones larmon_depth<-larmon[1,c(6:ncol(larmon))] tlarmon_depth<-t(larmon_depth) df_depth<-cbind(tlarmon_depth,sumperc_ord,df) larmon_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(larmon_plot,clust,nZone = 4, col="red") larmon_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(larmon_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 mean_larmon1<-mean(sumperc_ord[1:c(length(zone1))]) sd_larmon1<-sd(sumperc_ord[1:c(length(zone1))]) lower_th1<-mean_larmon1-sd_larmon1 df_lower_th1<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th1,]) df_lower_th1 #S136550 190.5 , S136551 200.5 #Zone2 mean_larmon2<-mean(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) sd_larmon2<-sd(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) lower_th2<-mean_larmon2-sd_larmon2 df_lower_th2<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th2,]) df_lower_th2 #no #Zone3 mean_larmon3<-mean(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) sd_larmon3<-sd(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) lower_th3<-mean_larmon3-sd_larmon3 df_lower_th3<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th3,]) df_lower_th3 #S136573 470.5, S136574 480.5 #Zone4 mean_larmon4<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) sd_larmon4<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) lower_th4<-mean_larmon4-sd_larmon4 df_lower_th4<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th4,]) df_lower_th4 #S136579 530.5 #__________________________________________________________________________ #same but not divided by zones: #mean_larmon<-mean(sumperc_ord) #sd_larmon<-sd(sumperc_ord) #lower_th<-mean_larmon-sd_larmon #df_lower_th<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th,]) #df_lower_th #___________________________________________________________________________ df_depth[,1:3] larmon_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(larmon_plot_short,clust,nZone = 4, col="red") addZone(larmon_plot,upper=185.5, lower=205.5,col=rgb(1,0,0,alpha=0.3)) addZone(larmon_plot,upper=460.5, lower=485.5,col=rgb(1,0,0,alpha=0.3)) addZone(larmon_plot,upper=525.5, lower=535.5,col=rgb(1,0,0,alpha=0.3)) #not as disturbance event, as just one point in disturbance valley #Choose disturbance events where more than one sample is below the local threshold within the same "valley of percentages" (is this justifiable?) #Identified one zone of disturbance for sample-depths at ca. 160-260, as it is in the middle of the sequence and a recovery can be observed. # Recovery rate calculation according to Cole et al. 2014: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec # Fmax: maximum forest percentage after recovery # Fmin: lowest percentage of forest # Fpre: Forest percentage pre-disturbance # Trec: Time from Fmin to Fmax # From Cole et al. 2014: #Definition recovery rate: "Rate of increase in forest abundance relative to degree of #disturbance-induced change, that is, % increase in forest pollen #abundance per year in relation to pre-disturbance level" #"RR is measured as an average across the entire period of forest pollen increase (representing forest #re-establishment), until the maximum increase in percentage of forest pollen #post-decline is achieved, before a stabilizing point or further decline. It therefore #incorporates any non-linear components of the trajectory, which may result, for #example, from low-level perturbations during the period of recovery; it is, however, #the overall rate at which the forest is responding to each disturbance event that is #the variable of interest here." # 1:Fpre: S136555 (36.073059%); Fmax: S136549 (22.222222%); Fmin: S136551 (1.158301%); Trec: 952-723=229 years # 2:Fpre: S136575 490.5 (7.773852%); Fmax: S136572 450.5 (16.666667%); Fmin: S136574 480.5 (3.030303%); Trec: 2678-2585=93 years (Tfmin-Tfmax) 2886/2678/2521-2772/2585/2420 #not as disturbance event, as just one point in disturbance valley# 3:Fpre: S136580 540.5 (6.185567%); Fmax: S136578 520.5 (9.395973%); Fmin: S136579 530.5 (1.075269%); Trec: 2836-2804=32 years, 3063/2836/2715-3024/2804/2667 larmonage1<-read.csv("dataset15099Larmon2013.csv",sep=",") #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec (((9.395973-1.075269)/(6.185567-1.075269))*100)/32 #RR=0.2634479, 3.091094, 5.088196