Pollution inhibits wild plant regeneration from seeds by one-third
2025-03-17
Authors and emails: Qing-Qing Xu (xuqingqing22@mails.ucas.ac.cn); Si-Chong Chen (chensichong0528@gmail.com)
Brief introduction:
We quantify the effects of pollution on four aspects of plant
regeneration, including germination percentage (final seed germination
percentage), germination speed (mean germination time and germination
speed index), seedling growth (weight or length of root, shoot or whole
seedling), and biomass allocation (root-shoot ratio).
Explanation for abbreviation and variables in code:
- Abbreviations in the code:
GP: Germination percentage.
GS: Germination speed.
SG: Seedling growth.
BA: Biomass allocation.
m_control, sd_control, n_control: mean, sd and sample
size of control.
m_treat, sd_treat, n_treat: mean, sd and sample size of
pollution treatment.
Ref_ID, Effectsize_ID: identity of papers and effect
sizes.
- Predict variables:
Plant characteristics:
- Adult lifespan: short-lived (annual, biennial) and long-lived (perennial).
- Growth form: herbaceous and woody.
- Dormancy: MD, MPD, PD, PY, and ND.
- Seed mass: mean values of seed mass, unit is “mg”.
- Seedling part: above-, below- and total growth for seedling growth.
Pollution:
- Pollutant type: heavy metals, inorganic pollutants, organic pollutants, and emerging pollutants.
- Pollutant concentration: concentration of heavy metals, unit is “mg/L”.
1 Prepare data
1.1 Load packages
library(tidyverse)
library(metafor)
library(knitr)
library(kableExtra)1.2 Load data
library(readxl)
meta <- read_excel("p_data.xlsx", sheet = "data")
meta <- meta %>%
select(Effectsize_ID, Ref_ID, Species,
m_control, sd_control, n_control,
m_treat, sd_treat, n_treat,
Aspect, AspectType,
Pollution, Pollutant, Concentration, Unit,
SeedMass, Dormancy, GrowthForm, LifeSpan,
Lon, Lat)
glimpse(meta)## Rows: 4,476
## Columns: 21
## $ Effectsize_ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
## $ Ref_ID <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Species <chr> "Aporosa cardiosperma", "Aporosa cardiosperma", "Aporosa…
## $ m_control <dbl> 1.468610, 1.468610, 1.468610, 75.313810, 75.313810, 75.3…
## $ sd_control <chr> "7.8479999999999994E-2", "7.8479999999999994E-2", "7.847…
## $ n_control <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ m_treat <dbl> 1.3789200, 0.6614300, 0.4820600, 79.8326400, 63.2636000,…
## $ sd_treat <chr> "0.25785000000000002", "0.1009", "0.19058", "4.142260000…
## $ n_treat <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ Aspect <chr> "ger_index", "ger_index", "ger_index", "ger_per", "ger_p…
## $ AspectType <chr> "GS", "GS", "GS", "GP", "GP", "GP", "BE", "BE", "BE", "A…
## $ Pollution <chr> "Inorganic", "Inorganic", "Inorganic", "Inorganic", "Ino…
## $ Pollutant <chr> "gamma radiation", "gamma radiation", "gamma radiation",…
## $ Concentration <chr> "25", "50", "100", "25", "50", "100", "25", "50", "100",…
## $ Unit <chr> "Gy", "Gy", "Gy", "Gy", "Gy", "Gy", "Gy", "Gy", "Gy", "G…
## $ SeedMass <dbl> 46.0000, 46.0000, 46.0000, 46.0000, 46.0000, 46.0000, 46…
## $ Dormancy <chr> "ND", "ND", "ND", "ND", "ND", "ND", "ND", "ND", "ND", "N…
## $ GrowthForm <chr> "Woody", "Woody", "Woody", "Woody", "Woody", "Woody", "W…
## $ LifeSpan <chr> "Long-lived", "Long-lived", "Long-lived", "Long-lived", …
## $ Lon <chr> "74.809730555555547", "74.809730555555547", "74.80973055…
## $ Lat <chr> "13.497105555555555", "13.497105555555555", "13.49710555…1.3 Fill missing SD
# Select the corresponding column
meta$sd_control <- as.numeric(meta$sd_control)
meta$sd_treat <- as.numeric(meta$sd_treat)
# Summary of missingness:
library(metagear)
metadat_impute <- meta %>% dplyr::select(sd_treat, sd_control) %>% impute_missingness()##
## Summary of missingness:
##
## COLUMN PERCENT_MISSINGNESS IMPUTATIONS
## sd_treat 12 534
## sd_control 11 494
##
## Total missingness: 11% (1028 imputations needed)# Fill SD
meta <- impute_SD(as.data.frame(meta),
columnSDnames = c("sd_treat", "sd_control"),
columnXnames = c("m_treat", "m_control"),
method = "Bracken1992", range = 3, M = 1)1.4 Calculate yi and vi
meta <- escalc(data = meta,
measure = "ROM",
m1i = m_treat, sd1i = sd_treat, n1i = n_treat,
m2i = m_control, sd2i = sd_control, n2i = n_control)
# Convert sign for MGT
meta$yi <- ifelse(meta$Aspect == "MGT", -meta$yi, meta$yi)
# Log10-transformed seed mass
meta$Logmass <- log10(meta$SeedMass)
# Log10-transformed concentration
meta$Concentration <- as.numeric(meta$Concentration)
meta$Logcon <- log10(meta$Concentration)
# Convert log response ratio to percentage (%)
eff_per <- function(x) (exp(x) - 1) * 100
meta$per_yi <- eff_per(meta$yi)
# Save final data
write_csv(meta, "meta.csv")# Load final data
metadata <- read_csv("meta.csv")2 Data summary
2.1 Distributions of species and pollutants
dat_map <- metadata %>%
select(Species, Effectsize_ID, Ref_ID, Pollution, AspectType, GrowthForm, LifeSpan, Lon, Lat)
# Classify for Aspect
dat_map$AspectType[dat_map$AspectType %in% c("AB", "TB", "BE")] <- "SG"
# Summary for herbaceous and woody
dat_sum_growth <- dat_map %>%
group_by(GrowthForm) %>%
summarise(Species_count = n_distinct(Species),
Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
knitr::kable(dat_sum_growth, digits = 0,
caption = "Summary of growth form",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| GrowthForm | Species_count | Effectsize_count | Reference_count |
|---|---|---|---|
| Herbaceous | 78 | 2451 | 54 |
| Woody | 84 | 2025 | 58 |
# Summary for short- and long-lived
dat_sum_life <- dat_map %>%
group_by(LifeSpan) %>%
summarise(Species_count = n_distinct(Species),
Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
knitr::kable(dat_sum_life, digits = 0,
caption = "Summary of adult lifespan",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| LifeSpan | Species_count | Effectsize_count | Reference_count |
|---|---|---|---|
| Long-lived | 139 | 3461 | 92 |
| Short-lived | 23 | 1015 | 23 |
World map of pollution and species growth form
# Removal of duplicate Species for map
dat_map$Lon <- as.numeric(dat_map$Lon)
dat_map$Lat <- as.numeric(dat_map$Lat)
dat_fig_map <- distinct(dat_map, Species, Lon, Lat, Pollution, GrowthForm) %>% drop_na(Lon)
# World map of pollution and species growth form
library(ggthemes)
library(RColorBrewer)
world <- map_data("world")
worldmap <- ggplot() +
geom_map(data = world, map = world,
aes(long, lat, map_id = region),
fill = "#C4C4C4", alpha = 0.5, col = NA) +
geom_point(data = dat_fig_map,
aes(x = Lon, y = Lat, color = Pollution, shape = GrowthForm),
size = 3.5, alpha = 0.8, stroke = 2) +
scale_shape_manual(values = c(Herbaceous = 2, Woody = 1),
name = "Growth form") +
scale_color_manual(name = "Pollutant type", values = brewer.pal(4, "Set2")) +
coord_fixed(ratio = 0.5) +
coord_quickmap() +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank(),
panel.background = element_rect(fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
plot.margin = unit(c(0.05, 0.05, 0.05, 0.05), "cm")) +
theme_map() +
theme(legend.position = c(0.04, 0.11),
legend.key.height = unit(1.1, "cm"),
legend.key.width = unit(0.6, "cm"),
legend.text = element_text(family = "serif", size = 28),
legend.title = element_text(family = "serif", size = 28))
# Save figure
ggsave("figure/pollmap.jpg", width = 15, height = 8, units = "in", dpi = 300)
worldmap
2.2 Counts for regeneration aspects
# Summary of pollution
dat_bar <- metadata %>% select(AspectType, Pollution, Effectsize_ID)
# Classify for Aspect
dat_bar$AspectType[dat_bar$AspectType %in% c("AB", "TB", "BE")] <- "SG"
# Counts of pollution types
dat_bar <- dat_bar %>%
group_by(Pollution, AspectType) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
.groups = "drop") %>%
group_by(Pollution) %>%
mutate(Percentage = (Effectsize_count / sum(Effectsize_count)) * 100)
knitr::kable(dat_bar, digits = 2,
caption = "Summary of pollution types",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Pollution | AspectType | Effectsize_count | Percentage |
|---|---|---|---|
| Emerging | BA | 19 | 9.69 |
| Emerging | GP | 103 | 52.55 |
| Emerging | SG | 74 | 37.76 |
| Inorganic | BA | 92 | 17.90 |
| Inorganic | GP | 93 | 18.09 |
| Inorganic | GS | 77 | 14.98 |
| Inorganic | SG | 252 | 49.03 |
| Metal | BA | 591 | 17.47 |
| Metal | GP | 939 | 27.76 |
| Metal | GS | 200 | 5.91 |
| Metal | SG | 1653 | 48.86 |
| Organic | BA | 57 | 14.88 |
| Organic | GP | 137 | 35.77 |
| Organic | GS | 29 | 7.57 |
| Organic | SG | 160 | 41.78 |
Figure for sample sizes across pollutant types
dat_bar$Pollution <- factor(dat_bar$Pollution, levels = rev(c("Metal", "Inorganic", "Organic", "Emerging")))
dat_bar$AspectType <- factor(dat_bar$AspectType, levels = rev(c("GP", "GS", "SG", "BA")))
aspect_bar <- ggplot(dat_bar, aes(x = Pollution, y = Percentage / 100, fill = AspectType)) +
geom_bar(stat = "identity", position = "fill",
width = 0.6, colour = "black", linewidth = 1.2) +
geom_text(aes(x = Pollution, label = Effectsize_count),
position = position_stack(vjust = 0.5),
hjust = 0.5, size = 7.7, color="black", family = "serif") +
scale_fill_manual(
name = "Plant regeneration aspects",
breaks = c("GP", "GS", "SG", "BA"),
values = c("GP" = "#9F79EE","GS" = "#00B2EE",
"SG" = "#43CD80", "BA" = "#FF6347"),
labels = c("Germination percentage", "Germination speed",
"Seedling growth", "Biomass allocation")) +
scale_y_continuous(limits = c(0, 1.05),
breaks = seq(0, 1, 0.25),
labels = c("0", "25%", "50%", "75%", "100%"),
name = "Proportion of effect sizes",
expand = c(0, 0)) +
scale_x_discrete(expand = c(0.15, 0)) +
theme(axis.line = element_line(colour = "black", linewidth = 1.2),
axis.ticks = element_line(colour = "black", linewidth = 1.2),
axis.ticks.length = unit(0.2, "cm"),
axis.text = element_text(size = 30, colour = "black", family = "serif"),
axis.title.y = element_blank(),
axis.title.x = element_text(size = 30, family = "serif"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
legend.position = "right",
legend.key.height = unit(1.7, "cm"),
legend.key.width = unit(0.9, "cm"),
legend.title = element_text(size = 30, family = "serif"),
legend.text = element_text(size = 30, family = "serif"),
legend.background = element_rect(fill = "transparent"),
legend.box.background = element_rect(fill = "transparent", color = NA)) +
guides(fill = guide_legend(nrow = 4, byrow = T, title.position = "top")) +
coord_flip()
# Save figure
ggsave("figure/aspect.bar.jpg", width = 15, height = 5, units = "in", dpi = 300)
aspect_bar
2.3 Species counts
# library(devtools)
# devtools::install_github("jinyizju/V.PhyloMaker2")
library(V.PhyloMaker2)
library(plantlist)
# Data
dat_sp <- metadata %>%
select(Species, GrowthForm, AspectType, Logmass, Dormancy, LifeSpan) %>%
distinct(Species, .keep_all = T) %>% as.data.frame()
# TPL for species classification
sp.list <- subset(TPL(dat_sp$Species), select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
# One species (Arthrocaulon macrostachyum) could not retrieve the family
# We added it manually
sp.list[is.na(sp.list)] <- "Amaranthaceae"
colnames(sp.list) <- c("Species", "Genus", "Family")
# Family count
family.count <- sp.list %>%
group_by(Family) %>%
summarise(count = n()) %>%
filter(count >= 5) %>%
arrange(desc(count))
knitr::kable(family.count, digits = 0,
caption = "Summary of families with at least 5 species",
font_size = 20, align = "c", full_width = T, "html") %>%
kable_styling(position = "center")| Family | count |
|---|---|
| Fabaceae | 29 |
| Poaceae | 14 |
| Amaranthaceae | 12 |
| Asteraceae | 11 |
| Brassicaceae | 7 |
| Pinaceae | 6 |
| Cistaceae | 5 |
| Polygonaceae | 5 |
3 Meta-analyses
3.1 Prepare regeneration data
Subset data
# Aspects database
dat_gp <- filter(metadata, Aspect %in% "ger_per") %>% as.data.frame()
dat_gs <- filter(metadata, Aspect %in% c("ger_index", "MGT")) %>% as.data.frame()
dat_sg <- filter(metadata, Aspect %in% c("root_length", "root_dry", "root_fresh",
"shoot_length", "shoot_dry", "shoot_fresh",
"seedling_length", "seedling_dry", "seedling_fresh")) %>% as.data.frame()
dat_ba <- filter(metadata, Aspect %in% "rs_ratio") %>% as.data.frame()Heavy metal concentration for plant regeneration
# Chose heavy metals for subset data
dat_hm <- filter(metadata, Pollution == "Metal" & Unit == "mg/L")
# Summary for heavy metals
dat_sum_hm <- dat_hm %>%
group_by(Pollutant) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
arrange(desc(Effectsize_count))
knitr::kable(dat_sum_hm, digits = 0,
caption = "Summary of pollutant of heavy metals",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Pollutant | Effectsize_count | Reference_count |
|---|---|---|
| Zn | 576 | 19 |
| Cu | 558 | 18 |
| Pb | 520 | 26 |
| Cd | 388 | 24 |
| Ni | 284 | 11 |
| Mn | 163 | 4 |
| As | 139 | 5 |
| Fe | 98 | 2 |
| Hg | 92 | 5 |
| Al | 62 | 2 |
| Cr | 39 | 3 |
| Sr | 13 | 1 |
# Classify heavy metals
dat_hm$Pollutant[!(dat_hm$Pollutant %in% c("As", "Cd", "Cr", "Hg", "Pb"))] <- "others"
# Heavy metals for different aspect
dat_hm_gp <- filter(dat_hm, Aspect %in% "ger_per") %>% as.data.frame()
dat_hm_gs <- filter(dat_hm, Aspect %in% c("ger_index", "MGT")) %>% as.data.frame()
dat_hm_sg <- filter(dat_hm, Aspect %in% c("shoot_length", "shoot_dry", "shoot_fresh",
"root_length", "root_dry", "root_fresh",
"seedling_length", "seedling_dry", "seedling_fresh")) %>% as.data.frame()
dat_hm_ba <- filter(dat_hm, Aspect %in% "rs_ratio") %>% as.data.frame()Calculate VCV matrix
# In a experiments, the use of the same control is noted as the common ID
# Prepare for common ID
numbering <- function(col) {
num <- numeric(length(col))
cur_num <- 1
prev_val <- col[1]
for (i in 1:length(col)) {
if (col[i] != prev_val) {
cur_num <- cur_num + 1
prev_val <- col[i]
}
num[i] <- cur_num
}
return(num)
}
# Common ID for database
dat_gs$common_ID <- numbering(dat_gs[, "m_control"])
dat_gp$common_ID <- numbering(dat_gp[, "m_control"])
dat_sg$common_ID <- numbering(dat_sg[, "m_control"])
dat_ba$common_ID <- numbering(dat_ba[, "m_control"])
# Common ID for heavy metals database
dat_hm_gs$common_ID <- numbering(dat_hm_gs[, "m_control"])
dat_hm_gp$common_ID <- numbering(dat_hm_gp[, "m_control"])
dat_hm_sg$common_ID <- numbering(dat_hm_sg[, "m_control"])
dat_hm_ba$common_ID <- numbering(dat_hm_ba[, "m_control"])
# Formulation for VCV matrix
calc.v <- function(x) {
v <- matrix((x$sd_control[1]^2 / (x$n_control[1] * x$m_control[1]^2)),
nrow = nrow(x), ncol = nrow(x))
diag(v) <- x$vi
v
}
# Calculate VCV matrix
v.gs <- bldiag(lapply(split(dat_gs, dat_gs$common_ID), calc.v))
v.gp <- bldiag(lapply(split(dat_gp, dat_gp$common_ID), calc.v))
v.sg <- bldiag(lapply(split(dat_sg, dat_sg$common_ID), calc.v))
v.ba <- bldiag(lapply(split(dat_ba, dat_ba$common_ID), calc.v))
v.hm.gs <- bldiag(lapply(split(dat_hm_gs, dat_hm_gs$common_ID), calc.v))
v.hm.gp <- bldiag(lapply(split(dat_hm_gp, dat_hm_gp$common_ID), calc.v))
v.hm.sg <- bldiag(lapply(split(dat_hm_sg, dat_hm_sg$common_ID), calc.v))
v.hm.ba <- bldiag(lapply(split(dat_hm_ba, dat_hm_ba$common_ID), calc.v))3.2 Prepare seedling parts data
Seedling different parts, containing above-ground, below-ground and seedling total biomass
# Data
dat_ab <- filter(metadata, Aspect %in% c("shoot_length", "shoot_dry", "shoot_fresh")) %>% as.data.frame()
dat_be <- filter(metadata, Aspect %in% c("root_length", "root_dry", "root_fresh")) %>% as.data.frame()
dat_sl <- filter(metadata, Aspect %in% c("seedling_length", "seedling_dry", "seedling_fresh")) %>% as.data.frame()
# Heavy metal data
dat_hm_ab <- filter(dat_hm, Aspect %in% c("shoot_length", "shoot_dry", "shoot_fresh")) %>% as.data.frame()
dat_hm_be <- filter(dat_hm, Aspect %in% c("root_length", "root_dry", "root_fresh")) %>% as.data.frame()
dat_hm_sl <- filter(dat_hm, Aspect %in% c("seedling_length", "seedling_dry", "seedling_fresh")) %>% as.data.frame()
# Calculate VCV matrix
dat_ab$common_ID <- numbering(dat_ab[, "m_control"])
dat_be$common_ID <- numbering(dat_be[, "m_control"])
dat_sl$common_ID <- numbering(dat_sl[, "m_control"])
v.ab <- bldiag(lapply(split(dat_ab, dat_ab$common_ID), calc.v))
v.be <- bldiag(lapply(split(dat_be, dat_be$common_ID), calc.v))
v.sl <- bldiag(lapply(split(dat_sl, dat_sl$common_ID), calc.v))
dat_hm_ab$common_ID <- numbering(dat_hm_ab[, "m_control"])
dat_hm_be$common_ID <- numbering(dat_hm_be[, "m_control"])
dat_hm_sl$common_ID <- numbering(dat_hm_sl[, "m_control"])
v.hm.ab <- bldiag(lapply(split(dat_hm_ab, dat_hm_ab$common_ID), calc.v))
v.hm.be <- bldiag(lapply(split(dat_hm_be, dat_hm_be$common_ID), calc.v))
v.hm.sl <- bldiag(lapply(split(dat_hm_sl, dat_hm_sl$common_ID), calc.v))3.3 Overall effects of pollution on regeneration
Calculate mean effect size
# Quantify the effects of pollution on each plant regeneration Aspect separately
gp.mean.es <- rma.mv(yi, v.gp, mods = ~1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.mean.es <- rma.mv(yi, v.gs, mods = ~1, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.mean.es <- rma.mv(yi, v.sg, mods = ~1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.mean.es <- rma.mv(yi, v.ba, mods = ~1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary of results
mean.es.result.table <- data.frame(AspectType = c("GP", "GS", "SG", "BA"),
Mean = c(gp.mean.es$b, gs.mean.es$b,
sg.mean.es$b, ba.mean.es$b),
SE = c(gp.mean.es$se, gs.mean.es$se,
sg.mean.es$se, ba.mean.es$se),
Percentage = c(eff_per(gp.mean.es$b), eff_per(gs.mean.es$b),
eff_per(sg.mean.es$b), eff_per(ba.mean.es$b)),
Lower_95_CI = c(gp.mean.es$ci.lb, gs.mean.es$ci.lb,
sg.mean.es$ci.lb, ba.mean.es$ci.lb),
Upper_95_CI = c(gp.mean.es$ci.ub, gs.mean.es$ci.ub,
sg.mean.es$ci.ub, ba.mean.es$ci.ub),
Lower_95_CI_change = c(eff_per(gp.mean.es$ci.lb),
eff_per(gs.mean.es$ci.lb),
eff_per(sg.mean.es$ci.lb),
eff_per(ba.mean.es$ci.lb)),
Upper_95_CI_change = c(eff_per(gp.mean.es$ci.ub),
eff_per(gs.mean.es$ci.ub),
eff_per(sg.mean.es$ci.ub),
eff_per(ba.mean.es$ci.ub)),
P = c(gp.mean.es$pval, gs.mean.es$pval,
sg.mean.es$pval, ba.mean.es$pval),
Qt = c(gp.mean.es$QE, gs.mean.es$QE,
sg.mean.es$QE, ba.mean.es$QE),
Qt_p = c(gp.mean.es$QEp, gs.mean.es$QEp,
sg.mean.es$QEp, ba.mean.es$QEp),
N = c(gp.mean.es$k, gs.mean.es$k,
sg.mean.es$k, ba.mean.es$k))
knitr::kable(mean.es.result.table, digits = 3,
caption = "Summary of meta-analytical results",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| AspectType | Mean | SE | Percentage | Lower_95_CI | Upper_95_CI | Lower_95_CI_change | Upper_95_CI_change | P | Qt | Qt_p | N |
|---|---|---|---|---|---|---|---|---|---|---|---|
| GP | -0.404 | 0.053 | -33.222 | -0.507 | -0.300 | -39.780 | -25.949 | 0.000 | 4244188.1 | 0 | 1272 |
| GS | -0.343 | 0.115 | -29.009 | -0.568 | -0.118 | -43.308 | -11.102 | 0.003 | 693312.3 | 0 | 306 |
| SG | -0.539 | 0.076 | -41.666 | -0.689 | -0.389 | -49.785 | -32.235 | 0.000 | 1080965.9 | 0 | 2139 |
| BA | -0.315 | 0.111 | -27.044 | -0.533 | -0.097 | -41.339 | -9.267 | 0.005 | 10966066.3 | 0 | 759 |
# Save data
write_csv(mean.es.result.table, "table/mean.es.result.table.csv")3.3.1 Figure for mean effect size
Theme for figure
# Theme for figure
theme_plot <- theme(
axis.ticks = element_line(colour = "black", linewidth = 1.5),
axis.ticks.length = unit(0.2, "cm"),
axis.title = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
legend.position = "none",
plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
panel.background = element_blank(),
panel.border = element_rect(linewidth = 1.5, fill = NA))
# Order for pollutant types
Predictortype.1_level <- c("Metal", "Inorganic", "Organic", "Emerging")
# Order for plant characteristics
Predictortype.2_level <- c("ND", "PY", "PD", "MPD", "MD", "Short-lived", "Long-lived", "Woody", "Herbaceous")
# Curly brace for significance
library(pBrackets)
bracketsGrob <- function(...){
l <- list(...)
e <- new.env()
e$l <- l
grid:::recordGrob( {
do.call(grid.brackets, l)
}, e)
}Main figure
library(gghalves)
# Data
mean.effsize <- metadata %>%
select(Effectsize_ID, AspectType, yi, vi, per_yi)
mean.effsize$AspectType[mean.effsize$AspectType %in% c("AB", "TB", "BE")] <- "SG"
mean.effsize$AspectType <- factor(mean.effsize$AspectType, levels = rev(c("GP", "GS", "SG", "BA")))
effsize.fig <- ggplot(data = mean.effsize,
mapping = aes(x = AspectType, y = yi, fill = AspectType)) +
geom_half_violin(aes(fill = AspectType, colour = AspectType),
position = position_nudge(x = 0),
linewidth = 0.5, alpha = 0.5, trim = TRUE, side = "r") +
geom_jitter(aes(stage(start = AspectType, after_scale = x - 0.3), color = AspectType),
alpha = 0.3, shape = 1, stroke = 1.1, width = 0.2) +
geom_pointrange(data = mean.es.result.table,
mapping = aes(y = Mean, ymin = Lower_95_CI, ymax = Upper_95_CI),
position = position_nudge(x = 0), linewidth = 1.5, size = 0.7) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("GP" = "#9F79EE",
"GS" = "#00B2EE",
"SG" = "#43CD80",
"BA" = "#FF8C69")) +
scale_color_manual(values = c("GP" = "#9F79EE",
"GS" = "#00B2EE",
"SG" = "#43CD80",
"BA" = "#FF8C69")) +
scale_x_discrete(limits = rev(c("GP", "GS", "SG", "BA")),
breaks = c("GP", "GS", "SG", "BA"),
labels = c("Germination\n percentage", "Germination\n speed",
"Seedling\n growth", "Biomass\n allocation")) +
scale_y_continuous(name = "Effect size (Log response ratio)") +
theme(axis.title.x.bottom = element_text(size = 30, color = "black",
family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 30, colour = "black", family = "serif"),
axis.text.y = element_text(size = 30, color = "black", family = "serif")) +
theme_plot +
guides(fill = "none") +
coord_flip()
# Save figure
ggsave("figure/effsize.size.fig.jpg", width = 9, height = 8, units = "in", dpi = 300)
effsize.fig
Mean effect sizes shown in percentage
# Convert to percentage
mean.es.result.table$SE_per <- eff_per(mean.es.result.table$SE)
mean.es.result.table$L95_percent <- mean.es.result.table$Percentage - 1.96 *mean.es.result.table$SE_per
mean.es.result.table$U95_percent <- mean.es.result.table$Percentage + 1.96 *mean.es.result.table$SE_per
effsize.per.fig <- ggplot(data = mean.es.result.table,
mapping = aes(x = AspectType)) +
geom_pointrange(mapping = aes(y = Percentage, ymin = L95_percent, ymax = U95_percent),
position = position_nudge(x = 0), linewidth = 1.2, size = 0.9) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.2) +
scale_x_discrete(limits = rev(c("GP", "GS", "SG", "BA")),
breaks = c("GP", "GS", "SG", "BA"),
labels = c("Germination\n percentage", "Germination\n speed",
"Seedling\n growth", "Biomass\n allocation")) +
scale_y_continuous(name = "Change (%)") +
theme(axis.title.x.bottom = element_text(size = 30, color = "black", family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 30, colour = "black", family = "serif"),
axis.text.y = element_text(size = 30, color = "black", family = "serif"),
legend.position = "none") +
theme_plot +
guides(fill = "none") +
coord_flip()
effsize.per.fig
Mean effect sizes shown in percentage
alleffsize.per.fig <- ggplot(data = mean.effsize,
mapping = aes(x = AspectType, y = per_yi, fill = AspectType)) +
geom_jitter(aes(stage(start = AspectType, after_scale = x), color = AspectType),
alpha = 0.3, shape = 1, stroke = 1.1, width = 0.4) +
geom_pointrange(data = mean.es.result.table,
mapping = aes(y = Percentage, ymin = L95_percent, ymax = U95_percent),
position = position_nudge(x = 0), linewidth = 1.2, size = 0.9) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.2) +
scale_color_manual(values = c("GP" = "#9F79EE",
"GS" = "#00B2EE",
"SG" = "#43CD80",
"BA" = "#FF8C69")) +
scale_x_discrete(limits = rev(c("GP", "GS", "SG", "BA")),
breaks = c("GP", "GS", "SG", "BA"),
labels = c("Germination\n percentage", "Germination\n speed",
"Seedling\n growth", "Biomass\n allocation")) +
scale_y_continuous(limits = c(-100, 1300), breaks = c(-100, 0, 300, 600, 900, 1200),
name = "Change (%)") +
theme(axis.title.x.bottom = element_text(size = 30, color = "black",
family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 30, colour = "black", family = "serif"),
axis.text.y = element_text(size = 30, color = "black", family = "serif")) +
theme_plot +
guides(fill = "none") +
coord_flip()
alleffsize.per.fig
3.4 Phylogeny signal of Blomberg’s K
3.4.1 Different species effect sizes
Germination percentage
gp.result.sp <- rma.mv(yi, v.gp, mods = ~Species - 1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
dat_sum_gp_sp <- dat_gp %>%
pivot_longer(cols = c(Species),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
gp.sp.table <- data.frame(Predictor = "Species",
Predictortype = c(sort(unique(dat_gp$Species))),
Mean = c(gp.result.sp$b),
SE = c(gp.result.sp$se),
Lower_95_CI = c(gp.result.sp$ci.lb),
Upper_95_CI = c(gp.result.sp$ci.ub),
P = c(gp.result.sp$pval))
# Combine table
gp_sp_result <- full_join(dat_sum_gp_sp, gp.sp.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything()) %>% as.data.frame()
# Save data
write_csv(gp_sp_result, "table/gp_sp_result.csv")Germination speed
gs.result.sp <- rma.mv(yi, v.gs, mods = ~Species - 1, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
dat_sum_gs_sp <- dat_gs %>%
pivot_longer(cols = c(Species),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
gs.sp.table <- data.frame(Predictor = "Species",
Predictortype = c(sort(unique(dat_gs$Species))),
Mean = c(gs.result.sp$b),
SE = c(gs.result.sp$se),
Lower_95_CI = c(gs.result.sp$ci.lb),
Upper_95_CI = c(gs.result.sp$ci.ub),
P = c(gs.result.sp$pval))
# Combine table
gs_sp_result <- full_join(dat_sum_gs_sp, gs.sp.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything()) %>% as.data.frame()
# Save data
write_csv(gs_sp_result, "table/gs_sp_result.csv")Seedling growth
sg.result.sp <- rma.mv(yi, v.sg, mods = ~Species - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
dat_sum_sg_sp <- dat_sg %>%
pivot_longer(cols = c(Species),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
sg.sp.table <- data.frame(Predictor = "Species",
Predictortype = c(sort(unique(dat_sg$Species))),
Mean = c(sg.result.sp$b),
SE = c(sg.result.sp$se),
Lower_95_CI = c(sg.result.sp$ci.lb),
Upper_95_CI = c(sg.result.sp$ci.ub),
P = c(sg.result.sp$pval))
# Combine table
sg_sp_result <- full_join(dat_sum_sg_sp, sg.sp.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything()) %>% as.data.frame()
# Save data
write_csv(sg_sp_result, "table/sg_sp_result.csv")Biomass allocation
ba.result.sp <- rma.mv(yi, v.ba, mods = ~Species - 1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
dat_sum_ba_sp <- dat_ba %>%
pivot_longer(cols = c(Species),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
ba.sp.table <- data.frame(Predictor = "Species",
Predictortype = c(sort(unique(dat_ba$Species))),
Mean = c(ba.result.sp$b),
SE = c(ba.result.sp$se),
Lower_95_CI = c(ba.result.sp$ci.lb),
Upper_95_CI = c(ba.result.sp$ci.ub),
P = c(ba.result.sp$pval))
# Combine table
ba_sp_result <- full_join(dat_sum_ba_sp, ba.sp.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything()) %>% as.data.frame()
# Save data
write_csv(ba_sp_result, "table/ba_sp_result.csv")3.4.2 Test for phylogenetic signals
Germination percentage
library(phytools)
# Data
dat_gp_sp <- gp.sp.table %>% select(Predictortype) %>% as.data.frame()
# TPL for species classification
sp.gp.list <- subset(TPL(dat_gp_sp$Predictortype), select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
colnames(sp.gp.list) <- c("Species", "Genus", "Family")
gp.species.tree <- phylo.maker(sp.gp.list, scenario = "S3")
gp.species.tree <- gp.species.tree$scenario.3
gp.sp.mean <- gp.sp.table %>% select(Predictortype, Mean) %>%
mutate(Predictortype = gsub(" ", "_", Predictortype)) %>%
column_to_rownames(var = "Predictortype")
# Test phylogenetic signals
gp.sp.mean <- setNames(gp.sp.mean$Mean, rownames(gp.sp.mean))
phylosig(gp.species.tree, gp.sp.mean, test = TRUE)##
## Phylogenetic signal K : 0.034756
## P-value (based on 1000 randomizations) : 0.754Germination speed
# Data
dat_gs_sp <- gs.sp.table %>% select(Predictortype) %>% as.data.frame()
# TPL for species classification
sp.gs.list <- subset(TPL(dat_gs_sp$Predictortype), select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
colnames(sp.gs.list) <- c("Species", "Genus", "Family")
gs.species.tree <- phylo.maker(sp.gs.list, scenario = "S3")
gs.species.tree <- gs.species.tree$scenario.3
gs.sp.mean <- gs.sp.table %>% select(Predictortype, Mean) %>%
mutate(Predictortype = gsub(" ", "_", Predictortype)) %>%
column_to_rownames(var = "Predictortype")
# Test phylogenetic signals
gs.sp.mean <- setNames(gs.sp.mean$Mean, rownames(gs.sp.mean))
phylosig(gs.species.tree, gs.sp.mean, test = TRUE)##
## Phylogenetic signal K : 0.0397494
## P-value (based on 1000 randomizations) : 0.962Seedling growth
# Data
dat_sg_sp <- sg.sp.table %>% select(Predictortype) %>% as.data.frame()
# TPL for species classification
sp.sg.list <- subset(TPL(dat_sg_sp$Predictortype), select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
sp.sg.list[is.na(sp.sg.list)] <- "Amaranthaceae"
colnames(sp.sg.list) <- c("Species", "Genus", "Family")
sg.species.tree <- phylo.maker(sp.sg.list, scenario = "S3")
sg.species.tree <- sg.species.tree$scenario.3
sg.sp.mean <- sg.sp.table %>% select(Predictortype, Mean) %>%
mutate(Predictortype = gsub(" ", "_", Predictortype)) %>%
column_to_rownames(var = "Predictortype")
# Test phylogenetic signals
sg.sp.mean <- setNames(sg.sp.mean$Mean, rownames(sg.sp.mean))
phylosig(sg.species.tree, sg.sp.mean, test = TRUE)##
## Phylogenetic signal K : 0.0546722
## P-value (based on 1000 randomizations) : 0.545Biomass allocation
# Data
dat_ba_sp <- ba.sp.table %>% select(Predictortype) %>% as.data.frame()
# TPL for species classification
sp.ba.list <- subset(TPL(dat_ba_sp$Predictortype), select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
sp.ba.list[is.na(sp.ba.list)] <- "Amaranthaceae"
colnames(sp.ba.list) <- c("Species", "Genus", "Family")
ba.species.tree <- phylo.maker(sp.ba.list, scenario = "S3")
ba.species.tree <- ba.species.tree$scenario.3
ba.sp.mean <- ba.sp.table %>% select(Predictortype, Mean) %>%
mutate(Predictortype = gsub(" ", "_", Predictortype)) %>%
column_to_rownames(var = "Predictortype")
# Test phylogenetic signals
ba.sp.mean <- setNames(ba.sp.mean$Mean, rownames(ba.sp.mean))
phylosig(ba.species.tree, ba.sp.mean, test = TRUE)##
## Phylogenetic signal K : 0.0610259
## P-value (based on 1000 randomizations) : 0.5733.4.3 Figure for phylogeny
library(ggtree)
library(ggnewscale)
library(viridis)
# Phylogeny tree
sp.tree <- phylo.maker(sp.list, scenario = "S3")
sp.tree <- sp.tree$scenario.3
## Read tree
sp.tree <- fortify(sp.tree)
row.names(dat_sp) <- sub(" ", "_", dat_sp$Species)
row.names(gp_sp_result) <- sub(" ", "_", gp_sp_result$Predictortype)
row.names(gs_sp_result) <- sub(" ", "_", gs_sp_result$Predictortype)
row.names(sg_sp_result) <- sub(" ", "_", sg_sp_result$Predictortype)
row.names(ba_sp_result) <- sub(" ", "_", ba_sp_result$Predictortype)
## Figure
line.tree <- ggtree(sp.tree,
layout = "rectangular",
size = 0.6) +
xlim(0, max(sp.tree$branch.length) * 5.5)
line.gp <- gheatmap(line.tree, gp_sp_result[, "Mean", drop = FALSE], colnames = FALSE, width = 0.2, offset = -25) +
scale_fill_gradientn(name = "Log response ratio of\nplant regeneration",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = guide_colorbar(order = 1, ncol = 1)) +
new_scale_fill()
line.gs <-gheatmap(line.gp, gs_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 41) +
scale_fill_gradientn(name = "Germination speed",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.sg <-gheatmap(line.gs, sg_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 107) +
scale_fill_gradientn(name = "Seedling growth",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.ba <-gheatmap(line.sg, ba_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 173) +
scale_fill_gradientn(name = "Biomass allocation",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.g <- gheatmap(line.ba, dat_sp[, "GrowthForm", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 244) +
scale_fill_manual(name = "Growth form",
breaks = c("Herbaceous", "Woody"),
values = c("#228B22", "#8B4500"),
guide = guide_legend(order = 1, ncol = 1)) +
new_scale_fill()
line.l <- gheatmap(line.g, dat_sp[, "LifeSpan", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 310) +
scale_fill_manual(name = "Adult lifespan",
breaks = c("Short-lived", "Long-lived"),
values = c("#00CDCD", "#FF8C00"),
guide = guide_legend(order = 2, ncol = 1)) +
new_scale_fill()
lin.d <- gheatmap(line.l, dat_sp[, "Dormancy", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 376) +
scale_fill_manual(name = "Dormancy class",
breaks = c("ND", "PD", "PY", "MD", "MPD"),
values = scales::hue_pal()(5),
guide = guide_legend(order = 3, ncol = 1)) +
new_scale_fill()
line.s <- gheatmap(lin.d, dat_sp[, "Logmass", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 442) +
scale_fill_viridis_c(name = "Seed mass (mg)",
option = "viridis", direction = -1,
breaks = c(log10(0.1), log10(1), log10(10), 2, 3, 4),
labels = c("0.1", "1", "10", "100", "1000", "10000"),
guide = guide_colorbar(order = 4, ncol = 1)) +
new_scale_fill() +
theme(
legend.position = c(0.12, 0.6),
legend.box = "vertical",
legend.key.height = unit(0.6, "cm"),
legend.key.width = unit(0.35, "cm"),
legend.text = element_text(size = 18, family = "serif"),
legend.title = element_text(size = 18, family = "serif")
)
lintree.family <- line.s +
geom_strip("Mimosa_scabrella", "Hymenaea_courbaril", label = "Fabaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 45, family ="serif") + # 29 sp
geom_strip("Lolium_perenne", "Spartina_densiflora", label = "Poaceae",
barsize = 1, fontsize = 7,hjust = 0.5, offset = 540,
offset.text = 40, family = "serif")+ # 14 sp
geom_strip("Arthrocaulon_macrostachyum", "Salicornia_ramosissima", label = "Amaranthaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 72, family = "serif")+ # 12 sp
geom_strip("Taraxacum_campylodes", "Tagetes_minuta", label = "Asteraceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 53, family = "serif") + # 11 sp
geom_strip("Thlaspi_arvense", "Coronopus_didymus", label = "Brassicaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 62, family = "serif") + # 7 sp
geom_strip("Picea_mariana", "Pinus_banksiana", label = "Pinaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 43, family = "serif") + # 6 sp
geom_strip("Halimium_atriplicifolium", "Cistus_salviifolius", label = "Cistaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 47, family = "serif") + # 5 sp
geom_strip("Rheum_palmatum", "Rumex_maritimus", label = "Polygonaceae",
barsize = 1, fontsize = 7, hjust = 0.5, offset = 540,
offset.text = 66, family = "serif") # 5 sp
# Save figure
ggsave("figure/line.tree.jpg", width = 15, height = 10, units = "in", dpi = 500)
lintree.family
Tree with species names
line.tree <- ggtree(sp.tree,
layout = "rectangular",
size = 1) +
xlim(0, max(sp.tree$branch.length) * 5.5) +
geom_tiplab(aes(label = label), hjust = -0.05, size = 8, family ="serif")
line.gp <- gheatmap(line.tree, gp_sp_result[, "Mean", drop = FALSE], colnames = FALSE, width = 0.2, offset = 95) +
scale_fill_gradientn(name = "Log response ratio of\nplant regeneration",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = guide_colorbar(order = 1, ncol = 1)) +
new_scale_fill()
line.gs <-gheatmap(line.gp, gs_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 161) +
scale_fill_gradientn(name = "Germination speed",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.sg <-gheatmap(line.gs, sg_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 227) +
scale_fill_gradientn(name = "Seedling growth",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.ba <-gheatmap(line.sg, ba_sp_result[, "Mean", drop = FALSE],
colnames = FALSE, width = 0.2, offset = 293) +
scale_fill_gradientn(name = "Biomass allocation",
colors = viridis(100, option = "H", begin = 0, end = 1),
na.value = "white",
limits = c(-4.1, 2.6), breaks = c(-4, -2, 0, 2),
labels = c("-4", "-2", "0", "2"),
guide = "none") +
new_scale_fill()
line.g <- gheatmap(line.ba, dat_sp[, "GrowthForm", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 364) +
scale_fill_manual(name = "Growth form",
breaks = c("Herbaceous", "Woody"),
values = c("#228B22", "#8B4500"),
guide = guide_legend(order = 1, ncol = 1)) +
new_scale_fill()
line.l <- gheatmap(line.g, dat_sp[, "LifeSpan", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 430) +
scale_fill_manual(name = "Adult lifespan",
breaks = c("Short-lived", "Long-lived"),
values = c("#00CDCD", "#FF8C00"),
guide = guide_legend(order = 2, ncol = 1)) +
new_scale_fill()
lin.d <- gheatmap(line.l, dat_sp[, "Dormancy", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 496) +
scale_fill_manual(name = "Dormancy class",
breaks = c("ND", "PD", "PY", "MD", "MPD"),
values = scales::hue_pal()(5),
guide = guide_legend(order = 3, ncol = 1)) +
new_scale_fill()
line.s <- gheatmap(lin.d, dat_sp[, "Logmass", drop = FALSE], colnames = FALSE, width = 0.2,
offset = 562) +
scale_fill_viridis_c(name = "Seed mass (mg)",
option = "viridis", direction = -1,
breaks = c(log10(0.1), log10(1), log10(10), 2, 3, 4),
labels = c("0.1", "1", "10", "100", "1000", "10000"),
guide = guide_colorbar(order = 4, ncol = 1)) +
new_scale_fill() +
theme(
legend.position = c(0.12, 0.59),
legend.box = "vertical",
legend.spacing.x = unit(0.7, "cm"),
legend.key.height = unit(3.5, "cm"),
legend.key.width = unit(1.8, "cm"),
legend.text = element_text(size = 50, family = "serif"),
legend.title = element_text(size = 50, family = "serif")
)
lintree.family <- line.s +
geom_strip("Mimosa_scabrella", "Hymenaea_courbaril", label = "Fabaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 44, family ="serif") + # 29 sp
geom_strip("Lolium_perenne", "Spartina_densiflora", label = "Poaceae",
barsize = 2, fontsize = 21,hjust = 0.5, offset = 660,
offset.text = 41, family = "serif")+ # 14 sp
geom_strip("Arthrocaulon_macrostachyum", "Salicornia_ramosissima", label = "Amaranthaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 72, family = "serif")+ # 12 sp
geom_strip("Taraxacum_campylodes", "Tagetes_minuta", label = "Asteraceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 51, family = "serif") + # 11 sp
geom_strip("Thlaspi_arvense", "Coronopus_didymus", label = "Brassicaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 60, family = "serif") + # 7 sp
geom_strip("Picea_mariana", "Pinus_banksiana", label = "Pinaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 43, family = "serif") + # 6 sp
geom_strip("Halimium_atriplicifolium", "Cistus_salviifolius", label = "Cistaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 46, family = "serif") + # 5 sp
geom_strip("Rheum_palmatum", "Rumex_maritimus", label = "Polygonaceae",
barsize = 2, fontsize = 21, hjust = 0.5, offset = 660,
offset.text = 65, family = "serif") # 5 sp
# Save figure
ggsave("figure/linetree.txt.jpg", width = 45, height = 40, units = "in", dpi = 500)
lintree.family
3.5 Introduced predictors
3.5.1 Q-statistic for predictors
3.5.1.1 Germination percentage
gp.qm.result.growth <- rma.mv(yi, v.gp, mods = ~GrowthForm, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.qm.result.life <- rma.mv(yi, v.gp, mods = ~LifeSpan, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.qm.result.d <- rma.mv(yi, v.gp, mods = ~Dormancy, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.qm.result.mass <- rma.mv(yi, v.gp, mods = ~1 + Logmass, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.qm.result.pt <- rma.mv(yi, v.gp, mods = ~Pollution, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Heavy metals
gp.qm.result.con <- rma.mv(yi, v.hm.gp, mods = ~1 + Logcon, data = dat_hm_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Result
gp.qm.result.table <- data.frame(Aspect = "GP",
Predictor = c("Pollution",
"Dormancy", "Growth form", "Life span",
"Seed mass", "Concentration"),
QM = c(gp.qm.result.pt$QM, gp.qm.result.d$QM,
gp.qm.result.growth$QM, gp.qm.result.life$QM,
gp.qm.result.mass$QM, gp.qm.result.con$QM),
QM_p = c(gp.qm.result.pt$QMp, gp.qm.result.d$QMp,
gp.qm.result.growth$QMp, gp.qm.result.life$QMp,
gp.qm.result.mass$QMp, gp.qm.result.con$QMp),
N = c(gp.qm.result.pt$k, gp.qm.result.d$k,
gp.qm.result.growth$k, gp.qm.result.life$k,
gp.qm.result.mass$k, gp.qm.result.con$k))
knitr::kable(gp.qm.result.table, digits = 3,
caption = "Summary of germination percentage predict variable QM",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Aspect | Predictor | QM | QM_p | N |
|---|---|---|---|---|
| GP | Pollution | 2.462 | 0.482 | 1272 |
| GP | Dormancy | 1.730 | 0.785 | 1272 |
| GP | Growth form | 3.677 | 0.055 | 1272 |
| GP | Life span | 1.363 | 0.243 | 1272 |
| GP | Seed mass | 3.661 | 0.056 | 1272 |
| GP | Concentration | 79.758 | 0.000 | 865 |
3.5.1.2 Germination speed
gs.qm.result.growth <- rma.mv(yi, v.gs, mods = ~GrowthForm, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.qm.result.d <- rma.mv(yi, v.gs, mods = ~Dormancy, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.qm.result.mass <- rma.mv(yi, v.gs, mods = ~1 + Logmass, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.qm.result.pt <- rma.mv(yi, v.gs, mods = ~Pollution, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Heavy metals
gs.qm.result.con <- rma.mv(yi, v.hm.gs, mods = ~1 + Logcon, data = dat_hm_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Result
gs.qm.result.table <- data.frame(Aspect = "GS",
Predictor = c("Pollution",
"Dormancy", "Growth form",
"Seed mass", "Concentration"),
QM = c(gs.qm.result.pt$QM, gs.qm.result.d$QM,
gs.qm.result.growth$QM, gs.qm.result.mass$QM,
gs.qm.result.con$QM),
QM_p = c(gs.qm.result.pt$QMp, gs.qm.result.d$QMp,
gs.qm.result.growth$QMp, gs.qm.result.mass$QMp,
gs.qm.result.con$QMp),
N = c(gs.qm.result.pt$k, gs.qm.result.d$k,
gs.qm.result.growth$k, gs.qm.result.mass$k,
gs.qm.result.con$k))
knitr::kable(gs.qm.result.table, digits = 3,
caption = "Summary of germination speed predict variable QM",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Aspect | Predictor | QM | QM_p | N |
|---|---|---|---|---|
| GS | Pollution | 0.634 | 0.728 | 306 |
| GS | Dormancy | 1.662 | 0.436 | 306 |
| GS | Growth form | 0.009 | 0.925 | 306 |
| GS | Seed mass | 0.292 | 0.589 | 306 |
| GS | Concentration | 31.506 | 0.000 | 181 |
3.5.1.3 Seedling growth
sg.qm.result.growth <- rma.mv(yi, v.sg, mods = ~GrowthForm, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.qm.result.life <- rma.mv(yi, v.sg, mods = ~LifeSpan, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.qm.result.d <- rma.mv(yi, v.sg, mods = ~Dormancy, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.qm.result.mass <- rma.mv(yi, v.sg, mods = ~1 + Logmass, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.qm.result.pt <- rma.mv(yi, v.sg, mods = ~Pollution, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Heavy metals
sg.qm.result.con <- rma.mv(yi, v.hm.sg, mods = ~1 + Logcon, data = dat_hm_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Result
sg.qm.result.table <- data.frame(Aspect = "SG",
Predictor = c("Pollution",
"Dormancy", "Growth form", "Life span",
"Seed mass", "Concentration"),
QM = c(sg.qm.result.pt$QM, sg.qm.result.d$QM,
sg.qm.result.growth$QM, sg.qm.result.life$QM,
sg.qm.result.mass$QM, sg.qm.result.con$QM),
QM_p = c(sg.qm.result.pt$QMp, sg.qm.result.d$QMp,
sg.qm.result.growth$QMp, sg.qm.result.life$QMp,
sg.qm.result.mass$QMp, sg.qm.result.con$QMp),
N = c(sg.qm.result.pt$k, sg.qm.result.d$k,
sg.qm.result.growth$k, sg.qm.result.life$k,
sg.qm.result.mass$k, sg.qm.result.con$k))
knitr::kable(sg.qm.result.table, digits = 3,
caption = "Summary of seedling growth predict variable QM",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Aspect | Predictor | QM | QM_p | N |
|---|---|---|---|---|
| SG | Pollution | 8.298 | 0.040 | 2139 |
| SG | Dormancy | 6.846 | 0.144 | 2139 |
| SG | Growth form | 2.392 | 0.122 | 2139 |
| SG | Life span | 0.836 | 0.360 | 2139 |
| SG | Seed mass | 0.234 | 0.628 | 2139 |
| SG | Concentration | 325.867 | 0.000 | 1394 |
3.5.1.4 Biomass allocation
ba.qm.result.growth <- rma.mv(yi, v.ba, mods = ~GrowthForm, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.qm.result.life <- rma.mv(yi, v.ba, mods = ~LifeSpan, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.qm.result.d <- rma.mv(yi, v.ba, mods = ~Dormancy, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.qm.result.mass <- rma.mv(yi, v.ba, mods = ~1 + Logmass, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.qm.result.pt <- rma.mv(yi, v.ba, mods = ~Pollution, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Heavy metals
ba.qm.result.con <- rma.mv(yi, v.hm.ba, mods = ~1 + Logcon, data = dat_hm_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Result
ba.qm.result.table <- data.frame(Aspect = "BA",
Predictor = c("Pollution",
"Dormancy", "Growth form", "Life span",
"Seed mass", "Concentration"),
QM = c(ba.qm.result.pt$QM, ba.qm.result.d$QM,
ba.qm.result.growth$QM, ba.qm.result.life$QM,
ba.qm.result.mass$QM, ba.qm.result.con$QM),
QM_p = c(ba.qm.result.pt$QMp, ba.qm.result.d$QMp,
ba.qm.result.growth$QMp, ba.qm.result.life$QMp,
ba.qm.result.mass$QMp, ba.qm.result.con$QMp),
N = c(ba.qm.result.pt$k, ba.qm.result.d$k,
ba.qm.result.growth$k, ba.qm.result.life$k,
ba.qm.result.mass$k, ba.qm.result.con$k))
knitr::kable(ba.qm.result.table, digits = 3,
caption = "Summary of biomass alloaction predict variable QM",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Aspect | Predictor | QM | QM_p | N |
|---|---|---|---|---|
| BA | Pollution | 5.466 | 0.141 | 759 |
| BA | Dormancy | 5.391 | 0.249 | 759 |
| BA | Growth form | 2.436 | 0.119 | 759 |
| BA | Life span | 0.308 | 0.579 | 759 |
| BA | Seed mass | 2.658 | 0.103 | 759 |
| BA | Concentration | 51.878 | 0.000 | 492 |
3.5.2 Comparisons for categorical variables
3.5.2.1 Germination percentage
gp.result.pt <- rma.mv(yi, v.gp, mods = ~Pollution - 1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.result.growth <- rma.mv(yi, v.gp, mods = ~GrowthForm - 1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.result.life <- rma.mv(yi, v.gp, mods = ~LifeSpan - 1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.result.d <- rma.mv(yi, v.gp, mods = ~Dormancy - 1, data = dat_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gp.result.hm.p <- rma.mv(yi, v.hm.gp, mods = ~Pollutant - 1, data = dat_hm_gp,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for germination percentage
dat_sum_gp <- dat_gp %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_gp %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Summary for heavy metal
dat_sum_hm_gp <- dat_hm_gp %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
# Germination percentage meta analysis result table
gp.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Heavy metal"),
c(4, 5, 2, 2, 6)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(gp.result.pt$b, gp.result.d$b,
gp.result.growth$b, gp.result.life$b,
gp.result.hm.p$b),
SE = c(gp.result.pt$se, gp.result.d$se,
gp.result.growth$se, gp.result.life$se,
gp.result.hm.p$se),
Lower_95_CI = c(gp.result.pt$ci.lb, gp.result.d$ci.lb,
gp.result.growth$ci.lb, gp.result.life$ci.lb,
gp.result.hm.p$ci.lb),
Upper_95_CI = c(gp.result.pt$ci.ub, gp.result.d$ci.ub,
gp.result.growth$ci.ub, gp.result.life$ci.ub,
gp.result.hm.p$ci.ub),
P = c(gp.result.pt$pval, gp.result.d$pval,
gp.result.growth$pval, gp.result.life$pval,
gp.result.hm.p$pval))
# Combine table
gp_result <- full_join(dat_sum_gp, gp.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(gp_result, digits = 3,
caption = "Summary of predict variables of germination percentage",
font_size = 20, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 21 | 1 | -0.134 | 0.456 | -1.028 | 0.759 | 0.768 |
| Dormancy | MPD | 21 | 1 | -0.439 | 0.456 | -1.332 | 0.455 | 0.336 |
| Dormancy | ND | 342 | 34 | -0.396 | 0.085 | -0.563 | -0.229 | 0.000 |
| Dormancy | PD | 652 | 42 | -0.458 | 0.076 | -0.607 | -0.309 | 0.000 |
| Dormancy | PY | 236 | 16 | -0.306 | 0.115 | -0.530 | -0.081 | 0.008 |
| Growth form | Herbaceous | 738 | 42 | -0.501 | 0.073 | -0.644 | -0.357 | 0.000 |
| Growth form | Woody | 534 | 45 | -0.305 | 0.074 | -0.451 | -0.160 | 0.000 |
| Life span | Long-lived | 999 | 71 | -0.379 | 0.057 | -0.490 | -0.268 | 0.000 |
| Life span | Short-lived | 273 | 17 | -0.528 | 0.119 | -0.761 | -0.295 | 0.000 |
| Pollution | Emerging | 103 | 2 | 0.143 | 0.379 | -0.599 | 0.885 | 0.706 |
| Pollution | Inorganic | 93 | 7 | -0.507 | 0.171 | -0.843 | -0.171 | 0.003 |
| Pollution | Metal | 939 | 58 | -0.409 | 0.063 | -0.532 | -0.286 | 0.000 |
| Pollution | Organic | 137 | 15 | -0.396 | 0.123 | -0.638 | -0.154 | 0.001 |
| Heavy metal | As | 49 | 3 | -0.741 | 0.223 | -1.179 | -0.303 | 0.001 |
| Heavy metal | Cd | 129 | 19 | -0.354 | 0.093 | -0.535 | -0.172 | 0.000 |
| Heavy metal | Cr | 10 | 3 | -0.862 | 0.231 | -1.316 | -0.409 | 0.000 |
| Heavy metal | Hg | 34 | 5 | -0.354 | 0.159 | -0.665 | -0.043 | 0.026 |
| Heavy metal | Pb | 139 | 21 | -0.328 | 0.089 | -0.503 | -0.153 | 0.000 |
| Heavy metal | others | 504 | 28 | -0.394 | 0.077 | -0.544 | -0.244 | 0.000 |
# Save data
write_csv(gp_result, "table/gp_result_predict variables.csv")Different level for comparisons
# Pollutant type
anova(gp.result.pt, L = c(0, 1, -1, 0)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: -0.0987 0.1820 -0.5422 0.5877
##
## Test of Hypothesis:
## QM(df = 1) = 0.2939, p-val = 0.5877anova(gp.result.pt, L = c(0, 1, 0, -1)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.1111 0.2086 -0.5327 0.5943
##
## Test of Hypothesis:
## QM(df = 1) = 0.2837, p-val = 0.5943anova(gp.result.pt, L = c(0, 0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.0124 0.1374 -0.0905 0.9279
##
## Test of Hypothesis:
## QM(df = 1) = 0.0082, p-val = 0.9279# Growth form
anova(gp.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.1953 0.1018 -1.9176 0.0552
##
## Test of Hypothesis:
## QM(df = 1) = 3.6771, p-val = 0.0552# Adult lifespan
anova(gp.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: 0.1492 0.1278 1.1674 0.2431
##
## Test of Hypothesis:
## QM(df = 1) = 1.3627, p-val = 0.2431# Dormancy
anova(gp.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: 0.0616 0.1094 0.5631 0.5734
##
## Test of Hypothesis:
## QM(df = 1) = 0.3171, p-val = 0.5734anova(gp.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: -0.0907 0.1398 -0.6486 0.5166
##
## Test of Hypothesis:
## QM(df = 1) = 0.4207, p-val = 0.5166anova(gp.result.d, L = c(0, 0, 0, 1, -1)) ##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: -0.1523 0.1351 -1.1267 0.2599
##
## Test of Hypothesis:
## QM(df = 1) = 1.2695, p-val = 0.25993.5.2.2 Germination speed
gs.result.pt <- rma.mv(yi, v.gs, mods = ~Pollution - 1, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.result.growth <- rma.mv(yi, v.gs, mods = ~GrowthForm - 1, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.result.d <- rma.mv(yi, v.gs, mods = ~Dormancy - 1, data = dat_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
gs.result.hm.p <- rma.mv(yi, v.hm.gs, mods = ~Pollutant - 1, data = dat_hm_gs,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for germination speed
dat_sum_gs <- dat_gs %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_gs %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Germination speed meta analysis result table
gs.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form", "Heavy metal"),
c(3, 3, 2, 6)),
Predictortype = c("Inorganic", "Metal", "Organic",
"ND", "PD", "PY",
"Herbaceous", "Woody",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(gs.result.pt$b, gs.result.d$b,
gs.result.growth$b, gs.result.hm.p$b),
SE = c(gs.result.pt$se, gs.result.d$se,
gs.result.growth$se, gs.result.hm.p$se),
Lower_95_CI = c(gs.result.pt$ci.lb, gs.result.d$ci.lb,
gs.result.growth$ci.lb, gs.result.hm.p$ci.lb),
Upper_95_CI = c(gs.result.pt$ci.ub, gs.result.d$ci.ub,
gs.result.growth$ci.ub, gs.result.hm.p$ci.ub),
P = c(gs.result.pt$pval, gs.result.d$pval,
gs.result.growth$pval, gs.result.hm.p$pval))
# Combine table
gs_result <- full_join(dat_sum_gs, gs.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(gs_result, digits = 3,
caption = "Summary of predict variables of germination speed for differnet levels",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | ND | 103 | 11 | -0.452 | 0.168 | -0.781 | -0.123 | 0.007 |
| Dormancy | PD | 154 | 13 | -0.385 | 0.168 | -0.715 | -0.055 | 0.022 |
| Dormancy | PY | 49 | 5 | -0.073 | 0.245 | -0.553 | 0.407 | 0.766 |
| Growth form | Herbaceous | 75 | 6 | -0.366 | 0.268 | -0.892 | 0.160 | 0.173 |
| Growth form | Woody | 231 | 21 | -0.338 | 0.130 | -0.592 | -0.084 | 0.009 |
| Life span | Long-lived | 294 | 26 | NA | NA | NA | NA | NA |
| Life span | Short-lived | 12 | 1 | NA | NA | NA | NA | NA |
| Pollution | Inorganic | 77 | 3 | -0.290 | 0.209 | -0.700 | 0.120 | 0.165 |
| Pollution | Metal | 200 | 19 | -0.407 | 0.144 | -0.688 | -0.125 | 0.005 |
| Pollution | Organic | 29 | 4 | -0.172 | 0.303 | -0.766 | 0.422 | 0.570 |
| Heavy metal | As | 6 | 1 | -0.674 | 0.175 | -1.017 | -0.331 | 0.000 |
| Heavy metal | Cd | 37 | 7 | -0.373 | 0.107 | -0.583 | -0.163 | 0.001 |
| Heavy metal | Cr | 5 | 1 | -0.306 | 0.179 | -0.657 | 0.046 | 0.089 |
| Heavy metal | Hg | 9 | 2 | -0.294 | 0.164 | -0.615 | 0.027 | 0.072 |
| Heavy metal | Pb | 52 | 9 | -0.262 | 0.102 | -0.461 | -0.063 | 0.010 |
| Heavy metal | others | 72 | 8 | -0.468 | 0.099 | -0.663 | -0.273 | 0.000 |
# Save data
write_csv(gs_result, "table/gs_result_predict variables.csv")Germination speed comparisons
# Pollutant type
anova(gs.result.pt, L = c(1, -1, 0))##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.1169 0.2296 0.5092 0.6106
##
## Test of Hypothesis:
## QM(df = 1) = 0.2593, p-val = 0.6106anova(gs.result.pt, L = c(1, 0, -1))##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.1176 0.3682 -0.3193 0.7495
##
## Test of Hypothesis:
## QM(df = 1) = 0.1019, p-val = 0.7495anova(gs.result.pt, L = c(0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.2345 0.3354 -0.6990 0.4846
##
## Test of Hypothesis:
## QM(df = 1) = 0.4886, p-val = 0.4846# Growth form
anova(gs.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.0280 0.2973 -0.0942 0.9249
##
## Test of Hypothesis:
## QM(df = 1) = 0.0089, p-val = 0.9249# Dormancy
anova(gs.result.d, L = c(1, -1, 0))##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: -0.0675 0.2093 -0.3223 0.7472
##
## Test of Hypothesis:
## QM(df = 1) = 0.1039, p-val = 0.7472anova(gs.result.d, L = c(1, 0, -1))##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: -0.3793 0.2958 -1.2822 0.1998
##
## Test of Hypothesis:
## QM(df = 1) = 1.6441, p-val = 0.1998anova(gs.result.d, L = c(0, 1, -1))##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: -0.3118 0.2969 -1.0501 0.2937
##
## Test of Hypothesis:
## QM(df = 1) = 1.1027, p-val = 0.29373.5.2.3 Seedling growth
sg.result.pt <- rma.mv(yi, v.sg, mods = ~Pollution - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.result.growth <- rma.mv(yi, v.sg, mods = ~GrowthForm - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.result.life <- rma.mv(yi, v.sg, mods = ~LifeSpan - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.result.d <- rma.mv(yi, v.sg, mods = ~Dormancy - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sg.result.hm.p <- rma.mv(yi, v.hm.sg, mods = ~Pollutant - 1, data = dat_hm_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Add categorical variable for seedling parts
sg.result.part <- rma.mv(yi, v.sg, mods = ~AspectType - 1, data = dat_sg,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for seedling growth
dat_sum_sg <- dat_sg %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan, AspectType),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
Predictor == "AspectType" ~ "Parts",
TRUE ~ Predictor
)) %>%
mutate(Predictortype = case_when(
Predictortype == "AB" ~ "Above-ground biomass",
Predictortype == "BE" ~ "Below-ground biomass",
Predictortype == "TB" ~ "Total biomass",
TRUE ~ Predictortype)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_sg %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Seedling growth meta analysis result table
sg.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Parts", "Heavy metal"),
c(4, 5, 2, 2, 3, 6)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"Above-ground biomass", "Below-ground biomass", "Total biomass",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(sg.result.pt$b, sg.result.d$b,
sg.result.growth$b, sg.result.life$b,
sg.result.part$b, sg.result.hm.p$b),
SE = c(sg.result.pt$se, sg.result.d$se,
sg.result.growth$se, sg.result.life$se,
sg.result.part$se, sg.result.hm.p$se),
Lower_95_CI = c(sg.result.pt$ci.lb, sg.result.d$ci.lb,
sg.result.growth$ci.lb, sg.result.life$ci.lb,
sg.result.part$ci.lb, sg.result.hm.p$ci.lb),
Upper_95_CI = c(sg.result.pt$ci.ub, sg.result.d$ci.ub,
sg.result.growth$ci.ub, sg.result.life$ci.ub,
sg.result.part$ci.ub, sg.result.hm.p$ci.ub),
P = c(sg.result.pt$pval, sg.result.d$pval,
sg.result.growth$pval, sg.result.life$pval,
sg.result.part$pval, sg.result.hm.p$pval))
# Combine table
sg_result <- full_join(dat_sum_sg, sg.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(sg_result, digits = 3,
caption = "Summary of predict variables of seedling growth for differnet levels",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 63 | 1 | -0.765 | 0.358 | -1.467 | -0.064 | 0.032 |
| Dormancy | MPD | 63 | 1 | -0.794 | 0.359 | -1.497 | -0.091 | 0.027 |
| Dormancy | ND | 638 | 33 | -0.362 | 0.104 | -0.566 | -0.159 | 0.000 |
| Dormancy | PD | 844 | 33 | -0.592 | 0.097 | -0.782 | -0.402 | 0.000 |
| Dormancy | PY | 531 | 19 | -0.747 | 0.132 | -1.007 | -0.488 | 0.000 |
| Growth form | Herbaceous | 1196 | 40 | -0.642 | 0.101 | -0.840 | -0.444 | 0.000 |
| Growth form | Woody | 943 | 42 | -0.438 | 0.100 | -0.634 | -0.243 | 0.000 |
| Life span | Long-lived | 1612 | 67 | -0.560 | 0.080 | -0.717 | -0.404 | 0.000 |
| Life span | Short-lived | 527 | 18 | -0.444 | 0.129 | -0.697 | -0.190 | 0.001 |
| Parts | Above-ground biomass | 859 | 55 | -0.339 | 0.079 | -0.493 | -0.184 | 0.000 |
| Parts | Below-ground biomass | 967 | 65 | -0.781 | 0.078 | -0.934 | -0.628 | 0.000 |
| Parts | Total biomass | 313 | 22 | -0.301 | 0.092 | -0.481 | -0.122 | 0.001 |
| Pollution | Emerging | 74 | 7 | -0.007 | 0.246 | -0.489 | 0.475 | 0.978 |
| Pollution | Inorganic | 252 | 6 | -0.250 | 0.262 | -0.764 | 0.264 | 0.341 |
| Pollution | Metal | 1653 | 56 | -0.643 | 0.086 | -0.812 | -0.474 | 0.000 |
| Pollution | Organic | 160 | 11 | -0.418 | 0.201 | -0.812 | -0.023 | 0.038 |
| Heavy metal | As | 58 | 4 | -0.858 | 0.226 | -1.301 | -0.414 | 0.000 |
| Heavy metal | Cd | 176 | 20 | -0.585 | 0.098 | -0.777 | -0.393 | 0.000 |
| Heavy metal | Cr | 16 | 2 | -0.348 | 0.207 | -0.753 | 0.058 | 0.093 |
| Heavy metal | Hg | 40 | 3 | 0.056 | 0.147 | -0.232 | 0.343 | 0.704 |
| Heavy metal | Pb | 259 | 20 | -0.527 | 0.097 | -0.717 | -0.337 | 0.000 |
| Heavy metal | others | 845 | 22 | -0.573 | 0.091 | -0.751 | -0.396 | 0.000 |
# Save data
write_csv(sg_result, "table/sg_result_predict variables.csv")Seedling growth comparisons
# Pollutant type
anova(sg.result.pt, L = c(0, 1, -1, 0))##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.3933 0.2757 1.4264 0.1538
##
## Test of Hypothesis:
## QM(df = 1) = 2.0346, p-val = 0.1538anova(sg.result.pt, L = c(0, 1, 0, -1))##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.1678 0.3307 0.5075 0.6118
##
## Test of Hypothesis:
## QM(df = 1) = 0.2576, p-val = 0.6118anova(sg.result.pt, L = c(0, 0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.2255 0.2189 -1.0301 0.3030
##
## Test of Hypothesis:
## QM(df = 1) = 1.0610, p-val = 0.3030anova(sg.result.pt, L = c(1, -1, 0, 0))##
## Hypothesis:
## 1: PollutionEmerging - PollutionInorganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.2430 0.3596 0.6757 0.4992
##
## Test of Hypothesis:
## QM(df = 1) = 0.4566, p-val = 0.4992anova(sg.result.pt, L = c(1, 0, -1, 0))##
## Hypothesis:
## 1: PollutionEmerging - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.6363 0.2528 2.5174 0.0118
##
## Test of Hypothesis:
## QM(df = 1) = 6.3375, p-val = 0.0118anova(sg.result.pt, L = c(1, 0, 0, -1))##
## Hypothesis:
## 1: PollutionEmerging - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.4108 0.3179 1.2922 0.1963
##
## Test of Hypothesis:
## QM(df = 1) = 1.6697, p-val = 0.1963# Growth form
anova(sg.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.2039 0.1318 -1.5465 0.1220
##
## Test of Hypothesis:
## QM(df = 1) = 2.3916, p-val = 0.1220# Adult lifespan
anova(sg.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: -0.1166 0.1276 -0.9145 0.3605
##
## Test of Hypothesis:
## QM(df = 1) = 0.8363, p-val = 0.3605# Dormancy
anova(sg.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: 0.2296 0.1222 1.8780 0.0604
##
## Test of Hypothesis:
## QM(df = 1) = 3.5270, p-val = 0.0604anova(sg.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.3849 0.1605 2.3984 0.0165
##
## Test of Hypothesis:
## QM(df = 1) = 5.7522, p-val = 0.0165anova(sg.result.d, L = c(0, 0, 0, 1, -1))##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.1553 0.1466 1.0594 0.2894
##
## Test of Hypothesis:
## QM(df = 1) = 1.1223, p-val = 0.28943.5.2.4 Seedling parts
Above-ground biomass
ab.result.pt <- rma.mv(yi, v.ab, mods = ~Pollution - 1, data = dat_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ab.result.growth <- rma.mv(yi, v.ab, mods = ~GrowthForm - 1, data = dat_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ab.result.life <- rma.mv(yi, v.ab, mods = ~LifeSpan - 1, data = dat_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ab.result.d <- rma.mv(yi, v.ab, mods = ~Dormancy - 1, data = dat_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ab.result.hm.p <- rma.mv(yi, v.hm.ab, mods = ~Pollutant - 1, data = dat_hm_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for results
dat_sum_ab <- dat_ab %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_ab %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Summary for heavy metal
dat_sum_hm_ab <- dat_hm_ab %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop")
# Above-ground biomass meta analysis result table
ab.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Heavy metal"),
c(4, 5, 2, 2, 6)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(ab.result.pt$b, ab.result.d$b,
ab.result.growth$b, ab.result.life$b,
ab.result.hm.p$b),
SE = c(ab.result.pt$se, ab.result.d$se,
ab.result.growth$se, ab.result.life$se,
ab.result.hm.p$se),
Lower_95_CI = c(ab.result.pt$ci.lb, ab.result.d$ci.lb,
ab.result.growth$ci.lb, ab.result.life$ci.lb,
ab.result.hm.p$ci.lb),
Upper_95_CI = c(ab.result.pt$ci.ub, ab.result.d$ci.ub,
ab.result.growth$ci.ub, ab.result.life$ci.ub,
ab.result.hm.p$ci.ub),
P = c(ab.result.pt$pval, ab.result.d$pval,
ab.result.growth$pval, ab.result.life$pval,
ab.result.hm.p$pval))
# Combine table
ab_result <- full_join(dat_sum_ab, ab.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(ab_result, digits = 3,
caption = "Summary of predict variables of above-ground biomass",
font_size = 20, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 21 | 1 | -0.617 | 0.343 | -1.289 | 0.056 | 0.072 |
| Dormancy | MPD | 21 | 1 | -0.703 | 0.349 | -1.388 | -0.019 | 0.044 |
| Dormancy | ND | 255 | 24 | -0.324 | 0.110 | -0.541 | -0.108 | 0.003 |
| Dormancy | PD | 353 | 23 | -0.389 | 0.105 | -0.595 | -0.184 | 0.000 |
| Dormancy | PY | 209 | 14 | -0.608 | 0.139 | -0.880 | -0.335 | 0.000 |
| Growth form | Herbaceous | 478 | 28 | -0.504 | 0.109 | -0.719 | -0.290 | 0.000 |
| Growth form | Woody | 381 | 29 | -0.326 | 0.109 | -0.539 | -0.113 | 0.003 |
| Life span | Long-lived | 657 | 47 | -0.458 | 0.086 | -0.627 | -0.289 | 0.000 |
| Life span | Short-lived | 202 | 11 | -0.204 | 0.136 | -0.470 | 0.063 | 0.134 |
| Pollution | Emerging | 22 | 4 | -0.030 | 0.306 | -0.629 | 0.569 | 0.923 |
| Pollution | Inorganic | 116 | 6 | -0.190 | 0.245 | -0.670 | 0.291 | 0.439 |
| Pollution | Metal | 651 | 39 | -0.435 | 0.097 | -0.625 | -0.245 | 0.000 |
| Pollution | Organic | 70 | 7 | -0.670 | 0.232 | -1.125 | -0.215 | 0.004 |
| Heavy metal | As | 26 | 3 | -0.398 | 0.218 | -0.825 | 0.029 | 0.067 |
| Heavy metal | Cd | 47 | 12 | -0.469 | 0.101 | -0.667 | -0.272 | 0.000 |
| Heavy metal | Cr | 8 | 2 | -0.328 | 0.192 | -0.706 | 0.049 | 0.088 |
| Heavy metal | Hg | 9 | 1 | -0.211 | 0.178 | -0.559 | 0.138 | 0.236 |
| Heavy metal | Pb | 94 | 14 | -0.412 | 0.096 | -0.600 | -0.223 | 0.000 |
| Heavy metal | others | 337 | 15 | -0.243 | 0.079 | -0.397 | -0.088 | 0.002 |
# Save data
write_csv(ab_result, "table/ab_result_predict variables.csv")Above-ground biomass comparisons
# Pollution
anova(ab.result.pt, L = c(0, 1, -1, 0)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.2451 0.2635 0.9304 0.3522
##
## Test of Hypothesis:
## QM(df = 1) = 0.8656, p-val = 0.3522anova(ab.result.pt, L = c(0, 1, 0, -1)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.4799 0.3377 1.4211 0.1553
##
## Test of Hypothesis:
## QM(df = 1) = 2.0194, p-val = 0.1553anova(ab.result.pt, L = c(0, 0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.2348 0.2515 0.9334 0.3506
##
## Test of Hypothesis:
## QM(df = 1) = 0.8713, p-val = 0.3506# Growth form
anova(ab.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.1782 0.1416 -1.2586 0.2082
##
## Test of Hypothesis:
## QM(df = 1) = 1.5840, p-val = 0.2082# Life span
anova(ab.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: -0.2540 0.1301 -1.9525 0.0509
##
## Test of Hypothesis:
## QM(df = 1) = 3.8122, p-val = 0.0509# Dormancy
anova(ab.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: 0.0649 0.1268 0.5122 0.6085
##
## Test of Hypothesis:
## QM(df = 1) = 0.2623, p-val = 0.6085anova(ab.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.2833 0.1669 1.6978 0.0895
##
## Test of Hypothesis:
## QM(df = 1) = 2.8826, p-val = 0.0895anova(ab.result.d, L = c(0, 0, 0, 1, -1)) ##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.2184 0.1507 1.4488 0.1474
##
## Test of Hypothesis:
## QM(df = 1) = 2.0991, p-val = 0.1474Below-ground biomass
be.result.pt <- rma.mv(yi, v.be, mods = ~Pollution - 1, data = dat_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.result.growth <- rma.mv(yi, v.be, mods = ~GrowthForm - 1, data = dat_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.result.life <- rma.mv(yi, v.be, mods = ~LifeSpan - 1, data = dat_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.result.d <- rma.mv(yi, v.be, mods = ~Dormancy - 1, data = dat_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.result.hm.p <- rma.mv(yi, v.hm.be, mods = ~Pollutant - 1, data = dat_hm_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for results
dat_sum_be <- dat_be %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_be %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Below-ground biomass meta analysis result table
be.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Heavy metal"),
c(4, 5, 2, 2, 6)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(be.result.pt$b, be.result.d$b,
be.result.growth$b, be.result.life$b,
be.result.hm.p$b),
SE = c(be.result.pt$se, be.result.d$se,
be.result.growth$se, be.result.life$se,
be.result.hm.p$se),
Lower_95_CI = c(be.result.pt$ci.lb, be.result.d$ci.lb,
be.result.growth$ci.lb, be.result.life$ci.lb,
be.result.hm.p$ci.lb),
Upper_95_CI = c(be.result.pt$ci.ub, be.result.d$ci.ub,
be.result.growth$ci.ub, be.result.life$ci.ub,
be.result.hm.p$ci.ub),
P = c(be.result.pt$pval, be.result.d$pval,
be.result.growth$pval, be.result.life$pval,
be.result.hm.p$pval))
# Combine table
be_result <- full_join(dat_sum_be, be.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(be_result, digits = 3,
caption = "Summary of predict variables of below-ground biomass",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 21 | 1 | -0.630 | 0.578 | -1.763 | 0.503 | 0.276 |
| Dormancy | MPD | 21 | 1 | -0.803 | 0.578 | -1.935 | 0.329 | 0.164 |
| Dormancy | ND | 285 | 24 | -0.408 | 0.157 | -0.714 | -0.101 | 0.009 |
| Dormancy | PD | 402 | 29 | -0.860 | 0.141 | -1.135 | -0.584 | 0.000 |
| Dormancy | PY | 238 | 18 | -0.975 | 0.190 | -1.348 | -0.602 | 0.000 |
| Growth form | Herbaceous | 539 | 32 | -0.915 | 0.142 | -1.193 | -0.637 | 0.000 |
| Growth form | Woody | 428 | 36 | -0.542 | 0.139 | -0.815 | -0.269 | 0.000 |
| Life span | Long-lived | 746 | 58 | -0.721 | 0.112 | -0.940 | -0.502 | 0.000 |
| Life span | Short-lived | 221 | 12 | -0.744 | 0.204 | -1.144 | -0.345 | 0.000 |
| Pollution | Emerging | 23 | 4 | 0.111 | 0.418 | -0.708 | 0.930 | 0.791 |
| Pollution | Inorganic | 114 | 5 | -0.360 | 0.348 | -1.042 | 0.322 | 0.301 |
| Pollution | Metal | 758 | 50 | -0.863 | 0.117 | -1.091 | -0.634 | 0.000 |
| Pollution | Organic | 72 | 7 | -0.395 | 0.313 | -1.009 | 0.219 | 0.208 |
| Heavy metal | As | 32 | 4 | -0.986 | 0.289 | -1.552 | -0.420 | 0.001 |
| Heavy metal | Cd | 85 | 17 | -0.635 | 0.153 | -0.935 | -0.335 | 0.000 |
| Heavy metal | Cr | 8 | 2 | -0.409 | 0.326 | -1.048 | 0.230 | 0.210 |
| Heavy metal | Hg | 22 | 3 | 0.507 | 0.234 | 0.049 | 0.965 | 0.030 |
| Heavy metal | Pb | 119 | 18 | -0.663 | 0.152 | -0.960 | -0.365 | 0.000 |
| Heavy metal | others | 385 | 20 | -0.971 | 0.138 | -1.241 | -0.701 | 0.000 |
# Save data
write_csv(be_result, "table/be_result_predict variables.csv")Below-ground biomass comparisons
# Pollution
anova(be.result.pt, L = c(0, 1, -1, 0)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.5029 0.3669 1.3704 0.1706
##
## Test of Hypothesis:
## QM(df = 1) = 1.8780, p-val = 0.1706anova(be.result.pt, L = c(0, 1, 0, -1)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.0346 0.4682 0.0740 0.9410
##
## Test of Hypothesis:
## QM(df = 1) = 0.0055, p-val = 0.9410anova(be.result.pt, L = c(0, 0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.4682 0.3339 -1.4022 0.1609
##
## Test of Hypothesis:
## QM(df = 1) = 1.9662, p-val = 0.1609# Growth form
anova(be.result.growth, L = c(1, -1)) ##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.3726 0.1880 -1.9820 0.0475
##
## Test of Hypothesis:
## QM(df = 1) = 3.9282, p-val = 0.0475# Life span
anova(be.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: 0.0233 0.2037 0.1146 0.9087
##
## Test of Hypothesis:
## QM(df = 1) = 0.0131, p-val = 0.9087# Dormancy
anova(be.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: 0.4521 0.1887 2.3953 0.0166
##
## Test of Hypothesis:
## QM(df = 1) = 5.7375, p-val = 0.0166anova(be.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.5676 0.2404 2.3608 0.0182
##
## Test of Hypothesis:
## QM(df = 1) = 5.5732, p-val = 0.0182anova(be.result.d, L = c(0, 0, 0, 1, -1)) ##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.1155 0.2173 0.5318 0.5948
##
## Test of Hypothesis:
## QM(df = 1) = 0.2829, p-val = 0.5948Seedling total biomass
sl.result.pt <- rma.mv(yi, v.sl, mods = ~Pollution - 1, data = dat_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.result.growth <- rma.mv(yi, v.sl, mods = ~GrowthForm - 1, data = dat_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.result.life <- rma.mv(yi, v.sl, mods = ~LifeSpan - 1, data = dat_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.result.d <- rma.mv(yi, v.sl, mods = ~Dormancy - 1, data = dat_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.result.hm.p <- rma.mv(yi, v.hm.sl, mods = ~Pollutant - 1, data = dat_hm_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for result
dat_sum_sl <- dat_sl %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_sl %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Seedling total biomass meta analysis result table
sl.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Heavy metal"),
c(4, 5, 2, 2, 4)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"Cd", "Hg", "others", "Pb"),
Mean = c(sl.result.pt$b, sl.result.d$b,
sl.result.growth$b, sl.result.life$b,
sl.result.hm.p$b),
SE = c(sl.result.pt$se, sl.result.d$se,
sl.result.growth$se, sl.result.life$se,
sl.result.hm.p$se),
Lower_95_CI = c(sl.result.pt$ci.lb, sl.result.d$ci.lb,
sl.result.growth$ci.lb, sl.result.life$ci.lb,
sl.result.hm.p$ci.lb),
Upper_95_CI = c(sl.result.pt$ci.ub, sl.result.d$ci.ub,
sl.result.growth$ci.ub, sl.result.life$ci.ub,
sl.result.hm.p$ci.ub),
P = c(sl.result.pt$pval, sl.result.d$pval,
sl.result.growth$pval, sl.result.life$pval,
sl.result.hm.p$pval))
# Combine table
sl_result <- full_join(dat_sum_sl, sl.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(sl_result, digits = 3,
caption = "Summary of predict variables of seedling total biomass",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 21 | 1 | -0.533 | 0.315 | -1.149 | 0.083 | 0.090 |
| Dormancy | MPD | 21 | 1 | -0.344 | 0.314 | -0.960 | 0.273 | 0.274 |
| Dormancy | ND | 98 | 7 | -0.507 | 0.213 | -0.924 | -0.090 | 0.017 |
| Dormancy | PD | 89 | 9 | -0.109 | 0.205 | -0.511 | 0.292 | 0.593 |
| Dormancy | PY | 84 | 8 | -0.642 | 0.241 | -1.114 | -0.170 | 0.008 |
| Growth form | Herbaceous | 179 | 8 | -0.597 | 0.245 | -1.077 | -0.116 | 0.015 |
| Growth form | Woody | 134 | 14 | -0.296 | 0.185 | -0.658 | 0.067 | 0.110 |
| Life span | Long-lived | 209 | 18 | -0.366 | 0.158 | -0.675 | -0.057 | 0.020 |
| Life span | Short-lived | 104 | 6 | -0.533 | 0.230 | -0.985 | -0.081 | 0.021 |
| Pollution | Emerging | 29 | 4 | 0.020 | 0.300 | -0.568 | 0.608 | 0.946 |
| Pollution | Inorganic | 22 | 1 | -0.139 | 0.635 | -1.383 | 1.106 | 0.827 |
| Pollution | Metal | 244 | 15 | -0.563 | 0.169 | -0.894 | -0.233 | 0.001 |
| Pollution | Organic | 18 | 3 | -0.153 | 0.401 | -0.939 | 0.632 | 0.703 |
| Heavy metal | Cd | 44 | 6 | -0.408 | 0.095 | -0.593 | -0.222 | 0.000 |
| Heavy metal | Hg | 9 | 1 | -0.275 | 0.156 | -0.581 | 0.032 | 0.079 |
| Heavy metal | Pb | 46 | 7 | -0.216 | 0.092 | -0.397 | -0.035 | 0.019 |
| Heavy metal | others | 123 | 6 | -0.351 | 0.082 | -0.512 | -0.190 | 0.000 |
# Save data
write_csv(sl_result, "table/sl_result_predict variables.csv")Seedling total biomass comparisons
# Pollution
anova(sl.result.pt, L = c(0, 1, -1, 0))##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.4247 0.6572 0.6463 0.5181
##
## Test of Hypothesis:
## QM(df = 1) = 0.4177, p-val = 0.5181anova(sl.result.pt, L = c(0, 1, 0, -1))##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.0145 0.7510 0.0193 0.9846
##
## Test of Hypothesis:
## QM(df = 1) = 0.0004, p-val = 0.9846anova(sl.result.pt, L = c(0, 0, 1, -1))##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.4102 0.4348 -0.9435 0.3454
##
## Test of Hypothesis:
## QM(df = 1) = 0.8901, p-val = 0.3454anova(sl.result.pt, L = c(1, -1, 0, 0))##
## Hypothesis:
## 1: PollutionEmerging - PollutionInorganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.1589 0.7025 0.2262 0.8211
##
## Test of Hypothesis:
## QM(df = 1) = 0.0512, p-val = 0.8211anova(sl.result.pt, L = c(1, 0, -1, 0))##
## Hypothesis:
## 1: PollutionEmerging - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.5836 0.3092 1.8877 0.0591
##
## Test of Hypothesis:
## QM(df = 1) = 3.5634, p-val = 0.0591anova(sl.result.pt, L = c(1, 0, 0, -1))##
## Hypothesis:
## 1: PollutionEmerging - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: 0.1734 0.5007 0.3463 0.7291
##
## Test of Hypothesis:
## QM(df = 1) = 0.1199, p-val = 0.7291# Growth form
anova(sl.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.3012 0.3072 -0.9803 0.3269
##
## Test of Hypothesis:
## QM(df = 1) = 0.9611, p-val = 0.3269# Life span
anova(sl.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: 0.1667 0.2302 0.7240 0.4690
##
## Test of Hypothesis:
## QM(df = 1) = 0.5242, p-val = 0.4690# Dormancy
anova(sl.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: -0.3976 0.2062 -1.9282 0.0538
##
## Test of Hypothesis:
## QM(df = 1) = 3.7180, p-val = 0.0538anova(sl.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.1350 0.2926 0.4614 0.6445
##
## Test of Hypothesis:
## QM(df = 1) = 0.2129, p-val = 0.6445anova(sl.result.d, L = c(0, 0, 0, 1, -1))##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.5326 0.3019 1.7641 0.0777
##
## Test of Hypothesis:
## QM(df = 1) = 3.1122, p-val = 0.07773.5.2.5 Biomass allocation
ba.result.pt <- rma.mv(yi, v.ba, mods = ~Pollution - 1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.result.growth <- rma.mv(yi, v.ba, mods = ~GrowthForm - 1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.result.life <- rma.mv(yi, v.ba, mods = ~LifeSpan - 1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.result.d <- rma.mv(yi, v.ba, mods = ~Dormancy - 1, data = dat_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
ba.result.hm.p <- rma.mv(yi, v.hm.ba, mods = ~Pollutant - 1, data = dat_hm_ba,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
# Summary for biomass allocation
dat_sum_ba <- dat_ba %>%
pivot_longer(cols = c(Pollution, Dormancy, GrowthForm, LifeSpan),
names_to = "Predictor", values_to = "Predictortype") %>%
mutate(Predictor = case_when(
Predictor == "GrowthForm" ~ "Growth form",
Predictor == "LifeSpan" ~ "Life span",
TRUE ~ Predictor
)) %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
rbind(dat_hm_ba %>%
pivot_longer(cols = c(Pollutant),
names_to = "Predictor", values_to = "Predictortype") %>%
group_by(Predictor, Predictortype) %>%
summarise(Effectsize_count = n_distinct(Effectsize_ID),
Reference_count = n_distinct(Ref_ID),
.groups = "drop") %>%
mutate(Predictor = case_when(
Predictor == "Pollutant" ~ "Heavy metal",
TRUE ~ Predictor
)))
# Biomass allocation meta analysis result table
ba.result.table <- data.frame(Predictor = rep(c("Pollution", "Dormancy", "Growth form",
"Life span", "Heavy metal"),
c(4, 5, 2, 2, 6)),
Predictortype = c("Emerging", "Inorganic", "Metal", "Organic",
"MD", "MPD", "ND", "PD", "PY",
"Herbaceous", "Woody",
"Long-lived", "Short-lived",
"As", "Cd", "Cr", "Hg", "others", "Pb"),
Mean = c(ba.result.pt$b, ba.result.d$b,
ba.result.growth$b, ba.result.life$b,
ba.result.hm.p$b),
SE = c(ba.result.pt$se, ba.result.d$se,
ba.result.growth$se, ba.result.life$se,
ba.result.hm.p$se),
Lower_95_CI = c(ba.result.pt$ci.lb, ba.result.d$ci.lb,
ba.result.growth$ci.lb, ba.result.life$ci.lb,
ba.result.hm.p$ci.lb),
Upper_95_CI = c(ba.result.pt$ci.ub, ba.result.d$ci.ub,
ba.result.growth$ci.ub, ba.result.life$ci.ub,
ba.result.hm.p$ci.ub),
P = c(ba.result.pt$pval, ba.result.d$pval,
ba.result.growth$pval, ba.result.life$pval,
ba.result.hm.p$pval))
# Combine table
ba_result <- full_join(dat_sum_ba, ba.result.table, by = c("Predictortype", "Predictor")) %>%
select(Predictor, everything())
knitr::kable(ba_result, digits = 3,
caption = "Summary of predict variables of biomass allocation for differnet levels",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| Predictor | Predictortype | Effectsize_count | Reference_count | Mean | SE | Lower_95_CI | Upper_95_CI | P |
|---|---|---|---|---|---|---|---|---|
| Dormancy | MD | 21 | 1 | -0.315 | 0.668 | -1.624 | 0.994 | 0.637 |
| Dormancy | MPD | 21 | 1 | -0.515 | 0.672 | -1.833 | 0.803 | 0.444 |
| Dormancy | ND | 227 | 18 | -0.083 | 0.170 | -0.416 | 0.250 | 0.626 |
| Dormancy | PD | 318 | 20 | -0.578 | 0.159 | -0.890 | -0.267 | 0.000 |
| Dormancy | PY | 172 | 12 | -0.246 | 0.218 | -0.674 | 0.182 | 0.259 |
| Growth form | Herbaceous | 442 | 24 | -0.473 | 0.147 | -0.760 | -0.185 | 0.001 |
| Growth form | Woody | 317 | 24 | -0.149 | 0.153 | -0.450 | 0.151 | 0.330 |
| Life span | Long-lived | 556 | 39 | -0.289 | 0.121 | -0.527 | -0.052 | 0.017 |
| Life span | Short-lived | 203 | 10 | -0.421 | 0.221 | -0.855 | 0.013 | 0.057 |
| Pollution | Emerging | 19 | 3 | 0.044 | 0.410 | -0.760 | 0.848 | 0.914 |
| Pollution | Inorganic | 92 | 4 | -0.149 | 0.323 | -0.782 | 0.484 | 0.644 |
| Pollution | Metal | 591 | 34 | -0.459 | 0.125 | -0.704 | -0.215 | 0.000 |
| Pollution | Organic | 57 | 6 | 0.164 | 0.287 | -0.398 | 0.726 | 0.567 |
| Heavy metal | As | 26 | 3 | -0.402 | 0.411 | -1.206 | 0.403 | 0.328 |
| Heavy metal | Cd | 46 | 11 | -0.289 | 0.183 | -0.648 | 0.070 | 0.115 |
| Heavy metal | Cr | 8 | 2 | -0.053 | 0.288 | -0.618 | 0.512 | 0.854 |
| Heavy metal | Hg | 9 | 1 | 0.057 | 0.280 | -0.491 | 0.605 | 0.839 |
| Heavy metal | Pb | 70 | 12 | -0.220 | 0.182 | -0.577 | 0.138 | 0.228 |
| Heavy metal | others | 333 | 15 | -0.786 | 0.156 | -1.091 | -0.480 | 0.000 |
# Save data
write_csv(ba_result, "table/ba_result_predict variables.csv")Biomass allocation comparisons
# Pollutant type
anova(ba.result.pt, L = c(0, 1, -1, 0)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionMetal = 0
##
## Results:
## estimate se zval pval
## 1: 0.3102 0.3463 0.8959 0.3703
##
## Test of Hypothesis:
## QM(df = 1) = 0.8026, p-val = 0.3703anova(ba.result.pt, L = c(0, 1, 0, -1)) ##
## Hypothesis:
## 1: PollutionInorganic - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.3134 0.4319 -0.7255 0.4682
##
## Test of Hypothesis:
## QM(df = 1) = 0.5263, p-val = 0.4682anova(ba.result.pt, L = c(0, 0, 1, -1)) ##
## Hypothesis:
## 1: PollutionMetal - PollutionOrganic = 0
##
## Results:
## estimate se zval pval
## 1: -0.6236 0.3110 -2.0048 0.0450
##
## Test of Hypothesis:
## QM(df = 1) = 4.0190, p-val = 0.0450# Growth form
anova(ba.result.growth, L = c(1, -1))##
## Hypothesis:
## 1: GrowthFormHerbaceous - GrowthFormWoody = 0
##
## Results:
## estimate se zval pval
## 1: -0.3234 0.2072 -1.5608 0.1186
##
## Test of Hypothesis:
## QM(df = 1) = 2.4360, p-val = 0.1186# Adult lifespan
anova(ba.result.life, L = c(1, -1))##
## Hypothesis:
## 1: LifeSpanLong-lived - LifeSpanShort-lived = 0
##
## Results:
## estimate se zval pval
## 1: 0.1319 0.2377 0.5548 0.5790
##
## Test of Hypothesis:
## QM(df = 1) = 0.3078, p-val = 0.5790# Dormancy
anova(ba.result.d, L = c(0, 0, 1, -1, 0)) ##
## Hypothesis:
## 1: DormancyND - DormancyPD = 0
##
## Results:
## estimate se zval pval
## 1: 0.4955 0.2219 2.2329 0.0256
##
## Test of Hypothesis:
## QM(df = 1) = 4.9856, p-val = 0.0256anova(ba.result.d, L = c(0, 0, 1, 0, -1)) ##
## Hypothesis:
## 1: DormancyND - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: 0.1634 0.2753 0.5934 0.5529
##
## Test of Hypothesis:
## QM(df = 1) = 0.3521, p-val = 0.5529anova(ba.result.d, L = c(0, 0, 0, 1, -1)) ##
## Hypothesis:
## 1: DormancyPD - DormancyPY = 0
##
## Results:
## estimate se zval pval
## 1: -0.3322 0.2548 -1.3037 0.1923
##
## Test of Hypothesis:
## QM(df = 1) = 1.6996, p-val = 0.19233.5.3 Continuous variables
3.5.3.1 Regeneration aspects
Seed mass
logmass.slope <- data.frame(AspectType = c("GP", "GS", "SG", "BA"),
Intercept = c(gp.qm.result.mass$b[1], gs.qm.result.mass$b[1],
sg.qm.result.mass$b[1], ba.qm.result.mass$b[1]),
Slop = c(gp.qm.result.mass$b[2], gs.qm.result.mass$b[2],
sg.qm.result.mass$b[2], ba.qm.result.mass$b[2]),
Slop_p = c(gp.qm.result.mass$pval[2], gs.qm.result.mass$pval[2],
sg.qm.result.mass$pval[2], ba.qm.result.mass$pval[2]))
knitr::kable(logmass.slope, digits = 3,
caption = "Summary of seed mass",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| AspectType | Intercept | Slop | Slop_p |
|---|---|---|---|
| GP | -0.457 | 0.082 | 0.056 |
| GS | -0.271 | -0.052 | 0.589 |
| SG | -0.521 | -0.025 | 0.628 |
| BA | -0.413 | 0.147 | 0.103 |
Concentration
hm.con.slope <- data.frame(AspectType = c("GP", "GS", "SG", "BA"),
Intercept = c(gp.qm.result.con$b[1], gs.qm.result.con$b[1],
sg.qm.result.con$b[1], ba.qm.result.con$b[1]),
Slop = c(gp.qm.result.con$b[2], gs.qm.result.con$b[2],
sg.qm.result.con$b[2], ba.qm.result.con$b[2]),
Slop_p = c(gp.qm.result.con$pval[2], gs.qm.result.con$pval[2],
sg.qm.result.con$pval[2], ba.qm.result.con$pval[2]))
knitr::kable(hm.con.slope, digits = 3,
caption = "Summary of pollution concentration",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| AspectType | Intercept | Slop | Slop_p |
|---|---|---|---|
| GP | -0.275 | -0.224 | 0 |
| GS | -0.242 | -0.147 | 0 |
| SG | -0.437 | -0.397 | 0 |
| BA | -0.462 | -0.256 | 0 |
3.5.3.2 Seedling parts
Seed mass
ab.qm.result.mass <- rma.mv(yi, v.ab, mods = ~1 + Logmass, data = dat_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.qm.result.mass <- rma.mv(yi, v.be, mods = ~1 + Logmass, data = dat_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.qm.result.mass <- rma.mv(yi, v.sl, mods = ~1 + Logmass, data = dat_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
logmass.slope <- data.frame(AspectType = c("AB", "BE", "SL"),
Intercept = c(ab.qm.result.mass$b[1], be.qm.result.mass$b[1],
sl.qm.result.mass$b[1]),
Slop = c(ab.qm.result.mass$b[2], be.qm.result.mass$b[2],
sl.qm.result.mass$b[2]),
Slop_p = c(ab.qm.result.mass$pval[2], be.qm.result.mass$pval[2],
sl.qm.result.mass$pval[2]))
knitr::kable(logmass.slope, digits = 3,
caption = "Summary of seed mass of different parts of seedling",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| AspectType | Intercept | Slop | Slop_p |
|---|---|---|---|
| AB | -0.379 | -0.049 | 0.334 |
| BE | -0.736 | 0.017 | 0.834 |
| SL | -0.421 | 0.021 | 0.856 |
Concentration
ab.qm.result.con <- rma.mv(yi, v.hm.ab, mods = ~1 + Logcon, data = dat_hm_ab,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
be.qm.result.con <- rma.mv(yi, v.hm.be, mods = ~1 + Logcon, data = dat_hm_be,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
sl.qm.result.con <- rma.mv(yi, v.hm.sl, mods = ~1 + Logcon, data = dat_hm_sl,
random = list(~1|Ref_ID/Effectsize_ID, ~1|Species))
hm.con.slope <- data.frame(AspectType = c("AB", "BE", "SL"),
Intercept = c(ab.qm.result.con$b[1], be.qm.result.con$b[1],
sl.qm.result.con$b[1]),
Slop = c(ab.qm.result.con$b[2], be.qm.result.con$b[2],
sl.qm.result.con$b[2]),
Slop_p = c(ab.qm.result.con$pval[2], be.qm.result.con$pval[2],
sl.qm.result.con$pval[2]))
knitr::kable(hm.con.slope, digits = 3,
caption = "Summary of pollution concentration of different parts of seedling",
full_width = T, font_size = 20, align = "c", "html") %>%
kable_styling(position = "center")| AspectType | Intercept | Slop | Slop_p |
|---|---|---|---|
| AB | -0.268 | -0.226 | 0 |
| BE | -0.558 | -0.527 | 0 |
| SL | 0.162 | -0.485 | 0 |
3.5.4 Figure for predictors in regeneration
3.5.4.1 Figure for plant characteristics
Germination percentage
# Data
gp.fig.pc.data <- filter(gp_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
gp.fig.pc.data$Predictortype <- factor(gp.fig.pc.data$Predictortype, levels = Predictortype.2_level)
gp.pc.fig <- ggplot(data = gp.fig.pc.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MPD" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Woody" = 16, "Herbaceous" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#9F79EE", "PD" = "#9F79EE", "PY" = "#9F79EE",
"Woody" = "#9F79EE", "Herbaceous" = "#9F79EE",
"Long-lived" = "#9F79EE", "Short-lived" = "#9F79EE")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#9F79EE", "PD" = "#9F79EE", "PY" = "#9F79EE",
"Woody" = "#9F79EE", "Herbaceous" = "#9F79EE",
"Long-lived" = "#9F79EE", "Short-lived" = "#9F79EE")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult lifespan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Germination percentage") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.y = element_text(size = 39, colour = "black", family = "serif"),
axis.text.x = element_text(size = 39, colour = "black", family = "serif"),
plot.title = element_text(size = 39, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[3])) + # ND
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[4])) + # PD
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[5])) + # PY
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[1])) + # MD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[2])) + # MPD
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[9])) + # Short
annotate("text", x = 8, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[8])) + # Long
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[7])) + # Herb
annotate("text", x = 11, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pc.data$Effectsize_count[6])) + # Woody
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
gp.pc.fig
Germination speed
# Data
gs.fig.pc.data <- filter(gs_result, Predictor %in% c("Dormancy", "Growth form"))
gs.fig.pc.data$Predictortype <- factor(gs.fig.pc.data$Predictortype, levels = Predictortype.2_level)
gs.pc.fig <- ggplot(data = gs.fig.pc.data, mapping = aes(x= Predictortype, y = Mean),
fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("ND" = 16, "PD" = 16, "PY" = 16,
"Herbaceous" = 16, "Woody" = 16),
color = c("ND" = "#00B2EE", "PD" = "#00B2EE", "PY" = "#ABABAB",
"Herbaceous" = "#ABABAB", "Woody" = "#00B2EE")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("ND" = "#00B2EE", "PD" = "#00B2EE", "PY" = "#ABABAB",
"Herbaceous" = "#ABABAB", "Woody" = "#00B2EE")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
labs(title = "Germination speed") +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MD", "MPD", "",
"Long-lived", "Short-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MD", "MPD", "",
"Long-lived", "Short-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "", "", "Dormancy class",
"", "", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 39, colour = "black", family = "serif"),
axis.text.y = element_blank(),
plot.title = element_text(size = 39, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
guides(fill = "none", color = "none") +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(gs.fig.pc.data$Effectsize_count[1])) + # ND
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(gs.fig.pc.data$Effectsize_count[3])) + # PY
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(gs.fig.pc.data$Effectsize_count[2])) + # PD
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(gs.fig.pc.data$Effectsize_count[5])) + # Woody
annotate("text", x = 11, y = -2, family = "serif", size = 14,
label = paste(gs.fig.pc.data$Effectsize_count[4])) + # Herb
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N") +
coord_flip()
gs.pc.fig
Seedling growth
# Data
sg.fig.pc.data <- filter(sg_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
sg.fig.pc.data$Predictortype <- factor(sg.fig.pc.data$Predictortype, levels = Predictortype.2_level)
sg.pc.fig <- ggplot(data = sg.fig.pc.data, mapping = aes(x= Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MPD" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Herbaceous" = 16, "Woody" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#43CD80", "MPD" = "#43CD80", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#43CD80", "MPD" = "#43CD80", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Seedling growth") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 39, colour = "black", family = "serif"),
axis.text.y = element_blank(),
plot.title = element_text(size = 39, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
bracketsGrob(x1 = 0.7, y1 = 1.6/12, x2 = 0.7, y2 = 0.5/12,
h = 0.03, lwd = 5, col = "black", type = 4, ticks = NA) %>% annotation_custom +
annotate("text", x = 1.25, y = 0.5, label = "*", size = 30) +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[3])) + #ND
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[5])) + #PY
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[4])) + #PD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[1])) + #MD
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[2])) + #MPD
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[9])) +#short
annotate("text", x = 8, y = -1.9, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[8])) +#long
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[7])) + #woody
annotate("text", x = 11, y = -1.9, family = "serif", size = 14,
label = paste(sg.fig.pc.data$Effectsize_count[6])) + #herb
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
sg.pc.fig
Biomass allocation
# Data
ba.fig.pc.data <- filter(ba_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
ba.fig.pc.data$Predictortype <- factor(ba.fig.pc.data$Predictortype, levels = Predictortype.2_level)
ba.pc.fig <- ggplot(data = ba.fig.pc.data, mapping = aes(x= Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MDP" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Herbaceous" = 16, "Woody" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#ABABAB", "MDP" = "#ABABAB", "ND" = "#ABABAB", "PD" = "#FF8C69", "PY" = "#ABABAB",
"Herbaceous" = "#FF8C69", "Woody" = "#ABABAB",
"Long-lived" = "#FF8C69", "Short-lived" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#ABABAB", "PD" = "#FF8C69", "PY" = "#ABABAB",
"Herbaceous" = "#FF8C69", "Woody" = "#ABABAB",
"Long-lived" = "#FF8C69", "Short-lived" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.05), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Biomass allocation") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 39, colour = "black", family = "serif"),
axis.text.y = element_blank(),
plot.title = element_text(size = 39, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
bracketsGrob(x1 = 0.8, y1 = 2.7/12, x2 = 0.8, y2 = 0.5/12,
h = 0.03, lwd = 5, col = "black", type = 4, ticks = NA) %>% annotation_custom +
annotate("text", x = 1.9, y = 0.9, label = "*", size = 30) +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[3])) + # ND
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[4])) + # PD
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[5])) + # PY
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[1])) + # MD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[2])) + # MPD
annotate("text", x = 11, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[6])) + # Herb
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[7])) + # Woody
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[9])) + # Short
annotate("text", x = 8, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pc.data$Effectsize_count[8])) + # Long
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
ba.pc.fig
Combine plant characteristics figures
library(patchwork)
# Data
plant.fig <- gp.pc.fig + gs.pc.fig + sg.pc.fig + ba.pc.fig +
plot_layout(nrow = 1, byrow = TRUE, widths = c(2, 2, 2, 2), heights = c(3, 2)) +
plot_annotation(caption = "Effect size (Log response ratio)",
theme = theme(plot.caption = element_text(size = 50, family = "serif", hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/plant.bind.jpg", width = 35, height = 12, units = "in", dpi = 300)
plant.fig
seed mass
Germination percentage
# Data
gp.logmass.preds <- data.frame(AspectType = "GP",
predict(gp.qm.result.mass, addx = TRUE))
gp.mass.fig <- ggplot() +
geom_point(dat_gp,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#9F79EE", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = gp.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = gp.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.7, 3.01), breaks = seq(-6, 3, 1.5),
name = "Log response ratio") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 32, colour = "black", family = "serif"),
axis.text.y = element_text(size = 32, colour = "black", family = "serif"),
axis.title.y = element_text(size = 32, colour = "black", family = "serif",
margin = margin(r = 20)),
axis.title = element_text(size = 32, colour = "black", family = "serif"))
gp.mass.fig
Germination speed
# Data
gs.logmass.preds <- data.frame(AspectType = "GS",
predict(gs.qm.result.mass, addx = TRUE))
gs.mass.fig <- ggplot() +
geom_point(dat_gs,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#00B2EE", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = gs.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = gs.logmass.preds,
mapping = aes(x = X.Logmass, y = pred),
linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.7, 3.01), breaks = seq(-6, 3, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 32, colour = "black", family = "serif"),
axis.text.y = element_blank())
gs.mass.fig
Seedling growth
# Data
sg.logmass.preds <- data.frame(AspectType = "SG",
predict(sg.qm.result.mass, addx = TRUE))
sg.mass.fig <- ggplot() +
geom_point(dat_sg,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.45, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = sg.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = sg.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.7, 3.01), breaks = seq(-6, 3, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 32, colour = "black", family = "serif"),
axis.text.y = element_blank())
sg.mass.fig
Biomass allocation
# Data
ba.logmass.preds <- data.frame(AspectType = "BA",
predict(ba.qm.result.mass, addx = TRUE))
ba.mass.fig <- ggplot() +
geom_point(dat_ba,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.45, shape = 1, color = "#FF8C69", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = ba.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = ba.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.7, 3.01), breaks = seq(-6, 3, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 32, colour = "black", family = "serif"),
axis.text.y = element_blank())
ba.mass.fig
Combine seed mass figures
mass.fig <- gp.mass.fig + gs.mass.fig + sg.mass.fig + ba.mass.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2, 2), heights = c(2, 2)) +
plot_annotation(caption = "Seed mass (mg) [log scale]",
theme = theme(plot.caption = element_text(size = 33, family = "serif", hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/smass.bind.jpg", width = 24, height = 6, units = "in", dpi = 300)
mass.fig
3.5.4.2 Figure for pollutant
Germination percentage
# Data
gp.fig.pt.data <- filter(gp_result, Predictor %in% "Pollution")
gp.fig.pt.data$Predictortype <- factor(gp.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
gp.pt.fig <- ggplot(data = gp.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10,
shape = c("Emerging" = 16, "Inorganic" = 16, "Organic" = 16, "Metal" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#9F79EE", "Organic" = "#9F79EE", "Metal" = "#9F79EE")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#9F79EE", "Organic" = "#9F79EE", "Metal" = "#9F79EE")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
labs(title = "Germination percentage") +
theme(axis.title.x.bottom = element_blank(),
axis.text.y = element_text(size = 37, colour = "black", family = "serif"),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold", colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pt.data$Effectsize_count[1])) + # Emerging
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pt.data$Effectsize_count[4])) + # Organic
annotate("text", x = 3, y = -2, family = "serif", size = 14,
label = paste(gp.fig.pt.data$Effectsize_count[2])) + # Inorganic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(gp.fig.pt.data$Effectsize_count[3])) + # Metal
coord_flip() +
guides(fill = "none", color = "none")
gp.pt.fig
Germination speed
# Data
gs.fig.pt.data <- filter(gs_result, Predictor %in% "Pollution")
gs.fig.pt.data$Predictortype <- factor(gs.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
gs.pt.fig <- ggplot(data = gs.fig.pt.data, mapping = aes(x = Predictortype, y = Mean),
fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Inorganic" = 16, "Metal" = 16, "Organic" = 16),
color = c("Inorganic" = "#ABABAB", "Metal" = "#00B2EE", "Organic" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Inorganic" = "#ABABAB", "Metal" = "#00B2EE", "Organic" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
scale_x_discrete(limits = rev(Predictortype.1_level),
breaks = c("Metal", "Organic", "Inorganic", "Emerging"),
labels = c("Metal", "Organic", "Inorganic", "")) +
labs(title = "Germination speed") +
theme(plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
theme_plot +
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(gs.fig.pt.data$Effectsize_count[3])) + #Organic
annotate("text", x = 3, y = -2, family = "serif", size = 14,
label = paste(gs.fig.pt.data$Effectsize_count[1])) + #Inorganic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(gs.fig.pt.data$Effectsize_count[2])) + #Metal
guides(fill = "none", color = "none") +
coord_flip()
gs.pt.fig
Seedling growth
# Data
sg.fig.pt.data <- filter(sg_result, Predictor %in% "Pollution")
sg.fig.pt.data$Predictortype <- factor(sg.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
sg.pt.fig <- ggplot(data = sg.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Emerging" = 16, "Inorganic" = 16, "Metal" = 16, "Organic" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#43CD80")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#43CD80")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
labs(title = "Seedling growth") +
theme(axis.text.y = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(sg.fig.pt.data$Effectsize_count[1])) + #Emerging
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pt.data$Effectsize_count[4])) + #Organic
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(sg.fig.pt.data$Effectsize_count[2])) + #Inorganic
annotate("text", x = 4, y = -1.9, family = "serif", size = 14,
label = paste(sg.fig.pt.data$Effectsize_count[3])) + #Metal
coord_flip() +
guides(fill = "none", color = "none")
sg.pt.fig
Biomass allocation
# Data
ba.fig.pt.data <- filter(ba_result, Predictor %in% "Pollution")
ba.fig.pt.data$Predictortype <- factor(ba.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
ba.pt.fig <- ggplot(data = ba.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Emerging" = 16, "Inorganic" = 16, "Organic" = 16, "Metal" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#FF8C69", "Organic" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#FF8C69", "Organic" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
labs(title = "Biomass allocation") +
theme(axis.text.y = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(ba.fig.pt.data$Effectsize_count[1])) + # Emerging
annotate("text", x = 3, y = -2, family = "serif", size = 14,
label = paste(ba.fig.pt.data$Effectsize_count[2])) + # Inorganic
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(ba.fig.pt.data$Effectsize_count[4])) + # Organic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(ba.fig.pt.data$Effectsize_count[3])) + # Metal
coord_flip() +
guides(fill = "none", color = "none")
ba.pt.fig
Combine pollutant type figures
poll.fig <- gp.pt.fig + gs.pt.fig + sg.pt.fig + ba.pt.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2, 2), heights = c(3, 2)) +
plot_annotation(caption = "Effect size (Log response ratio)",
theme = theme(plot.caption = element_text(size = 45, family = "serif",
hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/poll.bind.jpg", width = 28, height = 9, units = "in", dpi = 300)
poll.fig
Heavy metal concentration
Germination percentage
# Data
gp.con.data <- data.frame(AspectType = "GP",
predict(gp.qm.result.con, addx = TRUE))
gp.con.fig <- ggplot() +
geom_point(dat_hm_gp,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#9F79EE", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = gp.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
geom_line(data = gp.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.7, 2.2), breaks = seq(-6, 1.5, 1.5),
name = "Log response ratio") +
theme_plot +
theme(axis.text.y = element_text(size = 35, colour = "black", family = "serif"),
axis.title.y = element_text(size = 35, colour = "black", family = "serif",
margin = margin(r = 20)),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.title.x.bottom = element_blank()) +
annotate("text", x = -1.2, y = -5, family = "serif", size = 12,
label = "slope = -0.224") + #slope value
annotate("text", x = -0.8, y = -6, family = "serif", size = 12,
label = "p < 0.001") #slope p
gp.con.fig
Germination speed
# Data
gs.con.data <- data.frame(AspectType = "GS",
predict(gs.qm.result.con, addx = TRUE))
gs.con.fig <- ggplot() +
geom_point(dat_hm_gs,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#00B2EE", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = gs.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
geom_line(data = gs.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.7, 2.2), breaks = seq(-6, 1.5, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
annotate("text", x = -1.2, y = -5, family = "serif", size = 12,
label = "slope = -0.147") + #slope value
annotate("text", x = -0.8, y = -6, family = "serif", size = 12,
label = "p < 0.001") #slope p
gs.con.fig
Seedling growth
# Data
sg.con.data <- data.frame(AspectType = "SG",
predict(sg.qm.result.con, addx = TRUE))
sg.con.fig <- ggplot() +
geom_point(dat_hm_sg,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_line(data = sg.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
geom_ribbon(data = sg.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.7, 2.2), breaks = seq(-6, 1.5, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
annotate("text", x = -1.2, y = -5, family = "serif", size = 12,
label = "slope = -0.397") + #slope value
annotate("text", x = -0.8, y = -6, family = "serif", size = 12,
label = "p < 0.001") #slope p
sg.con.fig
Biomass allocation
# Data
ba.con.data <- data.frame(AspectType = "BA",
predict(ba.qm.result.con, addx = TRUE))
ba.con.fig <- ggplot() +
geom_point(dat_hm_ba,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#FF8C69", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_line(data = ba.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
geom_ribbon(data = ba.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.7, 2.2), breaks = seq(-6, 1.5, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
annotate("text", x = -1.2, y = -5, family = "serif", size = 12,
label = "slope = -0.256") + #slope value
annotate("text", x = -0.8, y = -6, family = "serif", size = 12,
label = "p < 0.001") #slope p
ba.con.fig
Combine pollution concentration figures
hm.fig <- gp.con.fig + gs.con.fig + sg.con.fig + ba.con.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2, 2), heights = c(2, 2)) +
plot_annotation(caption = "Concentration of metals (mg/L) [log scale]",
theme = theme(plot.caption = element_text(size = 40, family = "serif",
hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/hm.bind.jpg", width = 23, height = 6, units = "in", dpi = 300)
hm.fig
3.5.5 Figure for predictors in seedling parts
3.5.5.1 Figure for seedling parts
# Data
sg.fig.data <- filter(metadata, AspectType %in% c("TB", "AB", "BE")) %>%
mutate(AspectType = case_when(
AspectType == "AB" ~ "Above-ground biomass",
AspectType == "BE" ~ "Below-ground biomass",
AspectType == "TB" ~ "Total biomass",
TRUE ~ AspectType))
sg.fig.data.p <- filter(sg_result, Predictor %in% "Parts")
sg.fig.data.p <- sg.fig.data.p %>% rename_at("Predictortype", ~"AspectType")
sg.fig.p <- ggplot(data = sg.fig.data,
mapping = aes(x = AspectType, y = yi, fill = AspectType)) +
geom_half_violin(aes(fill = AspectType, colour = AspectType),
position = position_nudge(x = 0),
linewidth = 0.5, alpha = 0.5, trim = TRUE, side = "r") +
geom_jitter(aes(x = stage(start = AspectType, after_scale = x - 0.3), color = AspectType),
width = 0.2, alpha = 0.3, shape = 1, stroke = 1.1) +
geom_pointrange(data = sg.fig.data.p,
mapping = aes(y = Mean, ymin = Lower_95_CI, ymax = Upper_95_CI),
position = position_nudge(x = 0), linewidth = 1.5, size = 0.7) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("#43CD80", "#43CD80", "#43CD80")) +
scale_color_manual(values = c("#43CD80", "#43CD80", "#43CD80")) +
scale_y_continuous(name = "Effect size (Log response ratio)",
labels = scales::number_format(accuracy = 0.1)) +
scale_x_discrete(breaks = c("Total biomass", "Above-ground biomass", "Below-ground biomass"),
labels = c("Total\n growth", "Above-ground\n growth",
"Below-ground\n growth")) +
theme(axis.title.x.bottom = element_text(size = 30, color = "black",
family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 30, colour = "black", family = "serif"),
axis.text.y = element_text(size = 30, color = "black", family = "serif")) +
theme_plot +
guides(fill = "none") +
coord_flip()
# Save figure
ggsave("figure/sg.fig.p.jpg", width = 8, height = 6, units = "in", dpi = 300)
sg.fig.p
3.5.5.2 Figure for plant characteristics
Above-ground biomass
# Data
ab.fig.pc.data <- filter(ab_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
ab.fig.pc.data$Predictortype <- factor(ab.fig.pc.data$Predictortype, levels = Predictortype.2_level)
ab.pc.fig <- ggplot(data = ab.fig.pc.data, mapping = aes(x= Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MPD" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Woody" = 16, "Herbaceous" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#ABABAB", "MPD" = "#43CD80", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#ABABAB", "MPD" = "#43CD80", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Above-ground biomass") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[3])) + # ND
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[5])) + # PY
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[4])) + # PD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[1])) + # MD
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[2])) + # MPD
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[9])) + # Short
annotate("text", x = 8, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[8])) + # Long
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[7])) + # Woody
annotate("text", x = 11, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pc.data$Effectsize_count[6])) + # Herb
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
ab.pc.fig
Below-ground biomass
# Data
be.fig.pc.data <- filter(be_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
be.fig.pc.data$Predictortype <- factor(be.fig.pc.data$Predictortype, levels = Predictortype.2_level)
be.pc.fig <- ggplot(data = be.fig.pc.data, mapping = aes(x= Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MPD" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Woody" = 16, "Herbaceous" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#43CD80", "PD" = "#43CD80", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#43CD80",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Below-ground biomass") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_blank(),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
bracketsGrob(x1 = 0.85, y1 = 2.7/12, x2 = 0.85, y2 = 0.5/12,
h = 0.03, lwd = 5, col = "black", type = 4, ticks = NA) %>% annotation_custom +
annotate("text", x = 1.8, y = 0.96, label = "*", size = 30) +
bracketsGrob(x1 = 0.7, y1 = 1.6/12, x2 = 0.7, y2 = 0.5/12,
h = 0.03, lwd = 5, col = "black", type = 4, ticks = NA) %>% annotation_custom +
annotate("text", x = 1.22, y = 0.45, label = "*", size = 30) +
bracketsGrob(x1 = 0.7, y1 = 10.6/12, x2 = 0.7, y2 = 9.6/12,
h = 0.03, lwd = 5, col = "black", type = 4, ticks = NA) %>% annotation_custom +
annotate("text", x = 10.4, y = 0.45, label = "*", size = 30) +
annotate("text", x = 1, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[3])) + # ND
annotate("text", x = 2, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[5])) + # PY
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[4])) + # PD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[1])) + # MD
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[2])) + # MPD
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[9])) + # Short
annotate("text", x = 8, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[8])) + # Long
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[7])) + # Woody
annotate("text", x = 11, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pc.data$Effectsize_count[6])) + # Herb
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
be.pc.fig
Seedling total biomass
# Data
sl.fig.pc.data <- filter(sl_result, Predictor %in% c("Dormancy", "Growth form", "Life span"))
sl.fig.pc.data$Predictortype <- factor(sl.fig.pc.data$Predictortype, levels = Predictortype.2_level)
sl.pc.fig <- ggplot(data = sl.fig.pc.data, mapping = aes(x= Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("MD" = 16, "MPD" = 16, "ND" = 16, "PD" = 16, "PY" = 16,
"Woody" = 16, "Herbaceous" = 16,
"Long-lived" = 16, "Short-lived" = 16),
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#43CD80", "PD" = "#ABABAB", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#ABABAB",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("MD" = "#ABABAB", "MPD" = "#ABABAB", "ND" = "#43CD80", "PD" = "#ABABAB", "PY" = "#43CD80",
"Herbaceous" = "#43CD80", "Woody" = "#ABABAB",
"Long-lived" = "#43CD80", "Short-lived" = "#43CD80")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_vline(xintercept = 6, linetype = 3, linewidth = 1.5) +
geom_vline(xintercept = 9, linetype = 3, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.02), breaks = seq(-2, 2, 1)) +
scale_x_discrete(limits = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
breaks = c("ND", "PY", "PD", "MPD", "MD", "",
"Short-lived", "Long-lived", " ",
"Woody", "Herbaceous", " "),
labels = c("ND", "PY", "PD", "MPD", "MD", "Dormancy class",
"Short-lived", "Long-lived", "Adult LifeSpan",
"Woody", "Herbaceous", "Growth form")) +
labs(title = "Seedling biomass") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_blank(),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
coord_flip() +
guides(fill = "none", color = "none") +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[3])) + # ND
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[5])) + # PY
annotate("text", x = 3, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[4])) + # PD
annotate("text", x = 4, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[1])) + # MD
annotate("text", x = 5, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[2])) + # MPD
annotate("text", x = 7, y = -1.95, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[9])) + # Short
annotate("text", x = 8, y = -1.95, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[8])) + # Long
annotate("text", x = 10, y = -1.95, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[7])) + # Woody
annotate("text", x = 11, y = -1.95, family = "serif", size = 14,
label = paste(sl.fig.pc.data$Effectsize_count[6])) + # Herb
annotate("text", x = 12, y = -2.05, family = "serif", size = 14,
label = "N")
sl.pc.fig
Combine plant characteristics figures
# Data
plant.seedling.fig <- ab.pc.fig + be.pc.fig + sl.pc.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2), heights = c(4, 2)) +
plot_annotation(caption = "Effect size (Log response ratio)",
theme = theme(plot.caption = element_text(size = 48, family = "serif", hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/plant.seedling.fig.jpg", width = 25, height = 11, units = "in", dpi = 300)
plant.seedling.fig
Seed mass
Above-ground biomass
# Data
ab.logmass.preds <- data.frame(AspectType = "AB",
predict(ab.qm.result.mass, addx = TRUE))
ab.mass.fig <- ggplot() +
geom_point(dat_ab,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.45, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = ab.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = ab.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.2, 1.5), breaks = seq(-6, 3, 1.5),
name = "Log response ratio") +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_text(size = 35, colour = "black",
family = "serif"),
axis.title.y = element_text(size = 35, colour = "black", family = "serif",
margin = margin(r = 20)))
ab.mass.fig
Below-ground biomass
# Data
be.logmass.preds <- data.frame(AspectType = "BE",
predict(be.qm.result.mass, addx = TRUE))
be.mass.fig <- ggplot() +
geom_point(dat_be,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.45, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = be.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = be.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.2, 1.6), breaks = seq(-6, 3, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank())
be.mass.fig
Seedling total biomass
# Data
sl.logmass.preds <- data.frame(AspectType = "SL",
predict(sl.qm.result.mass, addx = TRUE))
sl.mass.fig <- ggplot() +
geom_point(dat_sl,
mapping = aes(x = Logmass, y = yi, size = 1/vi),
alpha = 0.45, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_ribbon(data = sl.logmass.preds,
aes(x = X.Logmass, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "gray", alpha = 0.7) +
geom_line(data = sl.logmass.preds,
mapping = aes(x = X.Logmass, y = pred), linetype = "longdash", color = "#575757", linewidth = 1.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-2.02, 4.5), breaks = c(-2, 0, 2, 4),
labels = c("0.01", "1", "100", "10000")) +
scale_y_continuous(limits = c(-6.2, 1.6), breaks = seq(-6, 3, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank())
sl.mass.fig
Combine seed mass figures
mass.seedling.fig <- ab.mass.fig + be.mass.fig + sl.mass.fig +
plot_layout(nrow = 1, byrow = TRUE, widths = c(2, 2, 2), heights = c(2, 2)) +
plot_annotation(caption = "Seed mass (mg) [log scale]",
theme = theme(plot.caption = element_text(size = 35, family = "serif", hjust = 0.5, colour = "black")))
# Save figure
ggsave("figure/mass.seedling.fig.jpg", width = 18, height = 6, units = "in", dpi = 300)
mass.seedling.fig
3.5.5.3 Figure for pollutant
Above-ground biomass
# Data
ab.fig.pt.data <- filter(ab_result, Predictor %in% "Pollution")
ab.fig.pt.data$Predictortype <- factor(ab.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
ab.pt.fig <- ggplot(data = ab.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Emerging" = 16, "Inorganic" = 16, "Metal" = 16, "Organic" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#43CD80")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#43CD80")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
labs(title = "Above-ground biomass") +
theme(axis.text.y = element_text(size = 37, colour = "black", family = "serif"),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(ab.fig.pt.data$Effectsize_count[1])) + # Emerging
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(ab.fig.pt.data$Effectsize_count[4])) + # Organic
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pt.data$Effectsize_count[2])) + # Inorganic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(ab.fig.pt.data$Effectsize_count[3])) + # Metal
coord_flip() +
guides(fill = "none", color = "none")
ab.pt.fig
Below-ground biomass
# Data
be.fig.pt.data <- filter(be_result, Predictor %in% "Pollution")
be.fig.pt.data$Predictortype <- factor(be.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
be.pt.fig <- ggplot(data = be.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Emerging" = 16, "Inorganic" = 16, "Metal" = 16, "Organic" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.3), breaks = seq(-2, 1, 1)) +
labs(title = "Below-ground biomass") +
theme(axis.text.y = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(be.fig.pt.data$Effectsize_count[1])) + # Emerging
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(be.fig.pt.data$Effectsize_count[4])) + # Organic
annotate("text", x = 3, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pt.data$Effectsize_count[2])) + # Inorganic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(be.fig.pt.data$Effectsize_count[3])) + # Metal
coord_flip() +
guides(fill = "none", color = "none")
be.pt.fig
Seedling total biomass
# Data
sl.fig.pt.data <- filter(sl_result, Predictor %in% "Pollution")
sl.fig.pt.data$Predictortype <- factor(sl.fig.pt.data$Predictortype, levels = rev(Predictortype.1_level))
sl.pt.fig <- ggplot(data = sl.fig.pt.data, mapping = aes(x = Predictortype, y = Mean), fill = Predictortype) +
geom_point(size = 10, stroke = 3,
shape = c("Emerging" = 16, "Inorganic" = 16, "Metal" = 16, "Organic" = 16),
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#ABABAB")) +
geom_errorbar(aes(ymin = Lower_95_CI, ymax = Upper_95_CI),
linewidth = 4, alpha = 1, width = 0,
color = c("Emerging" = "#ABABAB", "Inorganic" = "#ABABAB", "Metal" = "#43CD80", "Organic" = "#ABABAB")) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
scale_y_continuous(limits = c(-2.1, 1.5), breaks = seq(-2, 1, 1)) +
labs(title = "Seedling biomass") +
theme(axis.text.y = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
annotate("text", x = 1, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pt.data$Effectsize_count[1])) + # Emerging
annotate("text", x = 2, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pt.data$Effectsize_count[4])) + # Organic
annotate("text", x = 3, y = -2, family = "serif", size = 14,
label = paste(sl.fig.pt.data$Effectsize_count[2])) + # Inorganic
annotate("text", x = 4, y = -1.95, family = "serif", size = 14,
label = paste(sl.fig.pt.data$Effectsize_count[3])) + # Metal
coord_flip() +
guides(fill = "none", color = "none")
sl.pt.fig
Combine pollution type figures
poll.seedling.fig <- ab.pt.fig + be.pt.fig + sl.pt.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2), heights = c(3, 2)) +
plot_annotation(caption = "Effect size (Log response ratio)",
theme = theme(plot.caption = element_text(size = 45, family = "serif",
hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/poll.seedling.fig.jpg", width = 20, height = 9, units = "in", dpi = 300)
poll.seedling.fig
Heavy metal concentration
Above-ground biomass
# Data
ab.con.data <- data.frame(AspectType = "AB",
predict(ab.qm.result.con, addx = TRUE))
ab.con.fig <- ggplot() +
geom_point(dat_hm_ab,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_line(data = ab.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
geom_ribbon(data = ab.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.2, 1.6), breaks = seq(-6, 1.5, 1.5),
name = "Log response ratio") +
theme_plot +
theme(axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_text(size = 35, colour = "black", family = "serif"),
axis.title.y = element_text(size = 35, colour = "black", family = "serif",
margin = margin(r = 20))) +
annotate("text", x = -1.3, y = -5, family = "serif", size = 13,
label = "slope = -0.226") + # slope value
annotate("text", x = -0.9, y = -6, family = "serif", size = 13,
label = "p < 0.001") # slope p
ab.con.fig
Below-ground biomass
# Data
be.con.data <- data.frame(AspectType = "BE",
predict(be.qm.result.con, addx = TRUE))
be.con.fig <- ggplot() +
geom_point(dat_hm_be,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_line(data = be.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
geom_ribbon(data = be.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.2, 1.6), breaks = seq(-6, 1.5, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
annotate("text", x = -1.3, y = -5, family = "serif", size = 13,
label = "slope = -0.527") + # slope value
annotate("text", x = -0.9, y = -6, family = "serif", size = 13,
label = "p < 0.001") # slope p
be.con.fig
Seedling total biomass
# Data
sl.con.data <- data.frame(AspectType = "SL",
predict(sl.qm.result.con, addx = TRUE))
sl.con.fig <- ggplot() +
geom_point(dat_hm_sl,
mapping = aes(x = Logcon, y = yi, size = 1/vi),
alpha = 0.3, shape = 1, color = "#43CD80", stroke = 1.2) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 1.5) +
geom_line(data = sl.con.data,
mapping = aes(x = X.Logcon, y = pred),
linetype = 1, linewidth = 1.5, colour = "black") +
geom_ribbon(data = sl.con.data,
aes(x = X.Logcon, y = pred, ymin = ci.lb, ymax = ci.ub),
fill = "#575757", alpha = 0.5) +
scale_size(range = c(3, 12)) +
scale_x_continuous(limits = c(-3.1, 4.1),
breaks = seq(-3, 4, 2),
labels = c("0.001", "0.1", "10", "1000")) +
scale_y_continuous(limits = c(-6.2, 1.6), breaks = seq(-6, 1.5, 1.5)) +
theme_plot +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 35, colour = "black", family = "serif"),
axis.text.y = element_blank()) +
annotate("text", x = -1.3, y = -5, family = "serif", size = 13,
label = "slope = -0.485") + #slope value
annotate("text", x = -0.9, y = -6, family = "serif", size = 13,
label = "p < 0.001") #slope p
sl.con.fig
Combine pollution concentration figures
hm.seedling.fig <- ab.con.fig + be.con.fig + sl.con.fig +
plot_layout(nrow= 1, byrow = TRUE, widths = c(2, 2, 2), heights = c(2, 2)) +
plot_annotation(caption = "Concentration of metals (mg/L) [log scale]",
theme = theme(plot.caption = element_text(size = 40, family = "serif",
hjust = 0.5, colour = "black"))) +
theme(plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"))
# Save figure
ggsave("figure/hm.seedling.fig.jpg", width = 20, height = 7, units = "in", dpi = 300)
hm.seedling.fig
3.6 Relative importance of predict variables
# Packages
library(performance)
library(glmulti)3.6.1 Germination percentage
# Check for co-linearity
vif.gp <- lm(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_gp)
gp.result <- check_collinearity(vif.gp)
gp.result## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## GrowthForm 2.26 [2.08, 2.47] 1.50 0.44 [0.41, 0.48]
## LifeSpan 1.37 [1.28, 1.47] 1.17 0.73 [0.68, 0.78]
## Dormancy 1.81 [1.67, 1.96] 1.34 0.55 [0.51, 0.60]
## SeedMass 1.09 [1.04, 1.18] 1.04 0.92 [0.85, 0.96]
## Pollution 1.94 [1.80, 2.12] 1.39 0.51 [0.47, 0.56]
## Concentration 1.08 [1.04, 1.18] 1.04 0.92 [0.85, 0.96]# Model selection
rma.glmulti <- function(formula, data, ...) {
rma.mv(formula, vi, random = list(~1|Ref_ID/Effectsize_ID, ~1|Species),
data = data, method = "REML", ...)}
gp.selection<- glmulti(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_gp, level = 1, fitfunction = rma.glmulti,
crit = "aicc", confsetsize = 64)## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...##
## After 50 models:
## Best model: yi~1+GrowthForm+LifeSpan+Dormancy+SeedMass+Concentration
## Crit= 2782.12143884147
## Mean crit= 2819.86856889052
## Completed.# Aicc for models
gp.modelweights <- weightable(gp.selection)
write_csv(gp.modelweights, "table/gp.modelweights.csv")
knitr::kable(gp.modelweights, digits = 3,
caption = "Summary of model weights of germination percentage",
font_size = 18, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| model | aicc | weights |
|---|---|---|
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 2775.176 | 0.257 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution + Concentration | 2775.374 | 0.233 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution + Concentration | 2776.314 | 0.145 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution + Concentration | 2776.500 | 0.133 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 2778.605 | 0.046 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution + Concentration | 2778.708 | 0.044 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution + Concentration | 2779.734 | 0.026 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution + Concentration | 2779.889 | 0.024 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution + Concentration | 2780.906 | 0.015 |
| yi ~ 1 + Dormancy + SeedMass + Pollution + Concentration | 2780.985 | 0.014 |
| yi ~ 1 + GrowthForm + Pollution + Concentration | 2781.052 | 0.014 |
| yi ~ 1 + Dormancy + Pollution + Concentration | 2781.052 | 0.014 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Concentration | 2782.121 | 0.008 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Concentration | 2782.446 | 0.007 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Concentration | 2783.371 | 0.004 |
| yi ~ 1 + GrowthForm + Dormancy + Concentration | 2783.683 | 0.004 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution + Concentration | 2784.276 | 0.003 |
| yi ~ 1 + LifeSpan + Pollution + Concentration | 2784.330 | 0.003 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Concentration | 2785.201 | 0.002 |
| yi ~ 1 + LifeSpan + Dormancy + Concentration | 2785.411 | 0.002 |
| yi ~ 1 + SeedMass + Pollution + Concentration | 2786.951 | 0.001 |
| yi ~ 1 + Pollution + Concentration | 2786.979 | 0.001 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Concentration | 2787.075 | 0.001 |
| yi ~ 1 + GrowthForm + LifeSpan + Concentration | 2787.378 | 0.001 |
| yi ~ 1 + Dormancy + SeedMass + Concentration | 2787.753 | 0.000 |
| yi ~ 1 + Dormancy + Concentration | 2787.915 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Concentration | 2788.376 | 0.000 |
| yi ~ 1 + GrowthForm + Concentration | 2788.669 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Concentration | 2790.935 | 0.000 |
| yi ~ 1 + LifeSpan + Concentration | 2791.104 | 0.000 |
| yi ~ 1 + SeedMass + Concentration | 2793.794 | 0.000 |
| yi ~ 1 + Concentration | 2793.921 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution | 2829.247 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution | 2829.488 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution | 2830.366 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution | 2830.609 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution | 2831.659 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution | 2831.893 | 0.000 |
| yi ~ 1 + Dormancy + Pollution | 2833.715 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution | 2833.891 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass + Pollution | 2833.934 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution | 2834.146 | 0.000 |
| yi ~ 1 + GrowthForm + Pollution | 2835.042 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution | 2835.299 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy | 2835.352 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass | 2835.582 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy | 2836.553 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass | 2836.787 | 0.000 |
| yi ~ 1 + LifeSpan + Pollution | 2837.347 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution | 2837.575 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy | 2837.760 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass | 2837.987 | 0.000 |
| yi ~ 1 + Pollution | 2839.673 | 0.000 |
| yi ~ 1 + SeedMass + Pollution | 2839.878 | 0.000 |
| yi ~ 1 + Dormancy | 2840.014 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass | 2840.221 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan | 2840.294 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass | 2840.539 | 0.000 |
| yi ~ 1 + GrowthForm | 2841.536 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass | 2841.783 | 0.000 |
| yi ~ 1 + LifeSpan | 2843.446 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass | 2843.675 | 0.000 |
| yi ~ 1 | 2845.969 | 0.000 |
| yi ~ 1 + SeedMass | 2846.172 | 0.000 |
# Extracting model results
setOldClass("rma.mv")
setMethod('getfit', 'rma.mv', function(object, ...) {
if (object$test == "z") {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = Inf)} else {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = object$k-object$p)}})
# Results
gp.selection.reslut <- round(coef(gp.selection), 3) %>% as.data.frame()
gp.selection.reslut <- gp.selection.reslut %>%
rownames_to_column(var = "predict variables")
write_csv(gp.selection.reslut, "table/gp.selection.reslut.csv")
knitr::kable(gp.selection.reslut, digits = 3,
caption = "Summary of relative importance of predict variables in germination percentage",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| predict variables | Estimate | Uncond. variance | Nb models | Importance | +/- (alpha=0.05) |
|---|---|---|---|---|---|
| SeedMass | 0.000 | 0.000 | 32 | 0.523 | 0.000 |
| LifeSpanShort-lived | -0.042 | 0.013 | 32 | 0.655 | 0.225 |
| GrowthFormWoody | 0.199 | 0.020 | 32 | 0.871 | 0.276 |
| DormancyMPD | -0.277 | 0.300 | 32 | 0.912 | 1.073 |
| DormancyND | -0.261 | 0.212 | 32 | 0.912 | 0.903 |
| DormancyPD | -0.291 | 0.207 | 32 | 0.912 | 0.892 |
| DormancyPY | -0.236 | 0.220 | 32 | 0.912 | 0.919 |
| PollutionInorganic | -0.686 | 0.186 | 32 | 0.971 | 0.845 |
| PollutionMetal | -0.503 | 0.157 | 32 | 0.971 | 0.775 |
| PollutionOrganic | -0.626 | 0.172 | 32 | 0.971 | 0.812 |
| Concentration | 0.000 | 0.000 | 32 | 1.000 | 0.000 |
| intrcpt | 0.339 | 0.359 | 64 | 1.000 | 1.175 |
Figure
# Data for figure
gp.selection.fig.data <- data.frame(Aspect = "GP",
Predictor = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
Importance = c(gp.selection.reslut$Importance[1],
gp.selection.reslut$Importance[2],
gp.selection.reslut$Importance[3],
gp.selection.reslut$Importance[4],
gp.selection.reslut$Importance[8],
gp.selection.reslut$Importance[11]))
# Figure
gp.selection.fig <- ggplot(data = gp.selection.fig.data,
mapping = aes(x = Predictor, y = Importance, fill = Predictor)) +
geom_bar(stat = "identity",
width = 0.6, colour = "black", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("Concentration" = "#F39B7FCC",
"Dormancy" = "#4DBBD5CC",
"Growth form" = "#4DBBD5CC",
"Life span" = "#4DBBD5CC",
"Pollution" = "#F39B7FCC",
"Seed mass" = "#4DBBD5CC")) +
scale_x_discrete(
limits = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
breaks = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
labels = c("Seed mass", "Adult lifespan", "Growth form",
"Dormancy class", "Pollutant type", "Pollutant concentration")) +
scale_y_continuous(name = "Relative importance", expand = c(0, 0),
limits = c(0, 1.05), breaks = seq(0, 1, 0.2)) +
labs(title = "Germination percentage") +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_text(size = 37, color = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
guides(fill = "none") +
coord_flip()
gp.selection.fig
3.6.2 Germination speed
# Check for co-linearity
vif.gs <- lm(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_gs)
gs.result <- check_collinearity(vif.gs)
gs.result## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## GrowthForm 1.92 [1.66, 2.29] 1.39 0.52 [0.44, 0.60]
## LifeSpan 1.15 [1.06, 1.39] 1.07 0.87 [0.72, 0.94]
## Dormancy 2.19 [1.87, 2.62] 1.48 0.46 [0.38, 0.53]
## SeedMass 1.19 [1.09, 1.42] 1.09 0.84 [0.70, 0.92]
## Pollution 1.99 [1.71, 2.38] 1.41 0.50 [0.42, 0.58]
## Concentration 1.05 [1.01, 1.51] 1.03 0.95 [0.66, 0.99]# Model selection
rma.glmulti <- function(formula, data, ...) {
rma.mv(formula, vi, random = list(~1|Ref_ID/Effectsize_ID, ~1|Species),
data = data, method = "REML", ...)}
gs.selection<- glmulti(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_gs, level = 1, fitfunction = rma.glmulti,
crit = "aicc", confsetsize = 64)## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...##
## After 50 models:
## Best model: yi~1+GrowthForm+LifeSpan+Dormancy+SeedMass+Concentration
## Crit= 386.405900765363
## Mean crit= 428.517397048951
## Completed.# Aicc for models
gs.modelweights <- weightable(gs.selection)
write_csv(gs.modelweights, "table/gs.modelweights.csv")
knitr::kable(gs.modelweights, digits = 3,
caption = "Summary of model weights of germination speed",
font_size = 18, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| model | aicc | weights |
|---|---|---|
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 384.828 | 0.174 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 385.529 | 0.123 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Concentration | 386.406 | 0.079 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution + Concentration | 386.437 | 0.078 |
| yi ~ 1 + Dormancy + SeedMass + Pollution + Concentration | 387.138 | 0.055 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Concentration | 387.190 | 0.053 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution + Concentration | 387.301 | 0.051 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution + Concentration | 387.488 | 0.046 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Concentration | 388.030 | 0.035 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution + Concentration | 388.292 | 0.031 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution + Concentration | 388.309 | 0.031 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Concentration | 388.462 | 0.028 |
| yi ~ 1 + Dormancy + SeedMass + Concentration | 388.911 | 0.023 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution + Concentration | 388.927 | 0.022 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution + Concentration | 389.125 | 0.020 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Concentration | 389.206 | 0.019 |
| yi ~ 1 + LifeSpan + SeedMass + Concentration | 389.567 | 0.016 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution + Concentration | 389.896 | 0.014 |
| yi ~ 1 + Dormancy + Pollution + Concentration | 389.901 | 0.014 |
| yi ~ 1 + LifeSpan + Dormancy + Concentration | 390.070 | 0.013 |
| yi ~ 1 + GrowthForm + SeedMass + Concentration | 390.105 | 0.012 |
| yi ~ 1 + SeedMass + Pollution + Concentration | 390.275 | 0.011 |
| yi ~ 1 + GrowthForm + Dormancy + Concentration | 390.855 | 0.009 |
| yi ~ 1 + LifeSpan + Pollution + Concentration | 390.898 | 0.008 |
| yi ~ 1 + GrowthForm + LifeSpan + Concentration | 391.391 | 0.007 |
| yi ~ 1 + GrowthForm + Pollution + Concentration | 391.548 | 0.006 |
| yi ~ 1 + SeedMass + Concentration | 391.681 | 0.006 |
| yi ~ 1 + Dormancy + Concentration | 391.755 | 0.005 |
| yi ~ 1 + LifeSpan + Concentration | 392.403 | 0.004 |
| yi ~ 1 + Pollution + Concentration | 392.788 | 0.003 |
| yi ~ 1 + GrowthForm + Concentration | 393.059 | 0.003 |
| yi ~ 1 + Concentration | 394.348 | 0.001 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution | 442.769 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution | 443.782 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution | 444.623 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution | 444.724 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass + Pollution | 445.700 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution | 445.768 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass | 446.213 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution | 446.626 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution | 446.676 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass | 447.312 | 0.000 |
| yi ~ 1 + Dormancy + Pollution | 447.698 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution | 447.900 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy | 448.041 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass | 448.175 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution | 448.682 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass | 448.855 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution | 449.100 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy | 449.175 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass | 449.345 | 0.000 |
| yi ~ 1 + SeedMass + Pollution | 449.941 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy | 450.032 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass | 450.136 | 0.000 |
| yi ~ 1 + LifeSpan + Pollution | 450.421 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass | 450.960 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan | 451.079 | 0.000 |
| yi ~ 1 + GrowthForm + Pollution | 451.160 | 0.000 |
| yi ~ 1 + Dormancy | 451.208 | 0.000 |
| yi ~ 1 + SeedMass | 452.356 | 0.000 |
| yi ~ 1 + LifeSpan | 452.442 | 0.000 |
| yi ~ 1 + Pollution | 452.442 | 0.000 |
| yi ~ 1 + GrowthForm | 453.207 | 0.000 |
| yi ~ 1 | 454.575 | 0.000 |
# Extracting model results
setOldClass("rma.mv")
setMethod('getfit', 'rma.mv', function(object, ...) {
if (object$test == "z") {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = Inf)} else {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = object$k-object$p)}})
# Results
gs.selection.reslut <- round(coef(gs.selection), 3) %>% as.data.frame()
gs.selection.reslut <- gs.selection.reslut %>%
rownames_to_column(var = "predict variables")
write_csv(gs.selection.reslut, "table/gs.selection.reslut.csv")
knitr::kable(gs.selection.reslut, digits = 3,
caption = "Summary of relative importance of predict variables in germination speed",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| predict variables | Estimate | Uncond. variance | Nb models | Importance | +/- (alpha=0.05) |
|---|---|---|---|---|---|
| GrowthFormWoody | 0.017 | 0.037 | 32 | 0.603 | 0.378 |
| PollutionMetal | -0.094 | 0.044 | 32 | 0.686 | 0.410 |
| PollutionOrganic | -0.015 | 0.076 | 32 | 0.686 | 0.540 |
| LifeSpanShort-lived | -0.343 | 0.271 | 32 | 0.696 | 1.020 |
| DormancyPD | -0.019 | 0.034 | 32 | 0.776 | 0.362 |
| DormancyPY | 0.246 | 0.092 | 32 | 0.776 | 0.593 |
| SeedMass | 0.000 | 0.000 | 32 | 0.797 | 0.000 |
| Concentration | 0.000 | 0.000 | 32 | 1.000 | 0.000 |
| intrcpt | -0.247 | 0.078 | 64 | 1.000 | 0.548 |
Figure
# Data for figure
gs.selection.fig.data <- data.frame(Aspect = "GS",
Predictor = c("Growth form", "Pollution", "Life span",
"Dormancy", "Seed mass", "Concentration"),
Importance = c(gs.selection.reslut$Importance[1],
gs.selection.reslut$Importance[2],
gs.selection.reslut$Importance[4],
gs.selection.reslut$Importance[5],
gs.selection.reslut$Importance[7],
gs.selection.reslut$Importance[8]))
# Figure
gs.selection.fig <- ggplot(data = gs.selection.fig.data,
mapping = aes(x = Predictor, y = Importance, fill = Predictor)) +
geom_bar(stat = "identity",
width = 0.6, colour = "black", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("Concentration" = "#F39B7FCC",
"Dormancy" = "#4DBBD5CC",
"Growth form" = "#4DBBD5CC",
"Life span" = "#4DBBD5CC",
"Pollution" = "#F39B7FCC",
"Seed mass" = "#4DBBD5CC")) +
scale_x_discrete(limits = c("Growth form", "Pollution", "Life span",
"Dormancy", "Seed mass", "Concentration"),
breaks = c("Growth form", "Pollution", "Life span",
"Dormancy", "Seed mass", "Concentration"),
labels = c("Growth form", "Pollutant type", "Adult lifespan",
"Dormancy class", "Seed mass", "Pollutant concentration")) +
scale_y_continuous(name = "Relative importance", expand = c(0, 0),
limits = c(0, 1.05), breaks = seq(0, 1, 0.2)) +
labs(title = "Germination speed") +
theme(axis.title.x.bottom = element_blank(),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_text(size = 37, color = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
guides(fill = "none") +
coord_flip()
gs.selection.fig
3.6.3 Seedling growth
# Check for co-linearity
vif.sg <- lm(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_sg)
sg.result <- check_collinearity(vif.sg)
sg.result## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## GrowthForm 2.19 [2.05, 2.34] 1.48 0.46 [0.43, 0.49]
## LifeSpan 1.56 [1.48, 1.66] 1.25 0.64 [0.60, 0.68]
## Dormancy 1.90 [1.79, 2.02] 1.38 0.53 [0.49, 0.56]
## SeedMass 1.09 [1.05, 1.15] 1.04 0.92 [0.87, 0.95]
## Pollution 1.57 [1.49, 1.67] 1.25 0.64 [0.60, 0.67]
## Concentration 1.04 [1.01, 1.13] 1.02 0.96 [0.89, 0.99]# Model selection
rma.glmulti <- function(formula, data, ...) {
rma.mv(formula, vi, random = list(~1|Ref_ID/Effectsize_ID, ~1|Species),
data = data, method = "REML", ...)}
sg.selection<- glmulti(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_sg, level = 1, fitfunction = rma.glmulti,
crit = "aicc", confsetsize = 64)## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...##
## After 50 models:
## Best model: yi~1+GrowthForm+LifeSpan+Dormancy+SeedMass+Concentration
## Crit= 4648.40342326411
## Mean crit= 4836.26614180568
## Completed.# Aicc for models
sg.modelweights <- weightable(sg.selection)
write_csv(sg.modelweights, "table/sg.modelweights.csv")
knitr::kable(sg.modelweights, digits = 3,
caption = "Summary of model weights of seedling growth",
font_size = 18, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| model | aicc | weights |
|---|---|---|
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 4633.040 | 0.400 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution + Concentration | 4633.127 | 0.383 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution + Concentration | 4636.096 | 0.087 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution + Concentration | 4636.200 | 0.082 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution + Concentration | 4640.045 | 0.012 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 4640.109 | 0.012 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution + Concentration | 4641.450 | 0.006 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution + Concentration | 4641.533 | 0.006 |
| yi ~ 1 + Dormancy + Pollution + Concentration | 4641.896 | 0.005 |
| yi ~ 1 + Dormancy + SeedMass + Pollution + Concentration | 4641.935 | 0.005 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution + Concentration | 4645.203 | 0.001 |
| yi ~ 1 + GrowthForm + Pollution + Concentration | 4645.300 | 0.001 |
| yi ~ 1 + LifeSpan + Pollution + Concentration | 4648.285 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution + Concentration | 4648.337 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Concentration | 4648.403 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Concentration | 4648.442 | 0.000 |
| yi ~ 1 + Pollution + Concentration | 4650.339 | 0.000 |
| yi ~ 1 + SeedMass + Pollution + Concentration | 4650.352 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Concentration | 4651.387 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + Concentration | 4651.437 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + Concentration | 4655.939 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Concentration | 4656.025 | 0.000 |
| yi ~ 1 + Dormancy + Concentration | 4657.829 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass + Concentration | 4657.892 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Concentration | 4658.553 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Concentration | 4658.576 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Concentration | 4662.213 | 0.000 |
| yi ~ 1 + GrowthForm + Concentration | 4662.241 | 0.000 |
| yi ~ 1 + LifeSpan + Concentration | 4666.287 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Concentration | 4666.371 | 0.000 |
| yi ~ 1 + Concentration | 4668.315 | 0.000 |
| yi ~ 1 + SeedMass + Concentration | 4668.367 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution | 4913.586 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution | 4914.236 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution | 4916.551 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution | 4917.243 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution | 4919.785 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution | 4920.262 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass + Pollution | 4921.868 | 0.000 |
| yi ~ 1 + Dormancy + Pollution | 4922.369 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution | 4926.000 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass | 4926.012 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution | 4926.500 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy | 4926.517 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass | 4928.714 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy | 4929.244 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution | 4929.925 | 0.000 |
| yi ~ 1 + GrowthForm + Pollution | 4930.460 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution | 4931.027 | 0.000 |
| yi ~ 1 + LifeSpan + Pollution | 4931.467 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass | 4932.454 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy | 4932.867 | 0.000 |
| yi ~ 1 + SeedMass + Pollution | 4933.802 | 0.000 |
| yi ~ 1 + Pollution | 4934.268 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass | 4934.477 | 0.000 |
| yi ~ 1 + Dormancy | 4934.897 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass | 4939.693 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan | 4940.100 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass | 4943.238 | 0.000 |
| yi ~ 1 + GrowthForm | 4943.664 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass | 4945.314 | 0.000 |
| yi ~ 1 + LifeSpan | 4945.768 | 0.000 |
| yi ~ 1 + SeedMass | 4947.832 | 0.000 |
| yi ~ 1 | 4948.283 | 0.000 |
# Extracting model results
setOldClass("rma.mv")
setMethod('getfit', 'rma.mv', function(object, ...) {
if (object$test == "z") {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = Inf)} else {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = object$k-object$p)}})
# Results
sg.selection.reslut <- round(coef(sg.selection), 3) %>% as.data.frame()
sg.selection.reslut <- sg.selection.reslut %>%
rownames_to_column(var = "predict variables")
write_csv(sg.selection.reslut, "table/sg.selection.reslut.csv")
knitr::kable(sg.selection.reslut, digits = 3,
caption = "Summary of relative importance of predict variables in seedling growth",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| predict variables | Estimate | Uncond. variance | Nb models | Importance | +/- (alpha=0.05) |
|---|---|---|---|---|---|
| SeedMass | 0.000 | 0.000 | 32 | 0.511 | 0.000 |
| LifeSpanShort-lived | 0.135 | 0.015 | 32 | 0.819 | 0.242 |
| GrowthFormWoody | 0.305 | 0.018 | 32 | 0.966 | 0.264 |
| DormancyMPD | -0.026 | 0.087 | 32 | 0.986 | 0.578 |
| DormancyND | 0.229 | 0.090 | 32 | 0.986 | 0.589 |
| DormancyPD | 0.073 | 0.080 | 32 | 0.986 | 0.554 |
| DormancyPY | -0.022 | 0.097 | 32 | 0.986 | 0.609 |
| PollutionInorganic | -0.444 | 0.095 | 32 | 1.000 | 0.605 |
| PollutionMetal | -0.642 | 0.052 | 32 | 1.000 | 0.447 |
| PollutionOrganic | -0.204 | 0.090 | 32 | 1.000 | 0.588 |
| intrcpt | -0.211 | 0.135 | 64 | 1.000 | 0.720 |
| Concentration | 0.000 | 0.000 | 32 | 1.000 | 0.000 |
Figure
# Data for figure
sg.selection.fig.data <- data.frame(Aspect = "SG",
Predictor = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
Importance = c(sg.selection.reslut$Importance[1],
sg.selection.reslut$Importance[2],
sg.selection.reslut$Importance[3],
sg.selection.reslut$Importance[4],
sg.selection.reslut$Importance[8],
sg.selection.reslut$Importance[11]))
# Figure
sg.selection.fig <- ggplot(data = sg.selection.fig.data,
mapping = aes(x = Predictor, y = Importance, fill = Predictor)) +
geom_bar(stat = "identity",
width = 0.6, colour = "black", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("Concentration" = "#F39B7FCC",
"Dormancy" = "#4DBBD5CC",
"Growth form" = "#4DBBD5CC",
"Life span" = "#4DBBD5CC",
"Pollution" = "#F39B7FCC",
"Seed mass" = "#4DBBD5CC")) +
scale_x_discrete(limits = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
breaks = c("Seed mass", "Life span", "Growth form",
"Dormancy", "Pollution", "Concentration"),
labels = c("Seed mass", "Adult lifespan", "Growth form",
"Dormancy class", "Pollutant type", "Pollutant concentration"
)) +
scale_y_continuous(name = "Relative importance", expand = c(0, 0),
limits = c(0, 1.05), breaks = seq(0, 1, 0.2)) +
labs(title = "Seedling growth") +
theme(axis.title.x.bottom = element_text(size = 37, color = "black",
family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_text(size = 37, color = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
guides(fill = "none") +
coord_flip()
sg.selection.fig
3.6.4 Biomass allocation
# Check for co-linearity
vif.ba <- lm(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_ba)
ba.result <- check_collinearity(vif.ba)
ba.result## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## GrowthForm 2.03 [1.84, 2.27] 1.43 0.49 [0.44, 0.54]
## LifeSpan 1.59 [1.46, 1.77] 1.26 0.63 [0.57, 0.69]
## Dormancy 1.71 [1.56, 1.91] 1.31 0.58 [0.52, 0.64]
## SeedMass 1.14 [1.07, 1.26] 1.07 0.88 [0.79, 0.93]
## Pollution 1.51 [1.39, 1.68] 1.23 0.66 [0.60, 0.72]
## Concentration 1.04 [1.01, 1.27] 1.02 0.96 [0.79, 0.99]# Model selection
rma.glmulti <- function(formula, data, ...) {
rma.mv(formula, vi, random = list(~1|Ref_ID/Effectsize_ID, ~1|Species),
data = data, method = "REML", ...)}
ba.selection<- glmulti(yi ~ GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration,
data = dat_ba, level = 1, fitfunction = rma.glmulti,
crit = "aicc", confsetsize = 64)## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...##
## After 50 models:
## Best model: yi~1+GrowthForm+LifeSpan+Dormancy+SeedMass+Concentration
## Crit= 1510.10442452002
## Mean crit= 1542.51547622363
## Completed.# Aicc for models
ba.modelweights <- weightable(ba.selection)
write_csv(ba.modelweights, "table/ba.modelweights.csv")
knitr::kable(ba.modelweights, digits = 3,
caption = "Summary of model weights of biomass allocation",
font_size = 18, full_width = T, align = "c",
"html") %>%
kable_styling(position = "center")| model | aicc | weights |
|---|---|---|
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 1504.337 | 0.353 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution + Concentration | 1504.877 | 0.269 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution + Concentration | 1506.358 | 0.128 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution + Concentration | 1506.909 | 0.097 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution + Concentration | 1508.647 | 0.041 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution + Concentration | 1509.415 | 0.028 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Concentration | 1510.104 | 0.020 |
| yi ~ 1 + Dormancy + SeedMass + Pollution + Concentration | 1510.409 | 0.017 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Concentration | 1510.630 | 0.015 |
| yi ~ 1 + Dormancy + Pollution + Concentration | 1511.194 | 0.011 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Concentration | 1512.218 | 0.007 |
| yi ~ 1 + GrowthForm + Dormancy + Concentration | 1512.747 | 0.005 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Concentration | 1514.982 | 0.002 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution + Concentration | 1515.160 | 0.002 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution + Concentration | 1515.744 | 0.001 |
| yi ~ 1 + LifeSpan + Dormancy + Concentration | 1515.847 | 0.001 |
| yi ~ 1 + Dormancy + SeedMass + Concentration | 1516.784 | 0.001 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution + Concentration | 1517.176 | 0.001 |
| yi ~ 1 + Dormancy + Concentration | 1517.678 | 0.000 |
| yi ~ 1 + GrowthForm + Pollution + Concentration | 1517.765 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution + Concentration | 1519.495 | 0.000 |
| yi ~ 1 + LifeSpan + Pollution + Concentration | 1520.320 | 0.000 |
| yi ~ 1 + SeedMass + Pollution + Concentration | 1521.367 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Concentration | 1521.473 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Concentration | 1522.103 | 0.000 |
| yi ~ 1 + Pollution + Concentration | 1522.224 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Concentration | 1523.553 | 0.000 |
| yi ~ 1 + GrowthForm + Concentration | 1524.178 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Concentration | 1526.739 | 0.000 |
| yi ~ 1 + LifeSpan + Concentration | 1527.692 | 0.000 |
| yi ~ 1 + SeedMass + Concentration | 1528.700 | 0.000 |
| yi ~ 1 + Concentration | 1529.706 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass + Pollution | 1542.248 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + Pollution | 1542.858 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass + Pollution | 1543.916 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + Pollution | 1544.532 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass + Pollution | 1544.580 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + Pollution | 1545.346 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass + Pollution | 1546.361 | 0.000 |
| yi ~ 1 + Dormancy + Pollution | 1547.146 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy + SeedMass | 1550.318 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Dormancy | 1550.915 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy + SeedMass | 1552.076 | 0.000 |
| yi ~ 1 + GrowthForm + Dormancy | 1552.673 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass + Pollution | 1552.916 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy + SeedMass | 1553.125 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + Pollution | 1553.556 | 0.000 |
| yi ~ 1 + LifeSpan + Dormancy | 1553.937 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass + Pollution | 1554.633 | 0.000 |
| yi ~ 1 + Dormancy + SeedMass | 1554.973 | 0.000 |
| yi ~ 1 + GrowthForm + Pollution | 1555.274 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass + Pollution | 1555.749 | 0.000 |
| yi ~ 1 + Dormancy | 1555.815 | 0.000 |
| yi ~ 1 + LifeSpan + Pollution | 1556.570 | 0.000 |
| yi ~ 1 + SeedMass + Pollution | 1557.817 | 0.000 |
| yi ~ 1 + Pollution | 1558.681 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan + SeedMass | 1561.843 | 0.000 |
| yi ~ 1 + GrowthForm + LifeSpan | 1562.472 | 0.000 |
| yi ~ 1 + GrowthForm + SeedMass | 1563.611 | 0.000 |
| yi ~ 1 + GrowthForm | 1564.240 | 0.000 |
| yi ~ 1 + LifeSpan + SeedMass | 1565.391 | 0.000 |
| yi ~ 1 + LifeSpan | 1566.253 | 0.000 |
| yi ~ 1 + SeedMass | 1567.575 | 0.000 |
| yi ~ 1 | 1568.499 | 0.000 |
# Extracting model results
setOldClass("rma.mv")
setMethod('getfit', 'rma.mv', function(object, ...) {
if (object$test == "z") {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = Inf)} else {
cbind(estimate = coef(object), se = sqrt(diag(vcov(object))), df = object$k-object$p)}})
# Results
ba.selection.reslut <- round(coef(ba.selection), 3) %>% as.data.frame()
ba.selection.reslut <- ba.selection.reslut %>%
rownames_to_column(var = "predict variables")
write_csv(ba.selection.reslut, "table/ba.selection.reslut.csv")
knitr::kable(ba.selection.reslut, digits = 3,
caption = "Summary of relative importance of predict variables in biomass allocation",
font_size = 20, align = "c", full_width = T,
"html") %>%
kable_styling(position = "center")| predict variables | Estimate | Uncond. variance | Nb models | Importance | +/- (alpha=0.05) |
|---|---|---|---|---|---|
| SeedMass | 0.000 | 0.000 | 32 | 0.570 | 0.000 |
| LifeSpanShort-lived | 0.113 | 0.053 | 32 | 0.731 | 0.451 |
| GrowthFormWoody | 0.383 | 0.072 | 32 | 0.899 | 0.525 |
| PollutionInorganic | -0.394 | 0.279 | 32 | 0.949 | 1.035 |
| PollutionMetal | -0.447 | 0.176 | 32 | 0.949 | 0.823 |
| PollutionOrganic | -0.203 | 0.275 | 32 | 0.949 | 1.028 |
| DormancyMPD | -0.176 | 0.631 | 32 | 0.996 | 1.557 |
| DormancyND | 0.054 | 0.477 | 32 | 0.996 | 1.353 |
| DormancyPD | -0.356 | 0.453 | 32 | 0.996 | 1.318 |
| DormancyPY | -0.323 | 0.520 | 32 | 0.996 | 1.413 |
| Concentration | 0.000 | 0.000 | 32 | 1.000 | 0.000 |
| intrcpt | -0.002 | 0.613 | 64 | 1.000 | 1.535 |
Figure
# Data for figure
ba.selection.fig.data <- data.frame(Aspect = "BA",
Predictor = c("Seed mass", "Life span", "Growth form",
"Pollution" , "Dormancy", "Concentration"),
Importance = c(ba.selection.reslut$Importance[1],
ba.selection.reslut$Importance[2],
ba.selection.reslut$Importance[3],
ba.selection.reslut$Importance[4],
ba.selection.reslut$Importance[7],
ba.selection.reslut$Importance[11]))
# Figure
ba.selection.fig <- ggplot(data = ba.selection.fig.data,
mapping = aes(x = Predictor, y = Importance, fill = Predictor)) +
geom_bar(stat = "identity",
width = 0.6, colour = "black", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = 2, linewidth = 1.2) +
scale_fill_manual(values = c("Concentration" = "#F39B7FCC",
"Dormancy" = "#4DBBD5CC",
"Growth form" = "#4DBBD5CC",
"Life span" = "#4DBBD5CC",
"Pollution" = "#F39B7FCC",
"Seed mass" = "#4DBBD5CC")) +
scale_x_discrete(limits = c("Seed mass", "Life span", "Growth form",
"Pollution" , "Dormancy", "Concentration"),
breaks = c("Seed mass", "Life span", "Growth form",
"Pollution" , "Dormancy", "Concentration"),
labels = c("Seed mass", "Adult lifespan", "Growth form",
"Pollutant type", "Dormancy class", "Pollutant concentration")) +
scale_y_continuous(name = "Relative importance", expand = c(0, 0),
limits = c(0, 1.05), breaks = seq(0, 1, 0.2)) +
labs(title = "Biomass allocation") +
theme(axis.title.x.bottom = element_text(size = 37, color = "black",
family = "serif", hjust = 0.5),
axis.text.x = element_text(size = 37, colour = "black", family = "serif"),
axis.text.y = element_text(size = 37, color = "black", family = "serif"),
plot.title = element_text(size = 37, face = "bold",
colour = "black", family = "serif", hjust = 0.5)) +
theme_plot +
guides(fill = "none") +
coord_flip()
ba.selection.fig
Combine model selection figures
selection.fig <- gp.selection.fig + sg.selection.fig + gs.selection.fig + ba.selection.fig +
plot_layout(nrow = 2, byrow = 2, widths = c(2, 2), heights = c(2, 2)) +
theme(plot.margin = unit(c(0.35, 0.35, 0.35, 0.35), "cm"))
# Save figure
ggsave("figure/selection.bind.jpg", width = 28, height = 16, units = "in", dpi = 300)
selection.fig
4 Publication bias
4.1 Rosenberg fail-safe numbers
Germination percentage
fsn(yi, vi, data = dat_gp)##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 55278734Germination speed
fsn(yi, vi, data = dat_gs)##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 1315254Seedling growth
fsn(yi, vi, data = dat_sg)##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 90153187Biomass allocation
fsn(yi, vi, data = dat_ba)##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 394603404.2 Egger’s test and funnel plots
Germination percentage
gp.bia <- rma(yi, vi, data = dat_gp)
regtest(gp.bia)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = -9.7706, p < .0001
## Limit Estimate (as sei -> 0): b = -0.1360 (CI: -0.1906, -0.0814)trim.gp <- trimfill(gp.bia, side = "l")
trim.gp##
## Estimated number of missing studies on the left side: 388 (SE = 22.7476)
##
## Random-Effects Model (k = 1660; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.6143 (SE = 0.0217)
## tau (square root of estimated tau^2 value): 0.7838
## I^2 (total heterogeneity / total variability): 99.98%
## H^2 (total variability / sampling variability): 5953.12
##
## Test for Heterogeneity:
## Q(df = 1659) = 3694734.1707, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5049 0.0194 -26.0216 <.0001 -0.5429 -0.4669 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#jpeg("figure/trim.gp.fig.jpg", width = 1800, height = 1400, res = 300)
trim.gp.fig <- funnel(trim.gp, digits = 3, cex.lab = 1.5, family = "serif",
xlab = "Germination percentage observed outcome",
ylab = "Standard error")
#dev.off()Germination speed
gs.bia <- rma(yi, vi, data = dat_gs)
regtest(gs.bia)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = -6.1772, p < .0001
## Limit Estimate (as sei -> 0): b = -0.1379 (CI: -0.2282, -0.0475)trim.gs <- trimfill(gs.bia, side = "l")
trim.gs##
## Estimated number of missing studies on the left side: 52 (SE = 11.4258)
##
## Random-Effects Model (k = 358; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3969 (SE = 0.0307)
## tau (square root of estimated tau^2 value): 0.6300
## I^2 (total heterogeneity / total variability): 99.79%
## H^2 (total variability / sampling variability): 468.58
##
## Test for Heterogeneity:
## Q(df = 357) = 99815.7988, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.4763 0.0339 -14.0533 <.0001 -0.5428 -0.4099 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#jpeg("figure/trim.gs.fig.jpg", width = 1800, height = 1400, res = 300)
trim.gs.fig <- funnel(trim.gs, digits = 3, cex.lab = 1.5, family = "serif",
xlab = "Germination speed observed outcome",
ylab = "Standard error")
#dev.off()Seedling growth
sg.bia <- rma(yi, vi, data = dat_sg)
regtest(sg.bia)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = -3.8670, p = 0.0001
## Limit Estimate (as sei -> 0): b = -0.5384 (CI: -0.5953, -0.4815)trim.sg <- trimfill(sg.bia, side = "l")
trim.sg##
## Estimated number of missing studies on the left side: 418 (SE = 30.2169)
##
## Random-Effects Model (k = 2557; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.9771 (SE = 0.0281)
## tau (square root of estimated tau^2 value): 0.9885
## I^2 (total heterogeneity / total variability): 99.91%
## H^2 (total variability / sampling variability): 1087.48
##
## Test for Heterogeneity:
## Q(df = 2556) = 1941727.4018, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.8194 0.0198 -41.2860 <.0001 -0.8583 -0.7805 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#jpeg("figure/trim.sg.fig.jpg", width = 1800, height = 1400, res = 300)
trim.sg.fig <- funnel(trim.sg, digits = 3, cex.lab = 1.5, family = "serif",
xlab = "Seedling growth observed outcome",
ylab = "Standard error")
#dev.off()Biomass allocation
ba.bia <- rma(yi, vi, data = dat_ba)
regtest(ba.bia)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = 4.7594, p < .0001
## Limit Estimate (as sei -> 0): b = -0.6664 (CI: -0.7684, -0.5643)trim.ba <- trimfill(ba.bia)
trim.ba##
## Estimated number of missing studies on the left side: 146 (SE = 18.0075)
##
## Random-Effects Model (k = 905; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.9974 (SE = 0.0485)
## tau (square root of estimated tau^2 value): 0.9987
## I^2 (total heterogeneity / total variability): 99.96%
## H^2 (total variability / sampling variability): 2507.20
##
## Test for Heterogeneity:
## Q(df = 904) = 1777119.6938, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7052 0.0338 -20.8636 <.0001 -0.7715 -0.6390 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#jpeg("figure/trim.ba.fig.jpg", width = 1800, height = 1400, res = 300)
trim.ba.fig <- funnel(trim.ba, digits = 3, cex.lab = 1.5, family = "serif",
xlab = "Biomass allocation observed outcome",
ylab = "Standard error")
#dev.off()