All functions moved to 99

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

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){