Improved data and trait pre-selection

01_Mesobromion.R

@@ -4,9 +4,11 @@ library(foreign)
 source("99_HIDDEN_functions.R")
 
 
-#### 1. traits data 
+#### 1. traits data  #### 
 
 traits0 <- read_delim("_data/Mesobromion/traits3.txt", delim =";", col_names = T, locale = locale(encoding = 'latin1')) %>% 
+  dplyr::select(- c("PR_STAT_Indigen", "PR_STAT_Neophyt", "PR_STAT_Archaeophyt", "URBAN_urbanophob", "URBAN_maessig_urbanophob", 
+                    "URBAN_urbanoneutral", "URBAN_maessig_urbanophil", "URBAN_urbanophil", "WUH_von", "WUH_bis", "ARL_c_I_von", "ARL_c_I_bis")) %>% 
   mutate(species0=species) %>% 
   rowwise() %>% 
   # quick and dirty clean up names
@@ -14,18 +16,22 @@ traits0 <- read_delim("_data/Mesobromion/traits3.txt", delim =";", col_names = T
   mutate(species=gsub(pattern=" agg | x | spec$| agg$| s | Sec |  ", replacement=" ", x=species)) %>% 
   mutate(species=gsub(pattern=" $", replacement = "", x = species)) %>% 
   mutate(species=ifelse(is.na(word(species, 1, 2)), species, word(species, 1, 2)))
-str(traits0) #907 obs. of  85 variables:
+dim(traits0) #907 obs. of  75 variables:
+
+
 
 #keep only traits with >=88 completeness
 traits0 <- traits0 %>% 
-  dplyr::select_if(~mean(!is.na(.)) >= 0.88) #907 obs 79 variables
+  dplyr::select_if(~mean(!is.na(.)) >= 0.88) # 907 x 67
 
 
 ## Import traits from TRY and match to species
 load("/data/sPlot2.0/TRY.all.mean.sd.3.by.genus.species.tree.Rdata")
 alltry <- TRY.all.mean.sd.3.by.genus.species.tree %>% 
   dplyr::select(!ends_with(".sd")) %>% 
-  dplyr::select(StandSpeciesName, LeafArea.mean:Wood.vessel.length.mean)
+  dplyr::select(StandSpeciesName, LeafArea.mean:Wood.vessel.length.mean) %>% 
+  dplyr::select(-Wood.vessel.length.mean, -StemDens.mean, -Stem.cond.dens.mean) %>% 
+  rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.))
 
 traits <- traits0 %>% 
   ungroup() %>% 
@@ -33,18 +39,18 @@ traits <- traits0 %>%
   left_join(alltry %>% 
               rename(species=StandSpeciesName), 
             by="species") %>% 
-  filter(!is.na(LeafArea.mean))
-dim(traits) #[1] 805  97
+  filter(!is.na(LeafArea))
+dim(traits) #[1] 805  2
 
 
 
-##### 2. Header Data
+##### 2. Header Data #### 
 env0 <- read_delim("_data/Mesobromion/GVRD_MES2_site.csv", delim = ",")
 str(env0) #6868 obs. of  6 variables:
 
 set.seed(1984)
 header   <- "/data/sPlot/users/Francesco/Project_11/Germany/_data/tvhabita.dbf"
-env.all  <- env0 %>% 
+env  <- env0 %>% 
   left_join(foreign::read.dbf(header) %>% 
               as.data.frame() %>% 
               dplyr::select(RELEVE_NR, LAT, LON),
@@ -52,13 +58,9 @@ env.all  <- env0 %>%
   filter(!is.na(LAT)) %>%
   filter(!(LAT==0 | LON==0)) 
 
-env <- env.all
-#env <- env.all %>% 
-#  sample_frac(0.1) #Take 10% subset
 
-  
 
-### 2. Import species data
+### 3. Import species data #### 
 # columns in species correspond to those in env
 # there is no PlotObservationID (yet)
 species0 <- read.table("_data/Mesobromion/GVRD_Mes2_veg1.csv", sep = ",", header=T)
@@ -71,10 +73,12 @@ species <- env %>%
   left_join(species0 %>%
               mutate(RELEVE_NR=env0$RELEVE_NR), 
             by="RELEVE_NR") %>% 
-  column_to_rownames("RELEVE_NR")# %>% 
+  column_to_rownames("RELEVE_NR") %>% 
+  ## delete species not appearing in any plot
+  dplyr::select(colnames(.)[which(colSums(.)!=0)])
   #dplyr::select(traits$species0)
 
-dim(species) # [1] 5810 907 
+dim(species) # [1] 5810 881
 
 releve08trait <- species %>% 
   rownames_to_column("RELEVE_NR") %>% 
@@ -86,17 +90,17 @@ releve08trait <- species %>%
   ## attach traits 
   left_join(traits %>% 
               dplyr::select(-species), by="species0") %>% 
-  mutate_at(.vars = vars(PR_STAT_Indigen:Wood.vessel.length.mean), 
+  mutate_at(.vars = vars(LEB_F_Makrophanerophyt:Disp.unit.leng), 
             .funs = list(~if_else(is.na(.),0,1) * pres)) %>%
   group_by(RELEVE_NR) %>%
-  summarize_at(.vars= vars(PR_STAT_Indigen:Wood.vessel.length.mean),
+  summarize_at(.vars= vars(LEB_F_Makrophanerophyt:Disp.unit.leng),
                .funs = list(~mean(.))) %>%
   dplyr::select(RELEVE_NR, order(colnames(.))) %>%
   reshape2::melt(id.vars="RELEVE_NR", value.name="trait.coverage") %>% 
   group_by(RELEVE_NR) %>% 
-  summarize(ntraits08=sum(trait.coverage>=0.8)) %>% 
-  #select only those traits where at least 75 traits have a coverage of >0.8 of the species
-  #filter(ntraits08==95) %>% dim()
+  summarize(ntraits08=mean(trait.coverage>=0.8)) %>% 
+  #select only those releves where we have a coverage of >0.8 for all traits
+  filter(ntraits08==1) %>% 
   pull(RELEVE_NR)
 
 set.seed(1984)
@@ -105,9 +109,8 @@ species <- species %>%
   rownames_to_column("RELEVE_NR") %>% 
   filter(RELEVE_NR %in% releve08trait.samp) %>% 
   column_to_rownames("RELEVE_NR") %>% 
-  dplyr::select(traits %>% pull(species0)) %>% 
-  #select only species occurring it at least one plot
-  dplyr::select(names(which(colSums(.)!=0)) )
+  as.tbl() %>% 
+  dplyr::select(one_of(traits$species0))
 
 
 env <- env %>% 
@@ -118,9 +121,9 @@ traits <- traits %>%
   dplyr::select(species0, everything()) %>% 
   filter(species0 %in% colnames(species))
   
-dim(species) #581 508
-dim(traits) #508  96
-dim(env) #581   8
+dim(species) #558 783
+dim(traits) #783 81
+dim(env) #558 8
 
 
 ## Select only plots where >90% of species have trait info [TRY]
@@ -149,54 +152,42 @@ write_delim(species, path="_data/Mesobromion/species.out.10perc.txt", delim="\t"
 write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
 
 
+#### CORRELATION BETWEEN FUZZY WEIGHTED AND BEALS MATRICES  #### #### #### 
+#### WAS RUN IN THE CLUSTER WITH THE SCRIPT 01b_MesobromionCluster.R
 
 
-
-#############################
+#### ## Reimport data #### 
+############################
 ## calculate corr between species composition matrix and traits
 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") %>% 
-  rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.))# %>% 
+  column_to_rownames("species0") # %>% 
   #dplyr::select(PlantHeight, LeafC.perdrymass, LeafN, StemDens, Stem.cond.dens, Seed.num.rep.unit, SLA)
 
-## create list of indices for each combination of traits up to a max number of interactions
-n.traits <- ncol(traits)
-max.inter.t <- 1
-allcomb.t <- NULL
-for(n.inter in 1:max.inter.t){
-  allcomb.t <- c(allcomb.t, combn(1:n.traits, n.inter, simplify=F))
-}
-
-
-##Run on observed data
 
-    #Define get.corXY first!!
-require(parallel)
-require(doParallel)
-cl <- makeForkCluster(8, outfile="")
-registerDoParallel(cl)
+#### ## Import output #### 
 
+myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_[0-9]+.RData", full.names = T)
 corXY <- NULL
-corXY <- foreach(i=1:length(allcomb.t), .packages=c('SYNCSA', "vegan", "ade4", "tidyverse"), .combine=rbind) %do% {
-#for(i in 1:length(allcomb.t)){
-  tt <- unlist(allcomb.t[i])
-  corXY <- get.corXY(comm=species, trait=traits, trait.sel=tt, stat="RV")
-  print(i)
-  save(corXY, file=paste0("_data/Mesobromion/corXY/corXY_",i, ".RData"))
-  return(corXY)
+corXY.perm <- 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)
+  print(index)
 }
-stopCluster(cl)
-
-
-corXY %>% arrange(Test, desc(Coef))
-save(corXY, file="_data/Mesobromion/corXY/corXY.RData")
-
 
-aa <- matrix.x(comm=species, trait=traits, trait.sel=tt, stat="RV")
+corXY %>% 
+  arrange(Test, desc(Coef))
 
+aa <- data.frame(Trait.comb=paste0("t", 1:95), trait.name=colnames(traits)[-which(colnames(traits) %in% c("species", "species0"))])
+bb <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>% 
+  left_join(aa, by="Trait.comb") %>% 
+  arrange(desc(SES.np))
+print(bb, n=20)
 
 
 #### Map of plots