diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index e85854c33a8ef014903c6ca6fdd80fdecadb1ff0..4090a7de39b49347e5acce4f44811b458614b65a 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -163,7 +163,8 @@ species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
traits <- read_delim("_data/Mesobromion/traits.out.10perc.txt", delim="\t")
traits <- traits %>%
- column_to_rownames("species0") # %>%
+ column_to_rownames("species0") %>%
+ dplyr::select(LeafArea:Disp.unit.leng)
#dplyr::select(PlantHeight, LeafC.perdrymass, LeafN, StemDens, Stem.cond.dens, Seed.num.rep.unit, SLA)
@@ -172,22 +173,37 @@ traits <- traits %>%
myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_[0-9]+.RData", full.names = T)
corXY <- NULL
corXY.perm <- NULL
+corXY.bootstr <- NULL
for(ff in myfilelist){
index <- as.numeric(regmatches(ff, gregexpr("[[:digit:]]+", ff))[[1]])
load(ff)
corXY <- rbind(corXY, cor.obs)
corXY.perm <- rbind(corXY.perm, cor.perm)
+ corXY.bootstr <- rbind(corXY.bootstr, cor.bootstr)
print(index)
}
+
+diffs <- corXY.bootstr - corXY.perm
+sign <- apply(diffs, MARGIN=1, function(x){sum(x>0)>190})
+sign <- data.frame(Trait.comb=names(sign), sign=sign)
+
corXY %>%
arrange(Test, desc(Coef))
-aa <- data.frame(Trait.comb=paste0("t", 1:80),
+aa <- data.frame(Trait.comb=paste0("t", 1:15),
trait.name=colnames(traits))#[-which(colnames(traits) %in% c("species", "species0"))])
-bb <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>%
- separate(Trait.comb, c("trait1", "trait2")) %>%
+
+corXY.ci <- data.frame(Trait.comb=rownames(corXY.bootstr),
+ obs=corXY.bootstr[,1],
+ q025=apply(corXY.bootstr[,-1], MARGIN=1, quantile, probs=0.025),
+ q975=apply(corXY.bootstr[,-1], MARGIN=1, quantile, probs=0.975)) %>%
+ left_join(sign, by="Trait.comb") %>%
+ as.tbl() %>%
+ mutate(Trait.comb2=Trait.comb) %>%
+ separate(Trait.comb2, c("trait1", "trait2", "trait3")) %>%
mutate(trait2=paste0("t", trait2)) %>%
+ mutate(trait3=paste0("t", trait3)) %>%
left_join(aa %>%
rename(trait1=Trait.comb), by="trait1") %>%
dplyr::select(-trait1) %>%
@@ -195,12 +211,63 @@ bb <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>%
left_join(aa %>%
rename(trait2=Trait.comb), by="trait2") %>%
dplyr::select(-trait2) %>%
- rename(trait2=trait.name) %>%
- dplyr::select(trait1, trait2, q025:conf.p) %>%
+ rename(trait2=trait.name) %>%
+ left_join(aa %>%
+ rename(trait3=Trait.comb), by="trait3") %>%
+ dplyr::select(-trait3) %>%
+ rename(trait3=trait.name) %>%
+ mutate(ntraits= (!is.na(trait1)) + (!is.na(trait2)) + (!is.na(trait3))) %>%
+ arrange(ntraits, obs)
+
+### spider plot
+# faceted spider plot
+mydata <- corXY.ci %>%
+ filter(grepl(pattern = "t2", Trait.comb)) %>%
+ #filter(grepl(pattern = "_6", Trait.comb)) %>%
+ filter(sign==T) %>%
+ filter(ntraits<=2) %>%
+ mutate(Trait.comb=factor(Trait.comb)) %>%
+ mutate(seq=1:nrow(.))
+
+ggplot(data=mydata) +
+ geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq1, col=sign)) +
+ geom_point(aes(x=obs, y=seq), pch=15) +
+ scale_y_continuous(breaks=mydata$seq,
+ labels=mydata$Trait.comb)
+
+
+
+
+
+
+corXY.SES <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>%
+ separate(Trait.comb, c("trait1", "trait2", "trait3")) %>%
+ mutate(trait2=paste0("t", trait2)) %>%
+ mutate(trait3=paste0("t", trait3)) %>%
+ left_join(aa %>%
+ rename(trait1=Trait.comb), by="trait1") %>%
+ dplyr::select(-trait1) %>%
+ rename(trait1=trait.name) %>%
+ left_join(aa %>%
+ rename(trait2=Trait.comb), by="trait2") %>%
+ dplyr::select(-trait2) %>%
+ rename(trait2=trait.name) %>%
+ left_join(aa %>%
+ rename(trait3=Trait.comb), by="trait3") %>%
+ dplyr::select(-trait3) %>%
+ rename(trait3=trait.name) %>%
+ dplyr::select(trait1, trait2, trait3, q025:conf.p) %>%
arrange(desc(SES.np))
print(bb, n=20)
+
+
+break()
+
+
+
+
#### Calculate and plot PCA of CWMs
CWM.wide <- species %>%
rownames_to_column("RELEVE_NR") %>%
@@ -239,75 +306,32 @@ ggplot() +
ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
theme_bw() +
theme(panel.grid = element_blank()) +
- ggtitle("PCA of CWMs (Axis 1 and 2)")
-
-
-ggplot() +
- geom_point(data=as.data.frame(CWM.pca$CA$u[,c(1,3)]),
- aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(CWM.pca$CA$v*2) %>%
- rownames_to_column("Trait") %>%
- mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
- aes(x=0, xend=PC1, y=0, yend=PC3, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_text(data=as.data.frame(CWM.pca$CA$v*2.1) %>%
- rownames_to_column("Trait") %>%
- mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
- aes(x=PC1, y=PC3, label=Trait, col=in.try), size=3 ) +
- scale_color_identity() +
- xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PC3 (", varexpl[3], "%)", sep="")) +
- theme_bw() +
- theme(panel.grid = element_blank()) +
- ggtitle("PCA of CWMs (Axis 1 and 3)")
-
-
-
-
+ title("PCA of CWMs")
#PCA of individual traits
trait.pca <- vegan::rda(as.matrix(traits %>%
- filter(!is.na(rowSums(.)))), na.action=na.omit,
- scale=T)
+ filter(!is.na(rowSums(.)))), na.action=na.omit)
varexpl.t <- round(trait.pca$CA$eig/sum(trait.pca$CA$eig)*100,1)
ggplot() +
geom_point(data=as.data.frame(trait.pca$CA$u[,1:2]),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(trait.pca$CA$v*1.3) %>%
+ geom_segment(data=as.data.frame(trait.pca$CA$v*2) %>%
rownames_to_column("Trait") %>%
mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_text(data=as.data.frame(trait.pca$CA$v*1.4) %>%
+ geom_text(data=as.data.frame(trait.pca$CA$v*2.1) %>%
rownames_to_column("Trait") %>%
mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
aes(x=PC1, y=PC2, label=Trait, col=in.try), size=3 ) +
scale_color_identity() +
- xlab(paste("PC1 (", varexpl.t[1], "%)", sep="")) +
- ylab(paste("PC2 (", varexpl.t[2], "%)", sep="")) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
theme_bw() +
theme(panel.grid = element_blank()) +
- ggtitle("PCA of species-level traits") +
- coord_equal()
-
-trait.envf <- envfit(trait.pca, traits %>%
- filter(!is.na(rowSums(.))), choices = c(1,2,3), na.rm = T)
+ title("PCA of species-level traits")
-as.data.frame(trait.envf$vectors$arrows) %>%
- rownames_to_column("Trait") %>%
- arrange(desc(PC1))
-
-
-out <- NULL
-for(tt1 in 2:ncol(traits)){
- for(tt2 in 1:ncol(traits)) {
- #for(tt2 in 1:(tt1-1)) {
- tmp <- cor.test(traits[,tt1], traits[,tt2], use="pairwise.complete.obs")
- out <- rbind(out,
- data.frame(trait1=colnames(traits)[tt1] , trait2=colnames(traits)[tt2], r=tmp$estimate, pvalue.r=tmp$p.value))
- }
-}
-bb %>% left_join(out, by=c("trait1", "trait2"))
#### Map of plots
diff --git a/01b_MesobromionCluster.R b/01b_MesobromionCluster.R
index 7fef7c6d37173629f59b39eeb6520c540fee4bd1..c7a92aae76755b11c59d9b3bfd4f4b6c33b09f9c 100644
--- a/01b_MesobromionCluster.R
+++ b/01b_MesobromionCluster.R
@@ -14,148 +14,8 @@ require(ade4)
require(parallel)
require(doParallel)
-#### Function 1 - CorXY ####
-get.corXY <- function(comm, traits, trait.sel="all", stat=c("mantel", "RV", "procrustes")){
- if(identical(trait.sel, "all")) {trait.sel <- 1:ncol(traits)}
- ii <- trait.sel
- lab.tmp <- paste(ii, collapse="_")
- ### delete potential missing species
- if(any(colSums(comm)==0)){
- empty <- which(colSums(comm)==0)
- traits <- traits[-empty,]
- comm <- comm[,-empty]
- }
- syn.out.tmp <- matrix.x(comm=comm, traits=traits[,ii,drop=F], scale=T)$matrix.X
- W.beals <- as.data.frame(beals(comm, include=T, type=2))
- corXY <- NULL
- if("mantel" %in% stat){
- W.beals.d <- dist(W.beals)
- mantel.tmp <- mantel(W.beals.d, dist(syn.out.tmp[]))
- corXY <- rbind(corXY,
- data.frame(Trait.comb=lab.tmp, Test="Mantel", Coef=mantel.tmp$statistic, pvalue=mantel.tmp$signif))
- }
- if("RV" %in% stat){
- RV.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.tmp))
- corXY <- rbind(corXY,
- data.frame(Trait.comb=lab.tmp, Test="RV", Coef=RV.tmp$obs, pvalue=RV.tmp$pvalue))
- }
- if("procrustes" %in% stat){
- prot.tmp <- protest(W.beals, syn.out.tmp)
- corXY <- rbind(corXY,
- data.frame(Trait.comb=lab.tmp, Test="Procrustes", Coef=prot.tmp$t0, pvalue=prot.tmp$signif))
- }
- return(corXY)
-}
-
-
-
-#### Function 2 - CorTE ####
-get.corTE <- function(comm, traits, envir, trait.sel="all", env.sel="all", stat=c("mantel", "RV")){
- if(identical(trait.sel, "all")) {trait.sel <- 1:ncol(traits)}
- if(identical(env.sel, "all" )) {env.sel <- 1:ncol(envir)}
- ii <- trait.sel
- lab.tmp <- paste(ii, collapse="_")
- ## delete potential missing species
- if(any(colSums(comm)==0)){
- empty <- which(colSums(comm)==0)
- traits <- traits[-empty,]
- comm <- comm[,-empty]
- }
- syn.out.tmp <- matrix.t(comm=comm, traits=traits, scale=T)$matrix.T
- ee <- env.sel
- lab.env <- paste(ee, collapse="_")
- corTE <- NULL
- if("mantel" %in% stat){
- mantel.tmp <- mantel(dist(envir[,ee, drop=F]), dist(syn.out.tmp[,ii,drop=F]))
- corTE <- rbind(corTE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="Mantel", Coef=mantel.tmp$statistic, pvalue=mantel.tmp$signif))
- }
- if("RV" %in% stat){
- RV.tmp <- RV.rtest(as.data.frame(envir[,ee, drop=F]), as.data.frame(syn.out.tmp[,ii,drop=F]))
- corTE <- rbind(corTE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="RV", Coef=RV.tmp$obs, pvalue=RV.tmp$pvalue))
- }
- if("procrustes" %in% stat){
- prot.tmp <- protest(envir[,ee, drop=F], syn.out.tmp[,ii,drop=F])
- corTE <- rbind(corTE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="Procrustes", Coef=prot.tmp$t0, pvalue=prot.tmp$signif))
- }
- return(corTE)
-}
-
-
-#### Function 3 - CorXE ####
-get.corXE <- function(comm, traits, envir, trait.sel="all", env.sel="all", stat=c("mantel", "RV", "procrustes")){
- if(identical(trait.sel, "all")) {trait.sel <- 1:ncol(traits)}
- if(identical(env.sel, "all" )) {env.sel <- 1:ncol(envir)}
- ii <- trait.sel
- lab.tmp <- paste(ii, collapse="_")
- ### delete potential missing species
- if(any(colSums(comm)==0)){
- empty <- which(colSums(comm)==0)
- traits <- traits[-empty,]
- comm <- comm[,-empty]
- }
- syn.out.tmp <- matrix.x(comm=comm, traits=traits[,ii,drop=F], scale=T)$matrix.X
- ee <- env.sel
- lab.env <- paste(ee, collapse="_")
- corXE <- NULL
- if("mantel" %in% stat){
- mantel.tmp <- mantel(dist(envir[,ee, drop=F]), dist(syn.out.tmp[]))
- corXE <- rbind(corXE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="Mantel", Coef=mantel.tmp$statistic, pvalue=mantel.tmp$signif))
- }
- if("RV" %in% stat){
- RV.tmp <- RV.rtest(as.data.frame(envir[,ee, drop=F]), as.data.frame(syn.out.tmp))
- corXE <- rbind(corXE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="RV", Coef=RV.tmp$obs, pvalue=RV.tmp$pvalue))
- }
- if("procrustes" %in% stat){
- prot.tmp <- protest(envir[,ee, drop=F], syn.out.tmp)
- corXE <- rbind(corXE,
- data.frame(Trait.comb=lab.tmp, Env.comb=lab.env, Test="Procrustes", Coef=prot.tmp$t0, pvalue=prot.tmp$signif))
- }
- return(corXE)
-}
-
-
-
-### Get SES (both parametric and non parametric)
-#obs.df = output from one of the get functions above
-#perm.df = df traitCombination X permutations with the chosen statistic of correlation on permuteda data
-
-get.SES <- function(obs.df, perm.df, stat="RV") {
- nperm <- dim(perm.df)[2]
- SES.out <- as.data.frame.table(perm.df) %>%
- rename(Trait.comb=Var1, perm=Var2, coef=Freq ) %>%
- group_by(Trait.comb) %>%
- left_join(obs.df %>%
- filter(Test==stat) %>%
- dplyr::select(-pvalue, -Test) %>%
- mutate(Trait.comb=paste0("t", Trait.comb)) %>%
- rename(obs=Coef),
- by=c("Trait.comb")) %>%
- summarize(q15 = quantile(coef, probs = 0.15865),
- q50 = quantile(coef, .5),
- q84 = quantile(coef, 0.84135),
- mean.perm=mean(coef),
- sd.perm=sd(coef),
- greater.than.obs=sum( (abs(coef)>=abs(obs))/n())) %>%
- left_join(obs.df %>%
- filter(Test==stat) %>%
- dplyr::select(-pvalue, -Test) %>%
- mutate(Trait.comb=paste0("t", Trait.comb)) %>%
- rename(obs=Coef),
- by=c("Trait.comb")) %>%
- mutate(SES.np= ifelse(obs>=q50, (obs-q50)/(q84-q50), (obs-q50)/(q50-q15))) %>%
- mutate(SES=(obs-mean.perm)/sd.perm) %>%
- arrange(desc(SES.np)) %>%
- mutate(conf.m=obs-1.65*(sd.perm/sqrt(nperm))) %>%
- mutate(conf.p=obs+1.65*(sd.perm/sqrt(nperm)))
- return(SES.out)
-}
-
+source("99_HIDDEN_functions.R")
Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY", max.inter.t, chunkn, chunk.i, nperm=99){