### ordination analyses for meliefcleef 2008 Primavera I #in R: "obervations"-> rows, "variables"->columns meliefcleef <- read.csv("dataset21917Melief_Cleef2008Primavera_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(meliefcleef) # 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)) meliefcleef2<-meliefcleef[-c(1:2),-c(3:5)] #set all na values to zero values in the counts meliefcleef2[is.na(meliefcleef2)]<-0 str(meliefcleef2) #to calculate percentages, I need the whole terrestrial pollen dataset meliefcleef_terr<-with(meliefcleef2,meliefcleef2[meliefcleef2$group == "TRSH" | meliefcleef2$group == "UPHE" | meliefcleef2$group == "TRSHA" | meliefcleef2$group == "TRSHSP" | meliefcleef2$group == "TRSHP"| meliefcleef2$group == "PALM" | meliefcleef2$group == "MANG" | meliefcleef2$group == "SUCC",]) # 2. step: I need to calculate percentages, therefore I need to remove the non-numeric columns (taxa name and group) meliefcleef_terr_num<-meliefcleef_terr[,-c(1:2)] meliefcleefpcts = sapply(names(meliefcleef_terr_num), function(x) { (meliefcleef_terr[x] / sum(meliefcleef_terr_num[x]))*100 }) meliefcleefpcts<-as.data.frame(meliefcleefpcts) meliefcleefpcts[is.na(meliefcleefpcts)]<-0 #I need to reassign sample names, use column (variable-) names from meliefcleef2 names(meliefcleefpcts)<-names(meliefcleef2[,c(3:ncol(meliefcleef2))]) #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 meliefcleef_pcts_names<-cbind(meliefcleef_terr[,1:2],meliefcleefpcts) meliefcleef_arb_pcts<-with(meliefcleef_pcts_names,meliefcleef_pcts_names[meliefcleef_pcts_names$group == "TRSH" | meliefcleef_pcts_names$group == "PALM",]) meliefcleef_arb_pcts_num<-meliefcleef_arb_pcts[,-c(1:2)] # 2.2: I needed to transpose my dataframe df <- t(meliefcleefpcts) df<-as.data.frame(df) #add species names as columns names #get species names: meliefcleefpcts_ch<-as.character(meliefcleef_pcts_names$Melief_Cleef2008Primaveraname) ###!!! Manually edit this column name every time you change a site!!! colnames(df)<-c(meliefcleefpcts_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.28925 -> 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(meliefcleef_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) #5 zones meliefcleef_depth<-meliefcleef[1,c(6:ncol(meliefcleef))] tmeliefcleef_depth<-t(meliefcleef_depth) df_depth<-cbind(tmeliefcleef_depth,sumperc_ord,df) #windows(width=14, height=10) meliefcleef_plot<-strat.plot(df_depth[,c(2:ncol(df_depth))], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(meliefcleef_plot,clust,nZone = 5, col="red") meliefcleef_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(meliefcleef_plot_short,clust,nZone = 5, 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_meliefcleef1<-mean(sumperc_ord[1:c(length(zone1))]) sd_meliefcleef1<-sd(sumperc_ord[1:c(length(zone1))]) lower_th1<-mean_meliefcleef1-sd_meliefcleef1 df_lower_th1<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th1,]) df_lower_th1[,1:3] #S211495 55,S211496 64,S211510 150,S211511 156,S211523 223,S211530 261,S211531 274 zone1 #Zone2 rows mean_meliefcleef2<-mean(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) sd_meliefcleef2<-sd(sumperc_ord[c(length(zone1)+1):c(length(zone1)+length(zone2))]) lower_th2<-mean_meliefcleef2-sd_meliefcleef2 df_lower_th2<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th2,]) df_lower_th2[,1:3] #S211538 323,S211547 378,S211549 387 zone2 #Zone3 rows mean_meliefcleef3<-mean(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) sd_meliefcleef3<-sd(sumperc_ord[c(length(zone1)+length(zone2)+1):c(length(zone1)+length(zone2)+length(zone3))]) lower_th3<-mean_meliefcleef3-sd_meliefcleef3 df_lower_th3<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th3,]) df_lower_th3[,1:3] #S211576 553 zone3 #Zone4 rows mean_meliefcleef4<-mean(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) sd_meliefcleef4<-sd(sumperc_ord[c(length(zone1)+length(zone2)+length(zone3)+1):c(length(zone1)+length(zone2)+length(zone3)+length(zone4))]) lower_th4<-mean_meliefcleef4-sd_meliefcleef4 df_lower_th4<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th4,]) df_lower_th4[,1:3] #S211598 703 zone4 #Zone5 rows mean_meliefcleef5<-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_meliefcleef5<-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_meliefcleef5-sd_meliefcleef5 df_lower_th5<-with(df_depth,df_depth[df_depth$sumperc_ord < lower_th5,]) df_lower_th5[,1:3] #S211602 731,S211606 752,S211611 774,S211614 791 zone5 #Zone6 rows mean_meliefcleef6<-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_meliefcleef6<-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_meliefcleef6-sd_meliefcleef6 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_meliefcleef<-mean(sumperc_ord) #sd_meliefcleef<-sd(sumperc_ord) #lower_th<-mean_meliefcleef-sd_meliefcleef #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 #S211495 55,S211496 64,S211510 150,S211511 156,S211523 223,S211530 261,S211531 274 #S211538 323,S211547 378,S211549 387 #S211576 553 #S211598 703 #S211602 731,S211606 752,S211611 774,S211614 791 meliefcleef_plot_short<-strat.plot(df_depth[,2:10], yvar = df_depth[,1],y.rev = TRUE, scale.percent = TRUE) addClustZone(meliefcleef_plot_short,clust,nZone = 5, col="red") addZone(meliefcleef_plot_short,upper=54, lower=56,col=rgb(1,0,0,alpha=0.3)) addZone(meliefcleef_plot_short,upper=63, lower=65,col=rgb(1,0,0,alpha=0.3))#dist 1 addZone(meliefcleef_plot_short,upper=149, lower=157,col=rgb(1,0,0,alpha=0.3))#dist 2 addZone(meliefcleef_plot_short,upper=222, lower=224,col=rgb(1,0,0,alpha=0.3))#dist 3 addZone(meliefcleef_plot_short,upper=260, lower=262,col=rgb(1,0,0,alpha=0.3))#dist 4 addZone(meliefcleef_plot_short,upper=273, lower=275,col=rgb(1,0,0,alpha=0.3)) addZone(meliefcleef_plot_short,upper=322, lower=324,col=rgb(1,0,0,alpha=0.3))#dist 5 addZone(meliefcleef_plot_short,upper=377, lower=379,col=rgb(1,0,0,alpha=0.3))#dist 6 addZone(meliefcleef_plot_short,upper=386, lower=388,col=rgb(1,0,0,alpha=0.3)) addZone(meliefcleef_plot_short,upper=552, lower=554,col=rgb(1,0,0,alpha=0.3))#dist 7 addZone(meliefcleef_plot_short,upper=702, lower=704,col=rgb(1,0,0,alpha=0.3))#dist 8 addZone(meliefcleef_plot_short,upper=730, lower=732,col=rgb(1,0,0,alpha=0.3)) addZone(meliefcleef_plot_short,upper=751, lower=753,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, as just one point in disturbance valley addZone(meliefcleef_plot_short,upper=773, lower=775,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, as just one point in disturbance valley addZone(meliefcleef_plot_short,upper=790, lower=792,col=rgb(1,0,0,alpha=0.3))#not as disturbance event, 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 meliefcleefage1<-read.csv("dataset21917Melief_Cleef2008Primavera.csv",sep=",") #RR: RR=(((Fmax–Fmin)/(Fpre–Fmin))*100)/Trec #see excel fiel RecoveryRates_BiolLetters for calculations