Skip to content
Snippets Groups Projects
Commit 8c99549d authored by Francesco Sabatini's avatar Francesco Sabatini
Browse files

Minor updates, and trait labelling

parent ebcba6fd
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
01_Mesobromion.R
+ 377
− 112
View file @ 8c99549d
...@@ -11,7 +11,8 @@ source("99_HIDDEN_functions.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')) %>% 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", 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) %>% mutate(species0=species) %>%
rowwise() %>% rowwise() %>%
# quick and dirty clean up names # quick and dirty clean up names
...@@ -27,13 +28,17 @@ dim(traits0) #907 obs. of 75 variables: ...@@ -27,13 +28,17 @@ dim(traits0) #907 obs. of 75 variables:
traits0 <- traits0 %>% traits0 <- traits0 %>%
dplyr::select_if(~mean(!is.na(.)) >= 0.88) # 907 x 67 dplyr::select_if(~mean(!is.na(.)) >= 0.88) # 907 x 67
### Transform binary traits to 0-1 ### Transform binary traits to 0-1
traits.asym.binary <- c('LEB_F_Makrophanerophyt','LEB_F_Nanophanerophyt', traits.asym.binary <- c('LEB_F_Makrophanerophyt','LEB_F_Nanophanerophyt',
'LEB_F_Hemikryptophyt','LEB_F_Geophyt','LEB_F_Hemiphanerophyt','LEB_F_Therophyt', 'LEB_F_Hemikryptophyt','LEB_F_Geophyt','LEB_F_Hemiphanerophyt','LEB_F_Therophyt',
'LEB_F_Hydrophyt','LEB_F_Pseudophanerophyt','LEB_F_Chamaephyt', 'LEB_F_Hydrophyt','LEB_F_Pseudophanerophyt','LEB_F_Chamaephyt',
'LEB_D_plurienn_pollakanth','LBE_D_plurienn_hapaxanth','LEB_D_annuell', 'LEB_D_plurienn_pollakanth','LBE_D_plurienn_hapaxanth','LEB_D_annuell',
'LEB_D_bienn','V_VER_absent','V_VER_Wurzelspross','V_VER_Ausläufer',' 'LEB_D_bienn','V_VER_absent','V_VER_Wurzelspross','V_VER_Ausläufer',
V_VER_Rhizom','V_VER_Innovationsknopse.mit.Wurzelknolle', 'V_VER_Rhizom','V_VER_Innovationsknopse.mit.Wurzelknolle',
'V_VER_Innovationsknospe.mit.Speicherwurzel','V_VER_Ausläuferknolle', 'V_VER_Innovationsknospe.mit.Speicherwurzel','V_VER_Ausläuferknolle',
'V_VER_Brutsprösschen','V_VER_Fragmentation','V_VER_Turio','V_VER_Sprossknolle', 'V_VER_Brutsprösschen','V_VER_Fragmentation','V_VER_Turio','V_VER_Sprossknolle',
'V_VER_phyllogener_Spross','V_VER_Rhizompleiokorm','V_VER_Zwiebel', 'V_VER_phyllogener_Spross','V_VER_Rhizompleiokorm','V_VER_Zwiebel',
...@@ -87,6 +92,8 @@ env <- env0 %>% ...@@ -87,6 +92,8 @@ env <- env0 %>%
filter(!is.na(LAT)) %>% filter(!is.na(LAT)) %>%
filter(!(LAT==0 | LON==0)) filter(!(LAT==0 | LON==0))
env.all <- env
### 3. Import species data #### ### 3. Import species data ####
...@@ -160,25 +167,73 @@ traits <- traits %>% ...@@ -160,25 +167,73 @@ traits <- traits %>%
traits <- traits %>% traits <- traits %>%
as.tbl() %>% 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 %>% 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) %>% gather(key=Trait, value="value", -species0) %>%
separate(Trait, into = c("Trait", "Organ", "Level"), sep = "_", extra = "merge") %>% separate(Trait, into = c("Trait", "Organ", "Level"), sep = "_", extra = "merge") %>%
unite(Trait, Trait, Organ) %>% unite(Trait, Trait, Organ) %>%
filter(value==1) %>% filter(value==1) %>%
dplyr::select(-value) %>% dplyr::select(-value) %>%
spread(Trait, Level) %>% spread(Trait, Level) %>%
mutate_at(.vars=vars(BLU_KL:STRAT_T), mutate_at(.vars=vars(BL_AUSD:STRAT_T),
.funs=~as.factor(.)), .funs=~as.factor(.)),
by="species0") 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(species) #558 783
dim(traits) #783 81 dim(traits) #783 53
dim(env) #558 8 dim(env) #558 8
## Select only plots where >90% of species have trait info [TRY] ## Select only plots where >90% of species have trait info [TRY]
# releve08trait <- species %>% # releve08trait <- species %>%
# rownames_to_column("RELEVE_NR") %>% # rownames_to_column("RELEVE_NR") %>%
...@@ -203,6 +258,15 @@ dim(env) #558 8 ...@@ -203,6 +258,15 @@ dim(env) #558 8
##export for Valerio ##export for Valerio
write_delim(species, path="_data/Mesobromion/species.out.10perc.txt", delim="\t") 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(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") ...@@ -217,21 +281,47 @@ write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
## calculate corr between species composition matrix and traits ## calculate corr between species composition matrix and traits
species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t") species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
traits <- read_delim("_data/Mesobromion/traits.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 %>% traits <- traits %>%
as.data.frame() %>%
# filter(species0 %in% colnames(species)) %>% # filter(species0 %in% colnames(species)) %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0") 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(LeafArea:Disp.unit.leng)
#dplyr::select(PlantHeight, LeafC.perdrymass, LeafN, StemDens, Stem.cond.dens, Seed.num.rep.unit, SLA) #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, ...@@ -240,44 +330,51 @@ a <- get.corXY.bootstrap(comm=species, traits=traits,
#### ## Import output #### #### ## Import output ####
myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_round_[0-9]+.RData", full.names = T) myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_round_[0-9]+.RData", full.names = T)
#dataFiles = purrr::map(myfilelist, function(x){get(load(x))}) #dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
load("_derived/Mesobromion/HIDDEN_round_7.RData") load("_derived/Mesobromion/HIDDEN_round_11.RData")
corXY = bind_rows(dataFiles) %>%
as.tbl()
rm( dataFiles) #corXY = bind_rows(dataFiles) %>%
# as.tbl()
trait.labs <- data.frame(Trait.comb=1:ncol(traits), #rm( dataFiles)
trait.name=colnames(traits))#[-which(colnames(traits) %in% c("species", "species0"))])
#trait.labs <- data.frame(Trait.comb=1:ncol(traits),
## calculate confidence intervals for bootstrapped correlations # trait.name=colnames(traits)) %>% #[-which(colnames(traits) %in% c("species", "species0"))])
corXY.ci <- corXY %>% # # shorten some labels
filter(bootstr.n==0) %>% # mutate(trait.name = fct_recode(trait.name,
dplyr::select(-bootstr.n) %>% # V_VER_Wurzelknolle = "V_VER_Innovationsknopse.mit.Wurzelknolle",
group_by(Trait.comb) %>% # V_VER_Speicherwurzel = "V_VER_Innovationsknospe.mit.Speicherwurzel"))
slice(1) %>% #
left_join(corXY %>%
filter(bootstr.n!=0) %>% # ## calculate confidence intervals for bootstrapped correlations
group_by(Trait.comb) %>% # corXY.ci <- corXY %>%
summarize(q025=quantile(Coef.obs, 0.025), # filter(bootstr.n==0) %>%
q975=quantile(Coef.obs, 0.975), # dplyr::select(-bootstr.n) %>%
greater.than.perm=mean(Coef.obs>Coef.perm), # group_by(Trait.comb) %>%
n=n())) %>% # slice(1) %>%
#calculate significance using permuted correlations # left_join(corXY %>%
mutate(sign_plus=greater.than.perm>0.995) %>% # filter(bootstr.n!=0) %>%
mutate(sign_minus=greater.than.perm<0.005) %>% # group_by(Trait.comb) %>%
rowwise() %>% # summarize(q025=quantile(Coef.obs, 0.025),
mutate(ntraits= length(unlist(strsplit(Trait.comb, "_")))) %>% # q975=quantile(Coef.obs, 0.975),
ungroup() %>% # greater.than.perm=mean(Coef.obs>Coef.perm),
mutate_at(.vars=vars(starts_with("sign")), # n=n())) %>%
.funs=~factor(.*1, levels=c(0,1), labels=c("FALSE", "TRUE"))) # #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 %>% corXY.ci <- corXY.ci %>%
mutate(Trait.comb2=Trait.comb) %>% mutate(Trait.comb2=Trait.comb) %>%
#split strings of trait combinations and add labels #split strings of trait combinations and add labels
separate(Trait.comb2, into=paste0("trait", 1:22)) %>% separate(Trait.comb2, into=paste0("trait", 1:11)) %>%
mutate_at(.vars=vars(trait1:trait22), mutate_at(.vars=vars(trait1:trait11),
.funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$trait.name)) %>% .funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$Short_english_name)) %>%
arrange(ntraits, Coef.obs) %>% arrange(ntraits, Coef.obs) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything()) dplyr::select(Trait.comb, Test, n, ntraits, everything())
...@@ -302,24 +399,24 @@ mydata.byone <- mydata %>% ...@@ -302,24 +399,24 @@ mydata.byone <- mydata %>%
mutate(seq=1:n()) 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, geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq,
col=sign_plus)) + col=sign_plus)) +
geom_point(aes(x=Coef.obs, y=seq), pch=15) + geom_point(aes(x=Coef.obs, y=seq), pch=15) +
scale_y_continuous(breaks=mydata.byone$seq, scale_y_continuous(breaks=mydata.byone$seq,
labels=mydata.byone$trait1, name=NULL) + labels=mydata.byone$trait1, name=NULL) +
scale_x_continuous(name="RD correlation") + scale_x_continuous(name="RD correlation") +
labs(fill = "Significant") labs(color = "Significant") +
theme_bw() +
tt2 <- ttheme_minimal() theme(panel.grid.minor = element_blank(),
topfirst.leg <- cowplot::plot_grid(top.first, axis.text = element_text(size=7),
tableGrob(trait.labs %>% legend.position = c(0.8, 0.1))
dplyr::select(trait.name), theme=tt2), )
nrow=1, rel_widths=c(0.65, 0.35))
ggsave(filename = "_pics/TopFirstPredictors_CI.png", dpi=400, 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) + allpredictors <- top.first %+% (mydata) +
scale_y_continuous(breaks=mydata$seq, scale_y_continuous(breaks=mydata$seq,
...@@ -327,11 +424,11 @@ allpredictors <- top.first %+% (mydata) + ...@@ -327,11 +424,11 @@ allpredictors <- top.first %+% (mydata) +
# labels=mydata$Trait.comb) + # labels=mydata$Trait.comb) +
theme_classic() + theme_classic() +
theme(axis.text.y = element_text(size=3)) #+ theme(axis.text.y = element_text(size=3))# +
#facet_wrap(.~ntraits, scales = "free_y", nrow=3) # facet_wrap(.~ntraits, scales = "free_y", nrow=3)
ggsave(filename = "_pics/All_predictors_sign_individually_CI_faceted_0506.png", dpi=400, 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){ ...@@ -351,18 +448,19 @@ get.best <- function(x, N){
top.one.by.one <- get.best(mydata, N=1) top.one.by.one <- get.best(mydata, N=1)
maxtraits <- 7 maxtraits <- 11
for(nn in 1:maxtraits){ for(nn in 1:maxtraits){
if(nn==1) { if(nn==1) {
best.at.1 <- get.best(mydata, N=nn) best.at.1 <- get.best(mydata, N=nn)
newdata <- mydata %>% 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(.))) .vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.)))
new.best.row <- newdata %>% new.best.row <- newdata %>%
filter(Trait.comb==best.at.1$Trait.comb) filter(Trait.comb==best.at.1$Trait.comb)
upper <- new.best.row$q975 upper <- new.best.row$q975
lower <- new.best.row$q025 lower <- new.best.row$q025
print(paste("new best at nn", nn, best.at.1$trait.name)) print(paste("new best at nn", nn, best.at.1$trait.name))
best.progr <- best.at.1$Trait.comb
} }
if(nn>1){ if(nn>1){
better <- list() better <- list()
...@@ -392,28 +490,69 @@ for(nn in 1:maxtraits){ ...@@ -392,28 +490,69 @@ for(nn in 1:maxtraits){
lower <- new.best.row$q025 lower <- new.best.row$q025
print(paste("new best at nn", nn, paste(better$trait.name, collapse=" "))) print(paste("new best at nn", nn, paste(better$trait.name, collapse=" ")))
best <- better best <- better
best.progr <- c(best.progr, better$Trait.comb)
} }
} }
} }
mydata.best <- mydata %>% mydata.best <- mydata %>%
filter_at(.vars=vars(trait1:trait15), filter_at(.vars=vars(trait1:trait11),
.vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>% # .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>%
filter(ntraits <= length(best$trait.name)) %>% .vars_predicate = all_vars(. %in% traits.sign.alone | is.na(.))) %>%
#filter(ntraits <= length(best$trait.name)) %>%
filter(sign_plus==T) %>%
arrange(ntraits, Coef.obs) %>% 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 + write_csv(mydata.best %>%
scale_y_continuous(breaks=mydata.best$seq, dplyr::select(Trait.comb:sign_plus),
labels=mydata.best$Trait.comb) + path = "_derived/Mesobromion/BestSolutionTiers.csv")
theme(axis.text.y = element_text(size=3)))
ggsave(filename = "_pics/Best8_AllCombinations_CI.png", dpi=400, (top.all <- ggplot(data=(mydata.best %>%
width=6, height=10, top.all) 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 %>% mydata.out <- mydata %>%
#filter(ntraits<3)%>% #filter(ntraits<3)%>%
...@@ -458,7 +597,7 @@ break() ...@@ -458,7 +597,7 @@ break()
#### Calculate and plot PCA of CWMs ##### Calculate CWM and plot PCA ####
CWM.wide <- species %>% CWM.wide <- species %>%
rownames_to_column("RELEVE_NR") %>% rownames_to_column("RELEVE_NR") %>%
reshape2::melt(.id="RELEVE_NR") %>% reshape2::melt(.id="RELEVE_NR") %>%
...@@ -471,11 +610,62 @@ CWM.wide <- species %>% ...@@ -471,11 +610,62 @@ CWM.wide <- species %>%
rownames_to_column("species0"), rownames_to_column("species0"),
by="species0") %>% by="species0") %>%
group_by(RELEVE_NR) %>% 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))) %>% .funs = list(~weighted.mean(., pres, na.rm=T))) %>%
dplyr::select(RELEVE_NR, order(colnames(.))) %>% dplyr::select(RELEVE_NR, order(colnames(.))) %>%
column_to_rownames("RELEVE_NR") 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 #PCA of CWMs
CWM.pca <- vegan::rda(CWM.wide, scale=T) CWM.pca <- vegan::rda(CWM.wide, scale=T)
varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1) varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
...@@ -524,44 +714,13 @@ ggplot() + ...@@ -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) library(vegan)
fuzzy <- read.table("_data/Mesobromion/X_Dry-Grassland_581plots_488spp_R.txt", header=T)#, delim="\t") 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) %>% cwms <- read.table("_data/Mesobromion/X_scores_T_Dry-Grassland_581plots_5traits_Rao_R.txt", header=T) %>%
dplyr::select(-Axis1, -Axis2) #, delim="\t") dplyr::select(-Axis1, -Axis2) #, delim="\t")
pcoa.out <- wcmdscale(dist(fuzzy), eig=T) pcoa.out <- wcmdscale(dist(fuzzy), eig=T)
varexpl <- round((pcoa.out$eig)/sum(pcoa.out$eig)*100,1) varexpl <- round((pcoa.out$eig)/sum(pcoa.out$eig)*100,1)
cwms.envfit <- envfit(pcoa.out, cwms) cwms.envfit <- envfit(pcoa.out, cwms)
...@@ -576,8 +735,8 @@ ggplot() + ...@@ -576,8 +735,8 @@ ggplot() +
geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.0008), 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) + 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) %>% geom_text(data=as.data.frame(cwms.envfit$vectors$arrows*.00092) %>%
rownames_to_column("Trait"), rownames_to_column("Trait"),
aes(x=Dim1, y=Dim2, label=Trait), col="Dark blue", size=3 ) + aes(x=Dim1, y=Dim2, label=Trait), col="Dark blue", size=3 ) +
xlab(paste("PCoA1 (", varexpl[1], "%)", sep="")) + xlab(paste("PCoA1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PCoA2 (", varexpl[2], "%)", sep="")) + ylab(paste("PCoA2 (", varexpl[2], "%)", sep="")) +
theme_bw() + theme_bw() +
...@@ -589,3 +748,109 @@ ggsave("_data/Mesobromion/PCoA_traits.png", width=6, height=4, dpi=300, last_plo ...@@ -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"
99_HIDDEN_functions.R
+ 6
− 5
View file @ 8c99549d
...@@ -66,23 +66,24 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){ ...@@ -66,23 +66,24 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){
## caution ## caution
## ALL columns with only 0-1 values are AUTOMATICALLY considered as asym.bin sensu FD:gowdis ## ALL columns with only 0-1 values are AUTOMATICALLY considered as asym.bin sensu FD:gowdis
##get all columns with binary variables ##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) # 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 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)) 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)}) 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], 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 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 <- RV.rtest(W.beals, as.data.frame(syn.out.tmp), nrepet = 0)
RD.tmp <- dcor(W.beals, as.data.frame(syn.out.tmp))^2 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) #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 RD.perm.tmp <- dcor(W.beals, as.data.frame(syn.out.perm.tmp))^2
corXY <- rbind(corXY, corXY <- NULL
data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RD", corXY <- data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RD",
Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp)) Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp)
index.bootstr <- lapply(1:bootstrap, function(x){sample(1:n.sites, replace=T)}) index.bootstr <- lapply(1:bootstrap, function(x){sample(1:n.sites, replace=T)})
for(b in 1:bootstrap){ for(b in 1:bootstrap){
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment