diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index 22d7759830ee4e069be814e9c61d5ae537cb9247..79ac82d5f1e298d09abcaa8b63d9e7d29692839a 100644
--- a/01_Mesobromion.R
+++ b/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