diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index 79ac82d5f1e298d09abcaa8b63d9e7d29692839a..124995ab24861b27c7c56b1c16a4619dc394fbe0 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -183,13 +183,90 @@ for(ff in myfilelist){
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"))])
+aa <- data.frame(Trait.comb=paste0("t", 1:80),
+ 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") %>%
+ separate(Trait.comb, c("trait1", "trait2")) %>%
+ mutate(trait2=paste0("t", trait2)) %>%
+ 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) %>%
+ dplyr::select(trait1, trait2, q025:conf.p) %>%
arrange(desc(SES.np))
print(bb, n=20)
+#### Calculate and plot PCA of CWMs
+CWM.wide <- species %>%
+ rownames_to_column("RELEVE_NR") %>%
+ reshape2::melt(.id="RELEVE_NR") %>%
+ rename(species0=variable, pres=value) %>%
+ as.tbl() %>%
+ filter(pres>0) %>%
+ arrange(RELEVE_NR) %>%
+ ## attach traits
+ left_join(traits %>%
+ rownames_to_column("species0"),
+ by="species0") %>%
+ group_by(RELEVE_NR) %>%
+ summarize_at(.vars = vars(LEB_F_Makrophanerophyt:Disp.unit.leng),
+ .funs = list(~weighted.mean(., pres, na.rm=T))) %>%
+ dplyr::select(RELEVE_NR, order(colnames(.))) %>%
+ column_to_rownames("RELEVE_NR")
+
+#PCA of CWMs
+CWM.pca <- vegan::rda(CWM.wide, scale=T)
+varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
+
+ggplot() +
+ geom_point(data=as.data.frame(CWM.pca$CA$u[,1:2]),
+ aes(x=PC1, y=PC2), 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=PC2, 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=PC2, label=Trait, col=in.try), size=3 ) +
+ scale_color_identity() +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank()) +
+ title("PCA of CWMs")
+
+#PCA of individual traits
+trait.pca <- vegan::rda(as.matrix(traits %>%
+ 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*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*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[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank()) +
+ title("PCA of species-level traits")
+
+
+
+
#### Map of plots
library(rgdal)
library(sp)
diff --git a/98_SummarizeSimulations.R b/98_SummarizeSimulations.R
index 45bb8bdcf6ab7fad36838adc700534650315ff06..8fe4dc38b9d95434d95105f5584fb90844b5b09e 100644
--- a/98_SummarizeSimulations.R
+++ b/98_SummarizeSimulations.R
@@ -3,7 +3,7 @@
## It then plots the summarized output
library(tidyverse)
-mypath <- "_data/Experiment_27Feb2020"
+mypath <- "_data/Experiment_29Feb2020"
myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T)
output <- NULL
@@ -27,7 +27,7 @@ for(ff in myfiles){
}
outp.summary <- output %>%
- filter(!stat.type %in% c("XY", "XY.T", "XY.TR")) %>%
+ dplyr::filter(!stat.type %in% c("XY", "XY.T", "XY.TR")) %>%
group_by(main, inter, feedb, trait, envir, stat.type) %>%
summarize(stat.obs.med=median(stat.obs),
power=mean(pvalue<=0.05),
@@ -35,7 +35,7 @@ outp.summary <- output %>%
exp.med.med=median(exp.med),
nsim=n()) %>%
bind_rows(output %>%
- filter(stat.type %in% c("XY", "XY.T", "XY.TR")) %>%
+ dplyr::filter(stat.type %in% c("XY", "XY.T", "XY.TR")) %>%
group_by(main, inter, feedb, trait, stat.type) %>%
summarize(stat.obs.med=median(stat.obs),
power=mean(pvalue<=0.05),
@@ -46,33 +46,78 @@ outp.summary <- output %>%
arrange(stat.type, main, inter, feedb, trait, envir)
-## plotting XY
+## plotting power for XY with feedback
ggplot(data=outp.summary %>%
ungroup() %>%
- filter(stat.type=="XY") %>%
- filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ dplyr::filter(stat.type=="XY") %>%
+ filter(feedb==0) %>%
+ #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ dplyr::filter(trait %in% c("1", "2", "3")) %>%
mutate(inter=as.factor(inter))) +
geom_line(aes(x=main, y=power, group=trait, col=trait)) +
scale_colour_brewer(palette = "Dark2") +
- facet_grid(feedb~inter) +
+ facet_grid(.~inter) +
theme_bw() +
theme(panel.grid = element_blank())
+## plotting non-parametric SES for XY with feedback
+
+ggplot(data=outp.summary %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY") %>%
+ filter(feedb==0) %>%
+ #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ dplyr::filter(trait %in% c("1", "2", "3")) %>%
+ mutate(inter=as.factor(inter))) +
+ geom_line(aes(x=main, y=SES.med, group=trait, col=trait)) +
+ scale_colour_brewer(palette = "Dark2") +
+ facet_grid(.~inter) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+
+
## plotting XE
ggplot(data=outp.summary %>%
ungroup() %>%
- filter(stat.type=="XE") %>%
- filter(envir=="1 ") %>%
+ dplyr::filter(stat.type=="XE") %>%
+ #dplyr::filter(envir=="1") %>%
+ filter(feedb==0) %>%
+ filter(trait %in% c("1","2","3", "1 2")) %>%
mutate(inter=as.factor(inter))
#mutate(group0=paste("t", trait, " - e", envir))
) +
- geom_line(aes(x=main, y=power, group=trait, col=trait), lwd=1.5, alpha=0.3) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
+ scale_colour_brewer(palette = "Dark2") +
+ facet_grid(envir~inter) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+
+## plotting XY.T
+
+ggplot(data=outp.summary %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY.T") %>%
+ dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ mutate(inter=as.factor(inter))) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
scale_colour_brewer(palette = "Dark2") +
facet_grid(feedb~inter) +
theme_bw() +
theme(panel.grid = element_blank())
+## plotting XY.TR
+
+ggplot(data=outp.summary %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY.TR") %>%
+ dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ mutate(inter=as.factor(inter))) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
+ scale_colour_brewer(palette = "Dark2") +
+ facet_grid(feedb~inter) +
+ theme_bw() +
+ theme(panel.grid = element_blank())