diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index 69b1d1d7f0f667415ae20cc685c2d53ed8a94b91..99650c584a93ce1b1b9ded6164c8912e7d2d8d5e 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -140,7 +140,7 @@ species <- env %>%
#dplyr::select(traits$species0)
dim(species) # [1] 5810 873
-
+### 3.1 Subset plots based on trait coverage ####
releve08trait <- species %>%
rownames_to_column("RELEVE_NR") %>%
reshape2::melt(.id="RELEVE_NR") %>%
@@ -358,7 +358,7 @@ traits %>%
source("01b_MesobromionCluster.R")
#### 1. Traits individually significant for COVER data#### na.exclude=T ########
traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
-myfilelist <- list.files(path="_derived/Mesobromion/Cover", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T)
+myfilelist <- list.files(path="_derived/Mesobromion/Cover/onebyone", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
as_tibble() %>%
@@ -401,17 +401,17 @@ write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.cov.sign.txt",
8 35 0.241 0.128 0.217 0.289 0.970 999 TRUE FALSE 1 LeafP "
### COV - NONA2 - modified Matrix.x to deal with NAs INSIDE the analysis
-"# A tibble: 53 x 11
+"# A tibble: 48 x 11
Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 36 0.304 0.151 0.280 0.355 0.998 999 TRUE FALSE 1 PlantHeight
- 2 50 0.292 0.0641 0.265 0.346 0.988 999 TRUE FALSE 1 BL_ANAT
- 3 2 0.262 0.0633 0.223 0.312 0.989 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
- 4 30 0.259 0.145 0.230 0.307 0.976 999 TRUE FALSE 1 BL_DAU
- 5 20 0.251 0.109 0.205 0.324 0.984 999 TRUE FALSE 1 V_VER_Fragmentation
- 6 32 0.247 0.238 0.224 0.297 0.964 999 TRUE FALSE 1 SLA
- 7 35 0.241 0.200 0.221 0.287 0.965 999 TRUE FALSE 1 LeafP
- 8 49 0.238 0.0786 0.217 0.286 0.965 999 TRUE FALSE 1 STRAT_T "
+ 1 31 0.304 0.0805 0.280 0.357 0.994 999 TRUE FALSE 1 PlantHeight
+ 2 45 0.292 0.205 0.263 0.344 0.990 999 TRUE FALSE 1 BL_ANAT
+ 3 2 0.262 0.126 0.230 0.317 0.988 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 4 25 0.259 0.125 0.227 0.307 0.974 999 TRUE FALSE 1 BL_DAU
+ 5 16 0.251 0.0718 0.204 0.326 0.983 999 TRUE FALSE 1 V_VER_Fragmentation
+ 6 27 0.247 0.0761 0.227 0.297 0.965 999 TRUE FALSE 1 SLA
+ 7 30 0.241 0.207 0.213 0.286 0.951 999 TRUE FALSE 1 LeafP
+ 8 44 0.238 0.187 0.213 0.285 0.971 999 TRUE FALSE 1 STRAT_T "
### COV - without deleting NAs
"# A tibble: 53 x 11
@@ -423,7 +423,7 @@ write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.cov.sign.txt",
#### 2. Traits individually significant for Presence|absence data####
traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
-myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa-nona2_[0-9]+_.RData", full.names = T)
+myfilelist <- list.files(path="_derived/Mesobromion/PresAbs/onebyone", pattern="HIDDENpa-nona2_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
as_tibble()
@@ -466,22 +466,22 @@ write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.pa.sign.txt",
## Pres Abs NoNa2 -modified Matrix.x to deal with NAs INSIDE the analysis
## Data Filled up by HB used in this analysis
-" Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+"# A tibble: 48 x 11
+ Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 50 0.311 0.0624 0.284 0.362 0.996 999 TRUE FALSE 1 BL_ANAT
- 2 36 0.307 0.190 0.279 0.359 0.995 999 TRUE FALSE 1 PlantHeight
- 3 2 0.295 0.105 0.254 0.346 0.997 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
- 4 30 0.260 0.124 0.229 0.312 0.991 999 TRUE FALSE 1 BL_DAU
- 5 32 0.260 0.243 0.237 0.305 0.987 999 TRUE FALSE 1 SLA
- 6 35 0.253 0.131 0.230 0.299 0.983 999 TRUE FALSE 1 LeafP
- 7 33 0.240 0.0910 0.211 0.296 0.979 999 TRUE FALSE 1 LeafC.perdrymass
- 8 49 0.240 0.0954 0.212 0.287 0.975 999 TRUE FALSE 1 STRAT_T
- 9 34 0.233 0.294 0.210 0.281 0.963 999 TRUE FALSE 1 LeafN
-10 38 0.228 0.0907 0.203 0.278 0.958 999 TRUE FALSE 1 Seed.length
-11 39 0.226 0.106 0.198 0.282 0.963 999 TRUE FALSE 1 LDMC
-12 5 0.226 0.103 0.192 0.281 0.975 999 TRUE FALSE 1 LEB_F_Hemiphanerophyt
-13 1 0.202 0.0560 0.170 0.248 0.953 999 TRUE FALSE 1 LEB_F_Makrophanerophyt
-14 9 0.196 0.127 0.162 0.249 0.954 999 TRUE FALSE 1 LEB_F_Chamaephyt "
+ 1 45 0.311 0.139 0.282 0.364 0.999 999 TRUE FALSE 1 BL_ANAT
+ 2 31 0.307 0.160 0.279 0.358 0.999 999 TRUE FALSE 1 PlantHeight
+ 3 2 0.295 0.123 0.256 0.347 0.994 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 4 25 0.260 0.0982 0.231 0.307 0.983 999 TRUE FALSE 1 BL_DAU
+ 5 27 0.260 0.0963 0.237 0.305 0.984 999 TRUE FALSE 1 SLA
+ 6 30 0.253 0.158 0.229 0.300 0.987 999 TRUE FALSE 1 LeafP
+ 7 28 0.240 0.0705 0.212 0.295 0.980 999 TRUE FALSE 1 LeafC.perdrymass
+ 8 44 0.240 0.0618 0.213 0.293 0.981 999 TRUE FALSE 1 STRAT_T
+ 9 29 0.233 0.0978 0.210 0.285 0.966 999 TRUE FALSE 1 LeafN
+10 33 0.228 0.183 0.205 0.275 0.968 999 TRUE FALSE 1 Seed.length
+11 34 0.226 0.0736 0.196 0.281 0.965 999 TRUE FALSE 1 LDMC
+12 5 0.226 0.113 0.190 0.283 0.975 999 TRUE FALSE 1 LEB_F_Hemiphanerophyt
+13 40 0.209 0.128 0.185 0.259 0.957 999 TRUE FALSE 1 Disp.unit.leng "
## Pres Abs - Without excluding species with NA in traits
diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R
index 83d0d276ad9a961902b46b7150e9df95f576fd91..ab6dfe1ba13d336c08a68f0966d370fb60eac202 100644
--- a/02_Mesobromion_ExamineOutput.R
+++ b/02_Mesobromion_ExamineOutput.R
@@ -18,6 +18,7 @@ trait.labs <- data.frame(trait.name=colnames(traits)[-1]) %>%
rownames_to_column("Trait.comb") %>%
mutate_at(.vars=vars(Short_english_name:Long_English_name), ~ifelse(is.na(.), trait.name, {.}))
env <- read_delim("_data/Mesobromion/env.v2.10perc.cov.txt", delim="\t")
+species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
traits <- traits %>%
as.data.frame() %>%
@@ -81,6 +82,34 @@ trait.recap <- traits %>%
by=c("Variable" = "Short_english_name")) %>%
dplyr::select(Code=Variable, Variable=Long_English_name, everything())
+## calculate trait coverage across plots, both on p\a and cover
+trait.coverage <- species.cov %>%
+ pivot_longer(values_to="cover", names_to = "species", -RELEVE_NR) %>%
+ filter(cover >0) %>%
+ ## attach traits
+ left_join(traits %>%
+ rownames_to_column("species"),
+ by="species") %>%
+ mutate_at(.vars = vars(-RELEVE_NR, -species, -cover),
+ .funs = list(~if_else(is.na(.),0,1) * cover)) %>%
+ group_by(RELEVE_NR) %>%
+ summarize_at(.vars= vars(-species, -cover),
+ .funs = list(xxxcov=~sum(.),
+ xxxpa=~mean({.}>0))) %>%
+ dplyr::select(RELEVE_NR, order(colnames(.))) %>%
+ pivot_longer(-RELEVE_NR, values_to="coverage", names_to="trait" ) %>%
+ separate(trait, into = c("trait", "abundance"), sep="_xxx") %>%
+ group_by(trait, abundance) %>%
+ summarize(median.cov.perc=round(median(coverage)*100,1),
+ mean.cov.perc=round(mean(coverage)*100,1)) %>%
+ rename(median=median.cov.perc,
+ mean=mean.cov.perc) %>%
+ pivot_longer(cols=median:mean, names_to="metric") %>%
+ pivot_wider(names_from=metric:abundance, values_from=value)
+
+
+trait.recap <- trait.recap %>%
+ left_join(trait.coverage, by=c("Code"="trait"))
write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
@@ -165,7 +194,6 @@ traits.sign.alone <- corXY.ci %>%
### Plot combinations one by one
mydata.byone <- corXY.ci %>%
- filter(trait1 != c("V_VER_absent1", "V_VER_present1")) %>% ##delete traits that shouldn't be there
filter(run=="alone") %>%
filter(ntraits==1)%>%
arrange(Coef.obs) %>%
@@ -407,18 +435,49 @@ corXY.ci <- corXY.comb.ci %>%
#species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
+#traits <- traits.backup
+categorical.traits <- colnames(traits)[which(sapply(traits, "is.factor"))]
+traits$LeafPersistence <- factor(traits$LeafPersistence,
+ levels=c("immergrün","sommergrün","überwinternd_grün", "vorsommergrün"),
+ labels=c("eg", "sg", "wg", "se"))
+#traits$Pollination <- factor(traits$Pollination,
+# levels=c("NEKTAR_HONIG_INSEKTEN","POLLEN","WIND" ),
+# labels=c("insects", "pollen", "wind"))
+##spread categorical traits into multiple columns
+traits.dummy <- traits %>%
+ rownames_to_column("species0")
+while(any(sapply(traits.dummy, class)=="factor")){
+ id <- which(sapply(traits.dummy, class) == "factor")[1]
+ id.name <- names(id)
+ traits.dummy <- traits.dummy %>%
+ dplyr::select(-id.name) %>%
+ left_join(traits.dummy %>%
+ dplyr::select(species0, id) %>%
+ mutate(value = 1) %>%
+ filter(complete.cases(.)) %>%
+ pivot_wider(values_from = value,
+ names_from = id.name,
+ names_prefix = paste0(id.name, "_"),
+ values_fill = 0,
+ ),
+ by = "species0")
+ #print(id.name)
+}
+
+## match trait selection
+#e.g., traits.sign.alone <- c("SLA", "SeedMass", "CSR", "LeafPersistence")
+
CWM.wide <- species.cov %>%
pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "rel.cov") %>%
as.tbl() %>%
filter(rel.cov>0) %>%
arrange(RELEVE_NR) %>%
## attach traits
- left_join(traits %>%
- rownames_to_column("species0"),
+ left_join(traits.dummy,
by="species0") %>%
group_by(RELEVE_NR) %>%
- select_if(.predicate=~is.numeric(.)) %>% # CWM calculated ONLY for numeric traits
- summarize_at(.vars = vars(any_of(traits.sign.alone)), #GF_Macrophan:Rosette),
+ select_if(.predicate=~is.numeric(.)) %>%
+ summarize_at(.vars = vars(any_of(starts_with(as.character(traits.sign.alone)))),
.funs = list(~weighted.mean(., rel.cov, na.rm=T))) %>%
dplyr::select(RELEVE_NR, order(colnames(.))) %>%
column_to_rownames("RELEVE_NR")
@@ -443,6 +502,7 @@ circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
}
dat <- circleFun(diameter = 2, npoints = 100)
### 4.0 PCA of X (Fuzzy-weighted) matrix + CWMs ####
+#best.4traits <- c("Height", "BL_Dur", "VP_Fragm", "CSR")
source("01b_MesobromionCluster.R")
W.fuzzy <- matrix.x(comm=species.cov %>%
column_to_rownames("RELEVE_NR"),
@@ -468,7 +528,10 @@ myvectors <- as.data.frame(env.cor) %>%
mutate(category="Trait")) %>%
mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
- mutate(mycol=ifelse(category=="Trait", oilgreen, orange))
+ mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
+ mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
+ rowwise() %>%
+ mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab))
basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
theme_bw() +
@@ -477,27 +540,36 @@ basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
scale_x_continuous(limits=c(-1, 1)) + coord_equal() +
theme(panel.grid = element_blank())
-PCA.fuzz1_2 <- basemap0 +
+(PCA.fuzz1_2 <- basemap0 +
geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:3]*5),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors %>%
+ filter(categorical==0),
aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,# %>%
#mutate_if(~is.numeric(.), ~(.)*1.2),
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
+)
PCA.fuzz2_3 <- basemap0 +
geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:3]*5),
aes(x=PC2, y=PC3), pch="+", size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors%>%
+ filter(categorical==0),
aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC2, y=PC3, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,# %>%
#mutate_if(~is.numeric(.), ~(.)*1.2),
aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2,
@@ -520,10 +592,14 @@ tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>%
PCAfuzz1_2.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors %>%
+ filter(categorical==0),
aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
@@ -532,10 +608,14 @@ PCAfuzz1_2.sp <- basemap0 %+% tmp +
PCAfuzz2_3.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC2, y=PC3, label=labels), size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors %>%
+ filter(categorical==0),
aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC2, y=PC3, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7), segment.size = 0.5) +
@@ -550,8 +630,6 @@ ggsave("_pics/FigSXXXb_PCA_Fuzzy_2-3_wSpecies.png", width=8, height=8, dpi=300,
-
-
#### 4.1 PCA of Y (Bealls) matrix + CWM ####
W.beals <- as.data.frame(beals(species.cov %>%
column_to_rownames("RELEVE_NR") %>%
@@ -600,7 +678,10 @@ myvectors <- as.data.frame(env.cor) %>%
mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
- mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol))
+ mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol)) %>%
+ mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
+ rowwise() %>%
+ mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab))
@@ -615,10 +696,14 @@ basemap0 <- ggplot(data=as.data.frame(pca.out$CA$u[,1:4]*5)) +
PCA1_2 <- basemap0 +
geom_point(data=as.data.frame(pca.out$CA$u[,1:2]*5),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors %>%
+ filter(categorical==0),
aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
@@ -629,10 +714,14 @@ PCA1_2 <- basemap0 +
PCA3_4 <- basemap0 +
geom_point(data=as.data.frame(pca.out$CA$u[,3:4]*5),
aes(x=PC3, y=PC4), pch="+", size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors%>%
+ filter(categorical==0),
aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC3, y=PC4, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC3, y=PC4, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7), segment.size = 0.5) +
@@ -656,10 +745,14 @@ tmp <- as.data.frame(pca.out$CA$v[,1:4]*5) %>%
PCA1_2.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors%>%
+ filter(categorical==0),
aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
@@ -668,10 +761,14 @@ PCA1_2.sp <- basemap0 %+% tmp +
PCA3_4.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC3, y=PC4, label=labels), size=2, alpha=0.8) +
- geom_segment(data=myvectors,
+ geom_segment(data=myvectors%>%
+ filter(categorical==0),
aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
+ geom_point(data=myvectors %>%
+ filter(categorical==1),
+ aes(x=PC3, y=PC4, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
aes(x=PC3, y=PC4, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7), segment.size = 0.5) +
@@ -713,14 +810,11 @@ ggsave("_pics/Fig6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PC
# assign weights to traits [Growth forms and root structures are coded as binary, but should be downweighted]
#myweights <- c(rep(1/9, 9), rep(1/17, 17), rep(1, 26))
#traits.gowdis <- FD::gowdis(traits, w=myweights, ord="podani")
-data2 <- traits %>%
+data2 <- traits.dummy %>%
dplyr::select_if(.predicate = ~is.numeric(.)) %>% ### CAUTION -- ONLY NUMERIC TRAITS SHOWN
- dplyr::select(any_of(traits.sign.alone)) %>%
- rownames_to_column("species") %>%
- filter(species %in% colnames(species.cov)) %>%
- column_to_rownames("species")
-
-corr1 <- cor(data2, use="pairwise.complete.obs")
+ dplyr::select(any_of(starts_with(as.character(traits.sign.alone))))
+
+corr1 <- cor(data2, use="complete.obs")
pca2 <- princomp(covmat=corr1)
#plot(pca2$loadings[,2]~pca2$loadings[,1])
#text(pca2$loadings[,1], pca2$loadings[,2], dimnames(pca2$loadings)[[1]])
@@ -734,6 +828,13 @@ colnames(pca.scores)<-c('pca1','pca2','pca3','pca4','pca5') #just quickly naming
dat <- circleFun(diameter = 2, npoints = 100)
+myvector3 <- as.data.frame(pca2$loadings[,1:5]) %>%
+ rownames_to_column("mylab") %>%
+ mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
+ rowwise() %>%
+ mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab)) %>%
+ mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic"))
+
baseplot <- ggplot(data=as.data.frame(pca.scores*.2)) +
scale_color_identity() +
coord_equal(xlim = c(-1, 1), ylim=c(-1, 1)) +
@@ -744,28 +845,31 @@ baseplot <- ggplot(data=as.data.frame(pca.scores*.2)) +
PCA.t1 <- baseplot +
geom_point(aes(x=pca1, y=pca2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
- aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col=oilgreen),
+ geom_segment(data=myvector3 %>%
+ filter(categorical==0),
+ aes(x=0, xend=Comp.1, y=0, yend=Comp.2), col=oilgreen,
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
- geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
- aes(x=Comp.1, y=Comp.2, label=Trait, col=oilgreen, fontface=fontface0), size=2,
+ geom_point(data=myvector3 %>%
+ filter(categorical==1),
+ aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ geom_label_repel(data=myvector3,
+ aes(x=Comp.1, y=Comp.2, label=mylab, fontface=fontface0), col=oilgreen, size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
#ggtitle("PCoA of species-level traits")
PCA.t2 <- baseplot +
geom_point(aes(x=pca3, y=pca4), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
+ geom_segment(data=myvector3 %>%
+ filter(categorical==0),
aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col=oilgreen),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
- geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
- #geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
- aes(x=Comp.3, y=Comp.4, label=Trait, col=oilgreen, fontface=fontface0), size=2,
+ geom_point(data=myvector3 %>%
+ filter(categorical==1),
+ aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ geom_label_repel(data=myvector3,
+ aes(x=Comp.3, y=Comp.4, label=mylab, fontface=fontface0), col=oilgreen, size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
@@ -776,7 +880,7 @@ ggsave("_pics/S6_PCA_Traits_1-4_only7.png", width=10, height=5, dpi=300, PC_trai
##### 4.3b Alternative version of figS6, showing the species ####
tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
- mutate(species0=rownames(pca.scores)) %>%
+ mutate(species0=traits$species0) %>%
mutate(species=species0) %>%
separate(species0, sep="_", into=c("Gen", "Spe")) %>%
mutate(Gen=substr(Gen, 1, 3)) %>%
@@ -786,28 +890,31 @@ tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
PCA.t1.sp <- baseplot %+% tmp +
geom_text(aes(x=pca1, y=pca2, label=labels), size=2, alpha=0.7) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
+ geom_segment(data=myvector3 %>%
+ filter(categorical==0),
aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col=oilgreen),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
- geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
- aes(x=Comp.1, y=Comp.2, label=Trait, col=oilgreen, fontface=fontface0), size=2,
+ geom_point(data=myvector3 %>%
+ filter(categorical==1),
+ aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ geom_label_repel(data=myvector3,
+ aes(x=Comp.1, y=Comp.2, label=mylab, fontface=fontface0), col=oilgreen, size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
PCA.t2.sp <- baseplot %+% tmp +
geom_text(aes(x=pca3, y=pca4, label=labels), size=2, alpha=0.7) +
- geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
+ geom_segment(data=myvector3 %>%
+ filter(categorical==0),
aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col=oilgreen),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
- geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
- #geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
- rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
- aes(x=Comp.3, y=Comp.4, label=Trait, col=oilgreen, fontface=fontface0), size=2,
- position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
+ geom_label_repel(data=myvector3,
+ aes(x=Comp.3, y=Comp.4, label=mylab, fontface=fontface0),col=oilgreen, size=2,
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
+ geom_point(data=myvector3 %>%
+ filter(categorical==1),
+ aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
@@ -829,7 +936,7 @@ library(sp)
library(sf)
coords <- read_delim(file = "_data/Mesobromion/Mes2_somecoordinates.csv", delim = ";")
-env <- read_delim(file = "_data/Mesobromion/env.10perc.txt", delim = "\t")
+env <- read_delim(file = "_data/Mesobromion/env.v2.10perc.cov.txt", delim = "\t")
env0 <- read_delim("_data/Mesobromion/GVRD_MES2_site.csv", delim = ",")
#header <- "/data/sPlot/users/Francesco/Project_11/Germany/_data/tvhabita.dbf"
@@ -867,7 +974,7 @@ env.all.sf <- SpatialPointsDataFrame(coords=env.all %>% dplyr::select(LON, LAT),
theme_bw()
)
-ggsave(filename="_pics/PlotDistribution.png", height
+ggsave(filename="_pics/S2_PlotDistribution.png", height
=8, width=6, dpi=300, basemap2)
@@ -876,13 +983,17 @@ ggsave(filename="_pics/PlotDistribution.png", height
#### 6 Correlation matrices ####
#### 6.1 Trait level #####
library(corrplot)
-traits7 <- traits.cov %>%
- rownames_to_column("species") %>%
- filter(species %in% colnames(species.cov)) %>%
- column_to_rownames("species") %>%
+traits7 <- traits %>%
dplyr::select(any_of(traits.sign.alone)) %>%
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(all_of(best.2traits), everything())
+ relocate(all_of(best.4traits), everything()) %>%
+ select_if(.predicate=~is.numeric(.)) ## CAUTION only numerical variables shown
+
+## alternative using also dummy variables
+#traits7 <- traits.dummy %>%
+# dplyr::select(any_of(starts_with(as.character(traits.sign.alone)))) %>%
+# dplyr::select(sort(tidyselect::peek_vars())) %>%
+# relocate(all_of(starts_with(as.character(best.4traits))), everything())
res <- cor(traits7, use = "pairwise.complete.obs")
png(file="_pics/S7_Correlations_Trait.png", width=8, height=6.5, units = "in", res=300)
@@ -893,8 +1004,9 @@ dev.off()
#### 6.2 CWM level ####
res2 <- cor(CWM.wide %>%
+ dplyr::select(any_of(traits.sign.alone)) %>% ## caution selecting only numerical variables
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(all_of(best.2traits), everything()))
+ relocate(any_of(best.4traits), everything()))
png(file="_pics/S8_Correlations_CWMs.png", width=8, height=6.5, units = "in", res=300)
corrplot(res2, type = "upper",
tl.col = "black", tl.srt = 45, number.cex=0.6, addCoef.col = "black", diag=F)