diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index fe812adb2e89563f20cc5194612e11e1db0c6e0e..54ed6dc289abd0ad239fa31e708226269af355a0 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -11,7 +11,8 @@ source("99_HIDDEN_functions.R")
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")) %>%
+ "URBAN_urbanoneutral", "URBAN_maessig_urbanophil", "URBAN_urbanophil", "WUH_von", "WUH_bis", "ARL_c_I_von", "ARL_c_I_bis",
+ "BL_ANAT_hydromorph")) %>% #empty trait
mutate(species0=species) %>%
rowwise() %>%
# quick and dirty clean up names
@@ -27,13 +28,17 @@ dim(traits0) #907 obs. of 75 variables:
traits0 <- traits0 %>%
dplyr::select_if(~mean(!is.na(.)) >= 0.88) # 907 x 67
+
+
+
+
### Transform binary traits to 0-1
traits.asym.binary <- c('LEB_F_Makrophanerophyt','LEB_F_Nanophanerophyt',
'LEB_F_Hemikryptophyt','LEB_F_Geophyt','LEB_F_Hemiphanerophyt','LEB_F_Therophyt',
'LEB_F_Hydrophyt','LEB_F_Pseudophanerophyt','LEB_F_Chamaephyt',
'LEB_D_plurienn_pollakanth','LBE_D_plurienn_hapaxanth','LEB_D_annuell',
- 'LEB_D_bienn','V_VER_absent','V_VER_Wurzelspross','V_VER_Ausläufer','
- V_VER_Rhizom','V_VER_Innovationsknopse.mit.Wurzelknolle',
+ 'LEB_D_bienn','V_VER_absent','V_VER_Wurzelspross','V_VER_Ausläufer',
+ 'V_VER_Rhizom','V_VER_Innovationsknopse.mit.Wurzelknolle',
'V_VER_Innovationsknospe.mit.Speicherwurzel','V_VER_Ausläuferknolle',
'V_VER_Brutsprösschen','V_VER_Fragmentation','V_VER_Turio','V_VER_Sprossknolle',
'V_VER_phyllogener_Spross','V_VER_Rhizompleiokorm','V_VER_Zwiebel',
@@ -87,6 +92,8 @@ env <- env0 %>%
filter(!is.na(LAT)) %>%
filter(!(LAT==0 | LON==0))
+env.all <- env
+
### 3. Import species data ####
@@ -160,25 +167,73 @@ traits <- traits %>%
traits <- traits %>%
as.tbl() %>%
- dplyr::select(-starts_with("BL_FORM"), -starts_with("REPR_T"), -starts_with("BLU_KL"), -starts_with("STRAT_T")) %>%
+ dplyr::select(-starts_with("BL_FORM"), -starts_with("REPR_T"), -starts_with("BLU_KL"), -starts_with("STRAT_T"), -starts_with("BL_AUSD")) %>%
left_join(traits %>%
- dplyr::select(species0, REPR_T_Samen_Sporen:STRAT_T_SR) %>%
+ dplyr::select(species0, `BL_AUSD_immergrün`:`BL_AUSD_überwinternd_grün`, REPR_T_Samen_Sporen:STRAT_T_SR) %>%
gather(key=Trait, value="value", -species0) %>%
separate(Trait, into = c("Trait", "Organ", "Level"), sep = "_", extra = "merge") %>%
unite(Trait, Trait, Organ) %>%
filter(value==1) %>%
dplyr::select(-value) %>%
spread(Trait, Level) %>%
- mutate_at(.vars=vars(BLU_KL:STRAT_T),
+ mutate_at(.vars=vars(BL_AUSD:STRAT_T),
.funs=~as.factor(.)),
by="species0")
+## recode traits to numeric
+robust.mean <- function(x1,x2=NA,x3=NA,x4=NA){
+ x <- c(x1,x2,x3,x4)
+ if(any(!is.na(x))){mean(x, na.rm=T)} else {NA}
+}
+
+traits <- traits %>%
+ dplyr::select(-starts_with("BL_ANAT"), -starts_with("LEB_D"), -starts_with("ROS_T")) %>%
+ left_join(traits %>%
+ dplyr::select(species0, starts_with("BL_ANAT")) %>%
+ mutate(BL_ANAT_helomorph=ifelse(BL_ANAT_helomorph==1, 1, NA)) %>%
+ mutate(BL_ANAT_hygromorph=ifelse(BL_ANAT_hygromorph==1, 2, NA)) %>%
+ mutate(BL_ANAT_mesomorph=ifelse(BL_ANAT_mesomorph==1, 3, NA)) %>%
+ mutate(BL_ANAT_skleromorph=ifelse(BL_ANAT_skleromorph==1, 4, NA)) %>%
+ rowwise() %>%
+ mutate(BL_ANAT=robust.mean(BL_ANAT_helomorph, BL_ANAT_hygromorph, BL_ANAT_mesomorph, BL_ANAT_skleromorph)) %>%
+ ungroup() %>%
+ dplyr::select(species0, BL_ANAT, BL_ANAT_blattsukkulent),
+ by="species0") %>%
+ left_join(traits %>%
+ dplyr::select(species0, starts_with("LEB_D")) %>%
+ rowwise() %>%
+ mutate(LEB_D_plurienn=max(LEB_D_plurienn_pollakanth + LEB_D_plurienn_hapaxanth, na.rm=T)) %>%
+ ungroup() %>%
+ mutate(LEB_D_plurienn=ifelse(LEB_D_plurienn==1, 3, NA)) %>%
+ mutate(LEB_D_annuell=ifelse(LEB_D_annuell==1, 1, NA)) %>%
+ mutate(LEB_D_bienn =ifelse(LEB_D_bienn==1, 2, NA)) %>%
+ rowwise() %>%
+ mutate(LEB_D=robust.mean(LEB_D_annuell, LEB_D_bienn, LEB_D_plurienn)) %>%
+ ungroup() %>%
+ dplyr::select(species0, LEB_D),
+ by="species0") %>%
+ left_join(traits %>%
+ dplyr::select(species0, starts_with("ROS_T")) %>%
+ mutate(ROS_T=ROS_T_Ganzrosettenpflanzen) %>%
+ mutate(ROS_T=replace(ROS_T,
+ list=ROS_T_Halbrosettenpflanze==1,
+ values=0.5)) %>%
+ mutate(ROS_T=replace(ROS_T,
+ list=ROS_T_rosettenlose.Pflanzen==1,
+ values=0)) %>%
+ dplyr::select(species0, ROS_T),
+ by="species0")
+
+### ordered factors
+
dim(species) #558 783
-dim(traits) #783 81
+dim(traits) #783 53
dim(env) #558 8
+
+
## Select only plots where >90% of species have trait info [TRY]
# releve08trait <- species %>%
# rownames_to_column("RELEVE_NR") %>%
@@ -203,6 +258,15 @@ dim(env) #558 8
##export for Valerio
write_delim(species, path="_data/Mesobromion/species.out.10perc.txt", delim="\t")
write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
+write_delim(env, path="_data/Mesobromion/env.10perc.txt", delim="\t")
+
+## version without missing species
+ empty <- which(colSums(species)==0)
+ traits_nozero <- traits[-empty,]
+ species_nozero <- species[,-empty]
+
+write_delim(species_nozero , path="_data/Mesobromion/species.out.10perc_nozero.txt", delim="\t")
+write_delim(traits_nozero, path="_data/Mesobromion/traits.out.10perc_nozero.txt", delim="\t")
@@ -217,21 +281,47 @@ write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
## 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")
+trait.labs <- read_delim("_data/Mesobromion/TraitLabels.csv", delim=",") %>%
+ rownames_to_column("Trait.comb")
traits <- traits %>%
+ as.data.frame() %>%
# filter(species0 %in% colnames(species)) %>%
+ mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0")
-
+
+#rename trait based on short labels
+colnames(traits) <- trait.labs$Short_english_name
+
+
+##check traits
+is.binary <- function(x){(all(na.omit(x) %in% 0:1))}
+trait.recap <- traits %>%
+ mutate_if(.predicate = ~is.numeric(.), ~round(.,3)) %>%
+ summarize_all(.funs = list(xxxType.of.variable=~ifelse(is.binary(.), "binary",
+ ifelse(is.ordered(.), "ordered",
+ ifelse(is.numeric(.), "numeric",
+ ifelse(is.factor(.), "factor", NA)))),
+ xxxLevels=~(
+ ifelse(is.binary(.),
+ paste(names(table(.)), "=", table(.), collapse="; "),
+ ifelse(is.numeric(.),
+ paste(range(., na.rm=T), collapse=" - "),
+ ifelse(is.factor(.),
+ paste(levels(.), collapse=", "),
+ NA)))))) %>%
+ gather(key="Variable") %>%
+ separate(Variable, into = c("Variable", "feature"), sep="_xxx") %>%
+ spread(key=feature, value = value) %>%
+ rename(`Range/Levels`=Levels) %>%
+ mutate(Variable=factor(Variable, levels=colnames(traits))) %>%
+ arrange(Variable)
+
+## Shoudln't BL_ANAT better treated as ordered.factor?
+write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
#dplyr::select(LeafArea:Disp.unit.leng)
#dplyr::select(PlantHeight, LeafC.perdrymass, LeafN, StemDens, Stem.cond.dens, Seed.num.rep.unit, SLA)
-comm <- species
-traits <- traits
-trait.sel <- 6
-a <- get.corXY.bootstrap(comm=species, traits=traits,
- trait.sel=6, bootstrap=10)
-
-
@@ -240,44 +330,51 @@ a <- get.corXY.bootstrap(comm=species, traits=traits,
#### ## Import output ####
myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_round_[0-9]+.RData", full.names = T)
#dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
-load("_derived/Mesobromion/HIDDEN_round_7.RData")
-
-
-corXY = bind_rows(dataFiles) %>%
- as.tbl()
-rm( dataFiles)
-
-trait.labs <- data.frame(Trait.comb=1:ncol(traits),
- trait.name=colnames(traits))#[-which(colnames(traits) %in% c("species", "species0"))])
-
-## calculate confidence intervals for bootstrapped correlations
-corXY.ci <- corXY %>%
- filter(bootstr.n==0) %>%
- dplyr::select(-bootstr.n) %>%
- group_by(Trait.comb) %>%
- slice(1) %>%
- left_join(corXY %>%
- filter(bootstr.n!=0) %>%
- group_by(Trait.comb) %>%
- summarize(q025=quantile(Coef.obs, 0.025),
- q975=quantile(Coef.obs, 0.975),
- greater.than.perm=mean(Coef.obs>Coef.perm),
- n=n())) %>%
- #calculate significance using permuted correlations
- mutate(sign_plus=greater.than.perm>0.995) %>%
- mutate(sign_minus=greater.than.perm<0.005) %>%
- rowwise() %>%
- mutate(ntraits= length(unlist(strsplit(Trait.comb, "_")))) %>%
- ungroup() %>%
- mutate_at(.vars=vars(starts_with("sign")),
- .funs=~factor(.*1, levels=c(0,1), labels=c("FALSE", "TRUE")))
+load("_derived/Mesobromion/HIDDEN_round_11.RData")
+
+
+
+
+#corXY = bind_rows(dataFiles) %>%
+# as.tbl()
+#rm( dataFiles)
+
+#trait.labs <- data.frame(Trait.comb=1:ncol(traits),
+# trait.name=colnames(traits)) %>% #[-which(colnames(traits) %in% c("species", "species0"))])
+# # shorten some labels
+# mutate(trait.name = fct_recode(trait.name,
+# V_VER_Wurzelknolle = "V_VER_Innovationsknopse.mit.Wurzelknolle",
+# V_VER_Speicherwurzel = "V_VER_Innovationsknospe.mit.Speicherwurzel"))
+#
+
+# ## calculate confidence intervals for bootstrapped correlations
+# corXY.ci <- corXY %>%
+# filter(bootstr.n==0) %>%
+# dplyr::select(-bootstr.n) %>%
+# group_by(Trait.comb) %>%
+# slice(1) %>%
+# left_join(corXY %>%
+# filter(bootstr.n!=0) %>%
+# group_by(Trait.comb) %>%
+# summarize(q025=quantile(Coef.obs, 0.025),
+# q975=quantile(Coef.obs, 0.975),
+# greater.than.perm=mean(Coef.obs>Coef.perm),
+# n=n())) %>%
+# #calculate significance using permuted correlations
+# mutate(sign_plus=greater.than.perm>0.995) %>%
+# mutate(sign_minus=greater.than.perm<0.005) %>%
+# rowwise() %>%
+# mutate(ntraits= length(unlist(strsplit(Trait.comb, "_")))) %>%
+# ungroup() %>%
+# mutate_at(.vars=vars(starts_with("sign")),
+# .funs=~factor(.*1, levels=c(0,1), labels=c("FALSE", "TRUE")))
corXY.ci <- corXY.ci %>%
mutate(Trait.comb2=Trait.comb) %>%
#split strings of trait combinations and add labels
- separate(Trait.comb2, into=paste0("trait", 1:22)) %>%
- mutate_at(.vars=vars(trait1:trait22),
- .funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$trait.name)) %>%
+ separate(Trait.comb2, into=paste0("trait", 1:11)) %>%
+ mutate_at(.vars=vars(trait1:trait11),
+ .funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$Short_english_name)) %>%
arrange(ntraits, Coef.obs) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything())
@@ -302,24 +399,24 @@ mydata.byone <- mydata %>%
mutate(seq=1:n())
-top.first <- ggplot(data=mydata.byone) +
+(top.first <- ggplot(data=mydata.byone %>%
+ mutate(sign_plus=fct_rev(as.factor(sign_plus)))) +
geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq,
col=sign_plus)) +
geom_point(aes(x=Coef.obs, y=seq), pch=15) +
scale_y_continuous(breaks=mydata.byone$seq,
labels=mydata.byone$trait1, name=NULL) +
scale_x_continuous(name="RD correlation") +
- labs(fill = "Significant")
-
-tt2 <- ttheme_minimal()
-topfirst.leg <- cowplot::plot_grid(top.first,
- tableGrob(trait.labs %>%
- dplyr::select(trait.name), theme=tt2),
- nrow=1, rel_widths=c(0.65, 0.35))
+ labs(color = "Significant") +
+ theme_bw() +
+ theme(panel.grid.minor = element_blank(),
+ axis.text = element_text(size=7),
+ legend.position = c(0.8, 0.1))
+)
ggsave(filename = "_pics/TopFirstPredictors_CI.png", dpi=400,
- width=6, height=5, topfirst.leg)
+ width=4, height=6, plot = top.first)
allpredictors <- top.first %+% (mydata) +
scale_y_continuous(breaks=mydata$seq,
@@ -327,11 +424,11 @@ allpredictors <- top.first %+% (mydata) +
# labels=mydata$Trait.comb) +
theme_classic() +
- theme(axis.text.y = element_text(size=3)) #+
- #facet_wrap(.~ntraits, scales = "free_y", nrow=3)
+ theme(axis.text.y = element_text(size=3))# +
+ # facet_wrap(.~ntraits, scales = "free_y", nrow=3)
ggsave(filename = "_pics/All_predictors_sign_individually_CI_faceted_0506.png", dpi=400,
- width=6, height=30, allpredictors)
+ width=10, height=10, allpredictors)
@@ -351,18 +448,19 @@ get.best <- function(x, N){
top.one.by.one <- get.best(mydata, N=1)
-maxtraits <- 7
+maxtraits <- 11
for(nn in 1:maxtraits){
if(nn==1) {
best.at.1 <- get.best(mydata, N=nn)
newdata <- mydata %>%
- filter_at(.vars=vars(trait1:trait15),
+ filter_at(.vars=vars(trait1:trait11),
.vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.)))
new.best.row <- newdata %>%
filter(Trait.comb==best.at.1$Trait.comb)
upper <- new.best.row$q975
lower <- new.best.row$q025
print(paste("new best at nn", nn, best.at.1$trait.name))
+ best.progr <- best.at.1$Trait.comb
}
if(nn>1){
better <- list()
@@ -392,28 +490,69 @@ for(nn in 1:maxtraits){
lower <- new.best.row$q025
print(paste("new best at nn", nn, paste(better$trait.name, collapse=" ")))
best <- better
+ best.progr <- c(best.progr, better$Trait.comb)
}
}
}
-
-
-
mydata.best <- mydata %>%
- filter_at(.vars=vars(trait1:trait15),
- .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>%
- filter(ntraits <= length(best$trait.name)) %>%
+ filter_at(.vars=vars(trait1:trait11),
+ # .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>%
+ .vars_predicate = all_vars(. %in% traits.sign.alone | is.na(.))) %>%
+ #filter(ntraits <= length(best$trait.name)) %>%
+ filter(sign_plus==T) %>%
arrange(ntraits, Coef.obs) %>%
- mutate(seq=1:n())
+ group_by(ntraits) %>%
+ slice(n()) %>%
+ ungroup() %>%
+ bind_rows(mydata %>%
+ filter(ntraits==1) %>%
+ filter(trait1 %in% traits.sign.alone) %>%
+ arrange(Coef.obs) %>%
+ slice(1:n()-1)) %>%
+ arrange(ntraits, Coef.obs) %>%
+ mutate(seq=1:n()) %>%
+ mutate(sign_plus=factor(Trait.comb %in% best.progr))
-(top.all <- top.first %+% mydata.best +
- scale_y_continuous(breaks=mydata.best$seq,
- labels=mydata.best$Trait.comb) +
- theme(axis.text.y = element_text(size=3)))
+write_csv(mydata.best %>%
+ dplyr::select(Trait.comb:sign_plus),
+ path = "_derived/Mesobromion/BestSolutionTiers.csv")
-ggsave(filename = "_pics/Best8_AllCombinations_CI.png", dpi=400,
- width=6, height=10, top.all)
+(top.all <- ggplot(data=(mydata.best %>%
+ mutate(size0=.6+(as.numeric(sign_plus)-1)*.6))) +
+ geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, col="a",
+ lwd=size0)) +
+ geom_point(aes(x=Coef.obs, y=seq), pch=15) +
+ scale_y_continuous(breaks=mydata.best$seq,
+ labels=mydata.best$Trait.comb, name=NULL) +
+ scale_x_continuous(name="RD correlation") +
+ scale_size_identity() +
+ theme_bw() +
+ theme(panel.grid.minor = element_blank(),
+ axis.text = element_text(size=7),
+ legend.position = "none"))
+
+
+tt2 <- ttheme_minimal(
+ core = list(fg_params=list(cex = .7),
+ padding=unit(c(1, 1), "mm")),
+ colhead = list(fg_params=list(cex = .7)),
+ rowhead = list(fg_params=list(cex = .7)))
+
+(topall.leg <- cowplot::plot_grid(top.all,
+ tableGrob(trait.labs %>%
+ mutate(Code=1:n()) %>%
+ dplyr::select(Code, Trait=Short_english_name) %>%
+ filter(Trait %in% traits.sign.alone),
+ theme=tt2, rows = NULL),
+ #tableGrob(trait.labs %>%
+ # mutate(Code=1:n()) %>%
+ # dplyr::select(Code, Trait=trait.name) %>%
+ # slice(33:n()), theme=tt2),
+ nrow=1, rel_widths=c(0.70, 0.3)))
+ggsave(filename = "_pics/Best_AllCombinations_CI.png", dpi=400,
+ width=6, height=3, topall.leg)
mydata.out <- mydata %>%
#filter(ntraits<3)%>%
@@ -458,7 +597,7 @@ break()
-#### Calculate and plot PCA of CWMs
+##### Calculate CWM and plot PCA ####
CWM.wide <- species %>%
rownames_to_column("RELEVE_NR") %>%
reshape2::melt(.id="RELEVE_NR") %>%
@@ -471,11 +610,62 @@ CWM.wide <- species %>%
rownames_to_column("species0"),
by="species0") %>%
group_by(RELEVE_NR) %>%
- summarize_at(.vars = vars(LEB_F_Makrophanerophyt:Disp.unit.leng),
+ select_if(.predicate=~is.numeric(.)) %>% # CWM calculated ONLY for numeric traits
+ summarize_at(.vars = vars(any_of(traits.sign.alone)), #GF_Macrophan:Rosette),
.funs = list(~weighted.mean(., pres, na.rm=T))) %>%
dplyr::select(RELEVE_NR, order(colnames(.))) %>%
column_to_rownames("RELEVE_NR")
+### PCA of Y (Bealls) matrix +
+library(vegan)
+W.beals <- as.data.frame(beals(species, include=T, type=2))
+pca.out <- rda(W.beals)
+varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
+cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:3)
+
+
+PCA1_2 <- ggplot() +
+ geom_point(data=as.data.frame(pca.out$CA$u[,1:2]),
+ aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
+ aes(x=0, xend=PC1, y=0, yend=PC2), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
+ rownames_to_column("Trait"),
+ aes(x=PC1, y=PC2, label=Trait), col="Dark blue", size=2,
+ position = position_dodge(1) ) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ scale_y_continuous(limits=c(-0.25, 0.25)) +
+ scale_x_continuous(limits=c(-0.25, 0.25)) + coord_equal() +
+ theme(panel.grid = element_blank())
+
+(PCA1_3 <- ggplot() +
+ geom_point(data=as.data.frame(pca.out$CA$u[,1:3]),
+ aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
+ aes(x=0, xend=PC1, y=0, yend=PC3), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
+ rownames_to_column("Trait"),
+ aes(x=PC1, y=PC3, label=Trait), col="Dark blue", size=2,
+ position = position_dodge(1) ) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC3 (", varexpl[3], "%)", sep="")) +
+ theme_bw() +
+ scale_y_continuous(limits=c(-0.25, 0.25)) +
+ scale_x_continuous(limits=c(-0.25, 0.25)) +
+ coord_equal() +
+ theme(panel.grid = element_blank())
+)
+PC_beals <- cowplot::plot_grid(PCA1_2,PCA1_3, nrow=1)
+
+ggsave("_pics/PC_Beals.png", width=10, height=5, dpi=300, last_plot())
+
+
+
+#### Old down from here
+
+
#PCA of CWMs
CWM.pca <- vegan::rda(CWM.wide, scale=T)
varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
@@ -524,44 +714,13 @@ ggplot() +
-#### Map of plots
-library(rgdal)
-library(sp)
-library(sf)
-DE_NUTS1 <- readOGR(dsn="/data/sPlot/users/Francesco/Project_11/Germany/_data/Spatial", layer="Germany_LVL1_WGS84")
-DE_NUTS_sf <- DE_NUTS1 %>%
- st_as_sf()
-
-env.sf <- SpatialPointsDataFrame(coords=env %>% dplyr::select(LON, LAT),
- data=env %>% dplyr::select(-LON, -LAT),
- proj4string = raster::crs(DE_NUTS1)) %>%
- st_as_sf()
-
-env.all.sf <- SpatialPointsDataFrame(coords=env.all %>% dplyr::select(LON, LAT),
- data=env.all %>% dplyr::select(-LON, -LAT),
- proj4string = raster::crs(DE_NUTS1)) %>%
- st_as_sf()
-(basemap2 <- ggplot(data=DE_NUTS_sf) +
- geom_sf(col="darkgrey", lwd=0.5, fill="white", alpha=0.5)+
- geom_sf(data=env.all.sf, col=gray(0.5), size=1) +
- geom_sf(data=env.sf, col="red", alpha=0.5, size=1) +
- #geom_polygon(col="darkgrey", lwd=0.5, fill="white")+
- coord_sf() +
- theme_bw()
-)
-
-ggsave(filename="_data/Mesobromion/PlotDistribution.png", height=8, width=6, dpi=300, basemap2)
-
-
-
-## Graph of PCoA + passive projection of CWM
library(vegan)
fuzzy <- read.table("_data/Mesobromion/X_Dry-Grassland_581plots_488spp_R.txt", header=T)#, delim="\t")
cwms <- read.table("_data/Mesobromion/X_scores_T_Dry-Grassland_581plots_5traits_Rao_R.txt", header=T) %>%
dplyr::select(-Axis1, -Axis2) #, delim="\t")
-
+
pcoa.out <- wcmdscale(dist(fuzzy), eig=T)
varexpl <- round((pcoa.out$eig)/sum(pcoa.out$eig)*100,1)
cwms.envfit <- envfit(pcoa.out, cwms)
@@ -576,8 +735,8 @@ ggplot() +
geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.0008),
aes(x=0, xend=Dim1, y=0, yend=Dim2), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_text(data=as.data.frame(cwms.envfit$vectors$arrows*.00092) %>%
- rownames_to_column("Trait"),
- aes(x=Dim1, y=Dim2, label=Trait), col="Dark blue", size=3 ) +
+ rownames_to_column("Trait"),
+ aes(x=Dim1, y=Dim2, label=Trait), col="Dark blue", size=3 ) +
xlab(paste("PCoA1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PCoA2 (", varexpl[2], "%)", sep="")) +
theme_bw() +
@@ -589,3 +748,109 @@ ggsave("_data/Mesobromion/PCoA_traits.png", width=6, height=4, dpi=300, last_plo
+
+
+
+
+#### Map of plots ####
+library(rgdal)
+library(sp)
+library(sf)
+env <- read_delim(file = "_data/Mesobromion/env.10perc.txt", delim = "\t")
+DE_NUTS1 <- readOGR(dsn="/data/sPlot/users/Francesco/Project_11/Germany/_data/Spatial", layer="Germany_LVL1_WGS84")
+DE_NUTS_sf <- DE_NUTS1 %>%
+ st_as_sf()
+
+env.sf <- SpatialPointsDataFrame(coords=env %>% dplyr::select(LON, LAT),
+ data=env %>% dplyr::select(-LON, -LAT),
+ proj4string = raster::crs(DE_NUTS1)) %>%
+ st_as_sf()
+
+env.all.sf <- SpatialPointsDataFrame(coords=env.all %>% dplyr::select(LON, LAT),
+ data=env.all %>% dplyr::select(-LON, -LAT),
+ proj4string = raster::crs(DE_NUTS1)) %>%
+ st_as_sf()
+
+
+(basemap2 <- ggplot(data=DE_NUTS_sf) +
+ geom_sf(col="darkgrey", lwd=0.5, fill="white", alpha=0.5)+
+ geom_sf(data=env.all.sf, col=gray(0.5), size=1) +
+ geom_sf(data=env.sf, col="red", alpha=0.5, size=1) +
+ #geom_polygon(col="darkgrey", lwd=0.5, fill="white")+
+ coord_sf() +
+ theme_bw()
+)
+
+ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basemap2)
+
+
+
+
+
+
+### renaming variables
+"LF_Macroph = LEB_F_Makrophanerophyt
+LF_Nanoh = LEB_F_Nanophanerophyt
+LF_Hemicr = LEB_F_Hemikryptophyt
+LF_Geoph = LEB_F_Geophyt
+LF_Hemiph = LEB_F_Hemiphanerophyt
+LF_Theropf = LEB_F_Therophyt
+LF_Hydroph = LEB_F_Hydrophyt
+LF_Pseudoph = LEB_F_Pseudophanerophyt
+LF_Chamaeph = LEB_F_Chamaephyt
+LL_Pluri_poll = LEB_D_plurienn_pollakanth
+LL_Pluri_hapa = LEB_D_plurienn_hapaxanth
+LL_Annual = LEB_D_annuell
+LL_Bienn = LEB_D_bienn
+VS_Absent = V_V_VER_absent
+VS_RootSprout = V_VER_Wurzelspross
+VS_Runners = V_VER_Ausläufer
+VS_Rhizon = V_VER_Rhizom
+VS_Tuber = V_VER_Innovationsknopse.mit.Wurzelknolle
+VS_StorageRoot = V_VER_Innovationsknospe.mit.Speicherwurzel
+VS_RunningTuber = V_VER_Ausläuferknolle
+VS_Brutsprösschen = V_VER_Brutsprösschen ## How to translate?
+VS_Fragment = V_VER_Fragmentation
+VS_Turio = V_VER_Turio
+24 V_VER_Sprossknolle
+25 V_VER_phyllogener_Spross
+26 V_VER_Rhizompleiokorm
+27 V_VER_Zwiebel
+28 V_VER_Ausläuferrhizom
+29 V_VER_Ausläuferzwiebel
+30 V_VER_Bulbille
+31 ROS_T_rosettenlose.Pflanzen
+32 ROS_T_Halbrosettenpflanze
+33 ROS_T_Ganzrosettenpflanzen
+34 BL_AUSD_immergrün
+35 BL_AUSD_sommergrün
+36 BL_AUSD_vorsommergrün
+37 BL_AUSD_überwinternd_grün
+38 BL_ANAT_skleromorph
+39 BL_ANAT_mesomorph
+40 BL_ANAT_hygromorph
+41 BL_ANAT_hydromorph
+42 BL_ANAT_blattsukkulent
+43 BL_ANAT_helomorph
+44 BL_BEG
+45 BL_END
+46 BL_DAU
+47 LeafArea
+48 SLA
+49 LeafC.perdrymass
+50 LeafN
+51 LeafP
+52 PlantHeight
+53 SeedMass
+54 Seed.length
+55 LDMC
+56 LeafNperArea
+57 LeafNPratio
+58 Leaf.delta.15N
+59 Seed.num.rep.unit
+60 Leaffreshmass
+61 Disp.unit.leng
+62 BLU_KL
+63 REPR_T
+64 STRAT_T"
+
diff --git a/99_HIDDEN_functions.R b/99_HIDDEN_functions.R
index 1f206a113926f7e6ee9b757da38d12bb0feee9be..ef81a8006912f9ff823495935f21469d4dff945c 100644
--- a/99_HIDDEN_functions.R
+++ b/99_HIDDEN_functions.R
@@ -66,23 +66,24 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){
## caution
## ALL columns with only 0-1 values are AUTOMATICALLY considered as asym.bin sensu FD:gowdis
##get all columns with binary variables
+ ## CAUTION - function below doesn't work as expected
# binary.traits <- which(apply(traits[,ii,drop=F], MARGIN=2, function(x)( all(na.omit(x) %in% 0:1) ))==T)
syn.out.tmp <- matrix.x(comm=comm, traits=traits[,ii,drop=F], scale=T)$matrix.X #, asym.bin=binary.traits
W.beals <- as.data.frame(beals(comm, include=T, type=2))
- # permtute traits
+ # permute traits
index.traits <- lapply(1:(bootstrap+1), function(x){sample(1:n.species, replace=F)})
syn.out.perm.tmp <- matrix.x(comm=comm, traits=traits[index.traits[[bootstrap+1]],ii,drop=F],
scale=T)$matrix.X #, asym.bin=binary.traits
- corXY <- NULL
+
#RD.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.tmp), nrepet = 0)
RD.tmp <- dcor(W.beals, as.data.frame(syn.out.tmp))^2
#RV.perm.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.perm.tmp), nrepet = 0)
RD.perm.tmp <- dcor(W.beals, as.data.frame(syn.out.perm.tmp))^2
- corXY <- rbind(corXY,
- data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RD",
- Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp))
+ corXY <- NULL
+ corXY <- data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RD",
+ Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp)
index.bootstr <- lapply(1:bootstrap, function(x){sample(1:n.sites, replace=T)})
for(b in 1:bootstrap){