Explore P-fertilisation experiments in MESI
Evan Perkowski
2024-01-17
Script for exploring P fertilisation and N*P fertilisation experiments included in the MESI database as of Jan. 16, 2024. Initial exploration shows number of experiments for each experiment type and types of measurements included in the database. Then, script conducts a meta-analysis as an initial assessment of P and N*P fertilization effects on leaf and whole plant traits
# Libraries
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(readr)
library(metafor)
library(MAd)
# MESI data
df_mesi <- read_csv("../../../../../../eaperkowski/git/mesi-db/data/mesi_main.csv")Explore data availability in MESI dataset for N-fertilisation experiments
explore_mesi_nfert <- df_mesi %>%
# fertilisation experiments only
filter(treatment == "f") %>%
# N-fertilisation only (without P or K addition)
filter(npk == "_100")
head(explore_mesi_nfert)## # A tibble: 6 × 60
## db id duplicat…¹ citat…² respo…³ site study exp lat lon eleva…⁴
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 sichuan s1757 s1757 barger… agb 10.5… 10.5… 10.5… 10.5 -66.9 NA
## 2 hebei h962 h962 stape_… agb -11.… -11.… -11.… -12.0 -38.1 250
## 3 hebei h963 h963 stape_… agb -11.… -11.… -11.… -12.0 -38.1 250
## 4 hebei h964 h964 stape_… agb -11.… -11.… -11.… -12.0 -38.1 250
## 5 hebei h965 h965 stape_… agb -11.… -11.… -11.… -12.0 -38.1 250
## 6 sichuan s170 s170 ohallo… agb -15.… -15.… -15.… -15.4 23.2 NA
## # … with 49 more variables: mat <dbl>, map <dbl>, ecosystem_type <chr>,
## # vegetation_type <chr>, experiment_type <chr>, community_type <chr>,
## # dominant_species <chr>, growth_form <chr>, age <dbl>,
## # disturbance_type <chr>, treatment <chr>, npk <chr>, w_t1 <chr>, c_c <dbl>,
## # c_t <dbl>, d_t <dbl>, d_t2 <dbl>, n_c <dbl>, n_t <dbl>, p_c <dbl>,
## # p_t <dbl>, k_c <dbl>, k_t <dbl>, i_c <dbl>, i_t <dbl>, i_t2 <dbl>,
## # s_c <dbl>, s_t <dbl>, w_t2 <dbl>, w_t3 <dbl>, start_year <dbl>, …## How many experiments?
length(unique(explore_mesi_nfert$citation))## [1] 606## What traits are available?
unique(explore_mesi_nfert$response)## [1] "agb"
## [2] "agb_c"
## [3] "agb_group"
## [4] "agb_height"
## [5] "agb_n"
## [6] "agb_n_stock"
## [7] "agb_ndvi"
## [8] "agb_p"
## [9] "agb_panicle"
## [10] "anpp"
## [11] "anpp_grain"
## [12] "anpp_group"
## [13] "anpp_leaf"
## [14] "anpp_woody"
## [15] "asat"
## [16] "bgb"
## [17] "bgb_c"
## [18] "bgb_coarse"
## [19] "bgb_group"
## [20] "bgb_n"
## [21] "bgb_n_stock"
## [22] "bgb_p"
## [23] "bnpp"
## [24] "bpe"
## [25] "coarse_root_n"
## [26] "fine_root_biomass"
## [27] "fine_root_n"
## [28] "fine_root_production"
## [29] "fine_root_respiration"
## [30] "fine_root_turnover"
## [31] "gpp"
## [32] "gs"
## [33] "jmax"
## [34] "lai"
## [35] "leaf_area_eco"
## [36] "leaf_biomass_eco"
## [37] "leaf_biomass_plant"
## [38] "leaf_biomass_leaf"
## [39] "leaf_c"
## [40] "leaf_cn"
## [41] "leaf_k_mass"
## [42] "leaf_litter_decomposition"
## [43] "leaf_litterfall"
## [44] "leaf_n_area"
## [45] "leaf_n_mass"
## [46] "leaf_n_stock"
## [47] "leaf_np"
## [48] "leaf_p_mass"
## [49] "litter_biomass"
## [50] "litter_decomposition"
## [51] "litter_k"
## [52] "litter_n"
## [53] "litter_p"
## [54] "litterfall"
## [55] "lwp"
## [56] "mb"
## [57] "mbc"
## [58] "mbcn"
## [59] "mbn"
## [60] "nee"
## [61] "nep"
## [62] "npp"
## [63] "r_eco"
## [64] "r_root"
## [65] "r_soil"
## [66] "r_soilh"
## [67] "root_length"
## [68] "root_n_uptake"
## [69] "root_p_uptake"
## [70] "root_production"
## [71] "root_shoot_ratio"
## [72] "soc"
## [73] "soil_denitrification"
## [74] "soil_drp"
## [75] "soil_gross_nitrification"
## [76] "soil_n_leaching"
## [77] "soil_net_ammonification"
## [78] "soil_net_n_mineralization"
## [79] "soil_net_nitrification"
## [80] "soil_nh4"
## [81] "soil_nh4-n"
## [82] "soil_nh4-n_supply_rate"
## [83] "soil_nitrification"
## [84] "soil_no3"
## [85] "soil_no3_leaching"
## [86] "soil_no3-n"
## [87] "soil_no3-n_supply_rate"
## [88] "soil_np"
## [89] "soil_p"
## [90] "soil_p_supply_rate"
## [91] "soil_pbray"
## [92] "soil_phkcl"
## [93] "soil_phw"
## [94] "soil_potential_net_n_mineralization"
## [95] "soil_total_c"
## [96] "soil_total_c_min_layer"
## [97] "soil_total_c_org_layer"
## [98] "soil_total_c_profile"
## [99] "soil_total_cn"
## [100] "soil_total_cn_min_layer"
## [101] "soil_total_cn_org_layer"
## [102] "soil_total_n"
## [103] "soil_total_p"
## [104] "som"
## [105] "som_cn"
## [106] "stem_biomass"
## [107] "stem_n"
## [108] "swc"
## [109] "total_biomass"
## [110] "total_biomass_group"
## [111] "total_biomass_n_stock"
## [112] "total_biomass_p_stock"
## [113] "vcmax"
## [114] "wood_n"
## [115] "wue"
## [116] "wue_eco"Explore data availability in MESI dataset for P-fertilisation experiments
explore_mesi_pfert <- df_mesi %>%
# fertilisation experiments only
filter(treatment == "f") %>%
# P-fertilisation only (without N or K addition)
filter(npk == "_010")
head(explore_mesi_pfert)## # A tibble: 6 × 60
## db id duplicat…¹ citat…² respo…³ site study exp lat lon eleva…⁴
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 sichuan s172 s172 ohallo… agb -15.… -15.… -15.… -15.4 23.2 NA
## 2 sichuan s3497 s3497 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 3 sichuan s3500 s3500 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 4 sichuan s3503 s3503 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 5 sichuan s175 s175 ohallo… agb -18.… -18.… -18.… -18.7 25.5 NA
## 6 sichuan s168 s168 ries_a… agb -18.… -18.… -18.… -18.7 25.5 NA
## # … with 49 more variables: mat <dbl>, map <dbl>, ecosystem_type <chr>,
## # vegetation_type <chr>, experiment_type <chr>, community_type <chr>,
## # dominant_species <chr>, growth_form <chr>, age <dbl>,
## # disturbance_type <chr>, treatment <chr>, npk <chr>, w_t1 <chr>, c_c <dbl>,
## # c_t <dbl>, d_t <dbl>, d_t2 <dbl>, n_c <dbl>, n_t <dbl>, p_c <dbl>,
## # p_t <dbl>, k_c <dbl>, k_t <dbl>, i_c <dbl>, i_t <dbl>, i_t2 <dbl>,
## # s_c <dbl>, s_t <dbl>, w_t2 <dbl>, w_t3 <dbl>, start_year <dbl>, …## How many experiments?
length(unique(explore_mesi_pfert$citation))## [1] 126## What traits are available?
unique(explore_mesi_pfert$response)## [1] "agb" "agb_grain_n" "agb_grain_p"
## [4] "agb_group" "agb_height" "agb_n"
## [7] "agb_ndvi" "agb_p" "agb_p_stock"
## [10] "agb_pod_n" "agb_pod_p" "anpp"
## [13] "anpp_grain" "anpp_group" "bgb"
## [16] "bgb_n" "bgb_p" "bpe"
## [19] "fine_root_biomass" "fine_root_production" "gpp"
## [22] "gs" "lai" "leaf_biomass_eco"
## [25] "leaf_c" "leaf_n_mass" "leaf_np"
## [28] "leaf_p_mass" "litter_decomposition" "lwp"
## [31] "mbc" "mbn" "nee"
## [34] "r_root" "r_soil" "root_n_uptake"
## [37] "root_shoot_ratio" "soil_nh4-n_supply_rate" "soil_no3-n_supply_rate"
## [40] "soil_np" "soil_p" "soil_p_supply_rate"
## [43] "soil_total_c" "soil_total_n" "total_biomass"
## [46] "total_biomass_group" "total_biomass_n_stock" "total_biomass_np"
## [49] "total_biomass_p_stock" "wood_n" "wood_p"Note: no P fertilisation experiments include leaf-level physiology measurements beyond stomatal conductance. Leaf nutrient content is available
Select variables
use_response_p <- c("anpp",
"agb",
"bgb",
"fine_root_biomass",
"root_production",
"fine_root_production",
"root_shoot_ratio",
"leaf_n_mass",
"leaf_p_mass",
"lai",
"gpp",
"soil_no3-n",
"soil_nh4-n",
"soil_nh4",
"soil_no3",
"soil_solution_nh4",
"soil_solution_no3",
"root_n_uptake",
"root_nh4_uptake",
"root_no3_uptake"
)
pfert_mesi_responses <- explore_mesi_pfert %>%
filter(response %in% use_response_p) %>%
mutate(myvar = response) %>%
mutate(myvar = ifelse(myvar %in% c("agb",
"agb_coarse"),
"agb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("bgb",
"fine_root_biomass",
"coarse_root_c_stock",
"bgb_coarse"),
"bgb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("bgb",
"fine_root_biomass"),
"bgb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("root_production",
"fine_root_production",
"coarse_root_production"),
"root_production",
myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("root_n_uptake",
"root_nh4_uptake",
"root_no3_uptake"),
"n_uptake",
myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("soil_no3-n",
"soil_nh4-n",
"soil_nh4",
"soil_no3",
"soil_solution_nh4",
"soil_solution_no3"),
"n_inorg",
myvar))
use_vars_p <- unique(pfert_mesi_responses$myvar)Analysis
Calculate "ROM" - the log transformed ratio of means
(Hedges et al., 1999; Lajeunesse, 2011) for each observation pair
(ambient and elevated).
pfert_mesi_responses2 <- pfert_mesi_responses %>%
## keep only essential variables and drop rows containing missing values for
## essential variables
select(id, duplicate_id, exp, myvar, treatment, sampling_year,
x_t, x_c, sd_t, sd_c, rep_t, rep_c) %>%
## Get logarithm of response ratio and its variance
metafor::escalc(
measure = "ROM",
m1i = x_t, sd1i = sd_t, n1i = rep_t,
m2i = x_c, sd2i = sd_c, n2i = rep_c,
data = .,
append = TRUE,
var.names = c("logr", "logr_var")
) %>%
## to keep the output readable from the console output
as_tibble() %>%
## get standard error
mutate( logr_se = sqrt(logr_var) / sqrt(rep_t) )
head(pfert_mesi_responses2)## # A tibble: 6 × 15
## id duplicate_id exp myvar treat…¹ sampl…² x_t x_c sd_t sd_c rep_t
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 s172 s172 -15.… agb f 2 6.29 9.43 5.76 5.24 4
## 2 s3497 s3497 17.2… agb f 3 488. 106. 43.7 47.3 3
## 3 s3500 s3500 17.2… agb f 3 418. 200. 83.7 65.5 3
## 4 s3503 s3503 17.2… agb f 3 247. 80.1 43.7 29.1 3
## 5 s175 s175 -18.… agb f 2 3.67 4.19 5.76 6.29 4
## 6 s168 s168 -18.… agb f 1 80.0 51.1 19.6 17.1 4
## # … with 4 more variables: rep_c <dbl>, logr <dbl>, logr_var <dbl>,
## # logr_se <dbl>, and abbreviated variable names ¹treatment, ²sampling_year# Aggregate all measurements (multiple years, sampling dates and plots) by experiment (and response variable - although here only one) for meta analysis.
pfert_mesi_responses3 <- pfert_mesi_responses2 %>%
# suggested addition by Kevin, email 02.10.2023 10:03
dplyr::distinct(duplicate_id, x_t, x_c, .keep_all = TRUE) |>
filter(!is.na(logr_var) & !is.na(logr)) %>%
# re-create ID (common ID per experiment and response variable)
select(-id) %>%
mutate( id = paste(exp, myvar, sep = "_XXX_")) %>%
MAd::agg(
id = id,
es = logr,
var = logr_var,
cor = 1.0,
method = "BHHR",
data = .
) %>%
## to keep the output readable from the console output
as_tibble() %>%
# separate ID again for ease of data use
mutate( id = str_split(id, "_XXX_") ) %>%
mutate( exp = purrr::map_chr(id, 1),
myvar = purrr::map_chr(id, 2) ) %>%
## rename again
select(exp, myvar, logr = es, logr_var = var) %>%
## add number of observations (sum of plots and repeated samplings)
left_join(
pfert_mesi_responses2 %>%
group_by(exp, myvar, treatment) %>%
summarise(n_c = sum(rep_c), n_t = sum(rep_t)),
by = c("exp", "myvar")
) %>%
## get standard error. Verify if number available observations are identical
## for ambient and elevated. Use N from control here (n_c).
mutate( logr_se = sqrt(logr_var) / sqrt(n_c) )## `summarise()` has grouped output by 'exp', 'myvar'. You can override using the
## `.groups` argument.head(pfert_mesi_responses3)## # A tibble: 6 × 8
## exp myvar logr logr_var treatment n_c n_t logr_se
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 -15.44_23.25_f agb -0.405 0.287 f 4 4 0.268
## 2 17.25_-88.77_f agb 1.13 0.0572 f 9 9 0.0797
## 3 -18.66_25.5_f agb 0.158 0.418 f 8 8 0.229
## 4 -2.98_-47.52_f agb 0.269 0.000881 f 4 4 0.0148
## 5 22.13_-159.63_f agb 0.233 0.0368 f 11 11 0.0579
## 6 -22.283_117.666_f agb -0.194 0.00676 f 6 6 0.0336Meta-analysis
Aggregate log-ratios across multiple experiments, taking into account their respective variance and using the experiment identity as a grouping factor for random intercepts.
source("R/analyse_meta.R")
out_p <- purrr::map(as.list(use_vars_p),
~analyse_meta(pfert_mesi_responses3 %>%
rename(var = myvar), nam_target = .))
names(out_p) <- use_vars_p
df_box_p <- purrr::map_dfr(out_p, "df_box") |>
left_join(
pfert_mesi_responses3 |>
group_by(myvar) |>
summarise(logr_min = min(logr), logr_max = max(logr)) |>
rename(var = myvar),
by = "var"
)
saveRDS(df_box_p, file = paste0(here::here(), "df_box_mesi_pfert_eperkowski.rds"))Final data size
Number of data points (plot-level measurements) per variable:
pfert_mesi_responses3 %>%
group_by(myvar) %>%
summarise(n_plots = sum(n_c, na.rm = TRUE), n_exp = n()) %>%
rename_("Variable"="myvar", "N observations"="n_plots", "N experiments"="n_exp") %>%
knitr::kable()## Warning: `rename_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `rename()` instead.
## ℹ The deprecated feature was likely used in the dplyr package.
## Please report the issue at <]8;;https://github.com/tidyverse/dplyr/issueshttps://github.com/tidyverse/dplyr/issues]8;;>.| Variable | N observations | N experiments |
|---|---|---|
| agb | 839 | 103 |
| anpp | 82 | 13 |
| bgb | 189 | 30 |
| gpp | 57 | 3 |
| leaf_n_mass | 8 | 4 |
| leaf_p_mass | 8 | 4 |
| n_uptake | 24 | 3 |
| root_production | 24 | 1 |
| root_shoot_ratio | 15 | 4 |
Number of data points (plot-level measurements) per experiment:
pfert_mesi_responses3 %>%
group_by(exp) %>%
summarise(n_plots = sum(n_c), n_exp = n()) %>%
rename_("Experiment"="exp", "N observations"="n_plots", "N experiments"="n_exp") %>%
knitr::kable()| Experiment | N observations | N experiments |
|---|---|---|
| -15.44_23.25_f | 4 | 1 |
| -18.66_25.5_f | 8 | 1 |
| -2.98_-47.52_f | 4 | 1 |
| -22.283_117.666_f | 6 | 1 |
| -22.41_21.71_f | 4 | 1 |
| -22.78_31.25_f | 16 | 2 |
| -23.75_31.43_f | 16 | 2 |
| -24.17_21.89_f | 4 | 1 |
| -24.4_31.75_f | 16 | 2 |
| -25.12_31.23_f | 16 | 2 |
| -25.29_31.91_f | 16 | 2 |
| -3.5_36_f | 6 | 1 |
| -3.95_-79.03_f | 12 | 2 |
| -33.35_150.3_f | 6 | 2 |
| -34.58_-58.48_f | 9 | 1 |
| -35.73_-58.05_f | 12 | 1 |
| 0.28_37.88_f | 4 | 1 |
| 17.25_-88.77_f | 9 | 1 |
| 19.6_-155.33_fp | 3 | 1 |
| 22.13_-159.63_f | 11 | 1 |
| 26.52_109.78_fp | 3 | 1 |
| 31.3_-81.28_f | 4 | 1 |
| 31.37_90.02_f | 24 | 1 |
| 31.55_-81.78_f | 4 | 1 |
| 32.54_-116.7_fp | 12 | 1 |
| 34.92_102.88_f2p | 5 | 1 |
| 34.92_102.88_f3p | 5 | 1 |
| 34.92_102.88_fp | 5 | 1 |
| 37.25_-121.75_forb_fp | 3 | 1 |
| 37.25_-121.75_grass_fp | 3 | 1 |
| 37.48_101.2_fp | 12 | 2 |
| 37.55_-122.3_f | 3 | 1 |
| 37.6_101.32_fp | 16 | 2 |
| 37.87_-122.52_f | 18 | 2 |
| 39.15_-86.52_f | 3 | 1 |
| 39.25_-121.28_fp | 20 | 1 |
| 39.51_-108.19_fp | 4 | 1 |
| 39.75_-74.75_fp | 4 | 1 |
| 41.35_36.25_fp | 4 | 1 |
| 41.35_36.25_fp2 | 4 | 1 |
| 41.62_-71.32_fp | 8 | 2 |
| 42.28_-85.58_f | 6 | 2 |
| 42.58_122.21_fp | 12 | 2 |
| 44.8_116.03_fp | 6 | 2 |
| 47.57_7.6_f | 6 | 1 |
| 51.85_5.62_fp | 3 | 1 |
| 52.07_5.58_fp | 3 | 1 |
| 52.37_5.1_fp | 6 | 1 |
| 52.5_5.7_fp | 6 | 1 |
| 54.63_8.83_fp | 5 | 1 |
| 69.43_-133.02_fp | 10 | 1 |
| 69.43_-133.02_fp2 | 10 | 1 |
| 9.6_-79.5_f | 3 | 1 |
| alpflix_fp | 60 | 1 |
| andorra_mire_fp | 5 | 1 |
| bcnm_fp | 20 | 1 |
| bennekom_drained_fp | 5 | 1 |
| bennekom_undrained_fp | 5 | 1 |
| bordeaux_f | 12 | 1 |
| bordeaux_f2 | 12 | 1 |
| bordeaux_fp | 8 | 2 |
| buitengoor_1992_fp | 5 | 1 |
| buitengoor_1993_fp | 5 | 1 |
| daqinggou_fp | 12 | 2 |
| drentsche_aa_drained_fp | 5 | 1 |
| drentsche_aa_wet_fp | 5 | 1 |
| ewenke_f_p | 6 | 1 |
| gdfr_f | 10 | 1 |
| gusewell_s1_fp | 4 | 1 |
| gusewell_s2_fp | 4 | 1 |
| gusewell_t1_fp | 4 | 1 |
| gusewell_t2_fp | 4 | 1 |
| gusewell_v1_fp | 4 | 1 |
| gusewell_v2_fp | 4 | 1 |
| gusewell_v3_fp | 4 | 1 |
| gusewell_w1_p | 4 | 1 |
| gusewell_w2_p | 4 | 1 |
| gusewell_w3_p | 4 | 1 |
| gusewell_w4_p | 4 | 1 |
| hol_kortenhoef_fp | 5 | 1 |
| ibge_fp | 28 | 2 |
| imgers_ng_2006np_fp | 6 | 1 |
| indoneisa_f | 36 | 2 |
| katelijne_2016_fp | 10 | 2 |
| katelijne_2017_f | 5 | 1 |
| katelijne_2017_f2 | 10 | 1 |
| katelijne_2017_f3 | 15 | 2 |
| kisa_grahamiana_f | 20 | 1 |
| luneburg_field_2006_fp | 10 | 1 |
| luneburg_field_2008_fp | 10 | 1 |
| luneburg_gh_drought_molinia_fp | NA | 3 |
| luneburg_gh_fert_molinia_fp | 10 | 1 |
| michigan_underc_bog_fp | 39 | 3 |
| michigan_underc_intermfen_fp | 37 | 3 |
| michigan_underc_richfen_fp | 40 | 3 |
| molenpolder_fp | 5 | 1 |
| nashfield_pooled_fp | 8 | 2 |
| newyork_tnp_fp | 45 | 2 |
| niwot_ridge2_dm_fp | 10 | 2 |
| niwot_ridge2_wm_fp | 10 | 2 |
| pituffik_fp | 12 | 1 |
| sanjiang_mire_pfert_fp | 9 | 1 |
| sanjiang_mire_pfert_fp2 | 9 | 1 |
| sanjiang_mire_pfert_fp3 | 9 | 1 |
| sanpedro_fp | 4 | 2 |
| schiermonnikoog_old_fp | 12 | 1 |
| schiermonnikoog_old_fp2 | 12 | 1 |
| schiermonnikoog_young_fp | 12 | 1 |
| schiermonnikoog_young_fp2 | 12 | 1 |
| sedgwick_fp | 3 | 1 |
| shaaxi_330_bungeana_fp | 3 | 1 |
| shaaxi_330_davarica_fp | 3 | 1 |
| shaaxi_330_fp | 3 | 1 |
| shaaxi_6_bungeana_fp | 3 | 1 |
| shaaxi_6_davarica_fp | 3 | 1 |
| shaaxi_6_fp | 3 | 1 |
| tambopata_fp | 4 | 2 |
| teberda_f | 8 | 1 |
| tono_fp | 4 | 2 |
| toolik_nonacidic_fp | 9 | 3 |
| wayqecha_fp | 4 | 2 |
| westbroek_polder_fp | 5 | 1 |
| yucatan_marsh_highsalinity_fp | 30 | 1 |
| yucatan_marsh_lowsalinity_fp | 30 | 1 |
| yucatan_marsh_mediumsalinity_fp | 30 | 1 |
| zwarte_beek_drained_fp | 5 | 1 |
| zwarte_beek_wet_fp | 5 | 1 |
Some quick plots:
ggplot(data = pfert_mesi_responses3, aes(x = myvar, y = logr)) +
geom_jitter(color = rgb(0,0,0,0.3),
aes( size = 1/logr_se ),
position = position_jitter(w = 0.2, h = 0),
show.legend = FALSE) +
geom_crossbar( data = df_box_p %>% drop_na(var),
aes(x = var, y = middle, ymin = ymin, ymax = ymax),
fill = "royalblue",
color = "royalblue4",
alpha = 0.6,
width = 0.5 ) +
geom_hline( yintercept = 0.0, linewidth = 0.5, linetype = "dotted" ) +
#scale_x_discrete("", labels = mylabl) +
labs(x = "",
y = "Log response to P addition") +
coord_flip() +
theme_classic(base_size = 18)## Warning: Removed 3 rows containing missing values (`geom_point()`).
Explore data availability in MESI dataset NxP-fertilisation experiments
explore_mesi_npfert <- df_mesi %>%
# fertilisation experiments only
filter(treatment == "f") %>%
# NxP-fertilisation only (without K addition)
filter(npk == "_110")
head(explore_mesi_npfert)## # A tibble: 6 × 60
## db id duplicat…¹ citat…² respo…³ site study exp lat lon eleva…⁴
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 sichuan s171 s171 ohallo… agb -15.… -15.… -15.… -15.4 23.2 NA
## 2 sichuan s3499 s3499 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 3 sichuan s3502 s3502 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 4 sichuan s3505 s3505 rejman… agb 17.2… 17.2… 17.2… 17.2 -88.8 NA
## 5 sichuan s174 s174 ohallo… agb -18.… -18.… -18.… -18.7 25.5 NA
## 6 sichuan s169 s169 ries_a… agb -18.… -18.… -18.… -18.7 25.5 NA
## # … with 49 more variables: mat <dbl>, map <dbl>, ecosystem_type <chr>,
## # vegetation_type <chr>, experiment_type <chr>, community_type <chr>,
## # dominant_species <chr>, growth_form <chr>, age <dbl>,
## # disturbance_type <chr>, treatment <chr>, npk <chr>, w_t1 <chr>, c_c <dbl>,
## # c_t <dbl>, d_t <dbl>, d_t2 <dbl>, n_c <dbl>, n_t <dbl>, p_c <dbl>,
## # p_t <dbl>, k_c <dbl>, k_t <dbl>, i_c <dbl>, i_t <dbl>, i_t2 <dbl>,
## # s_c <dbl>, s_t <dbl>, w_t2 <dbl>, w_t3 <dbl>, start_year <dbl>, …## How many experiments?
length(unique(explore_mesi_npfert$citation))## [1] 133## What traits are measured?
unique(explore_mesi_npfert$response)## [1] "agb" "agb_c"
## [3] "agb_cn" "agb_cp"
## [5] "agb_group" "agb_height"
## [7] "agb_n" "agb_ndvi"
## [9] "agb_p" "agb_p_stock"
## [11] "anpp" "bgb"
## [13] "bgb_c" "bgb_n"
## [15] "bgb_p" "bpe"
## [17] "fine_root_biomass" "fine_root_production"
## [19] "gpp" "grain_c"
## [21] "grain_n" "gs"
## [23] "leaf_biomass_eco" "leaf_c"
## [25] "leaf_n_mass" "leaf_np"
## [27] "leaf_p_mass" "litter_decomposition"
## [29] "lwp" "mbc"
## [31] "nee" "r_root"
## [33] "r_soil" "root_n"
## [35] "root_n_uptake" "root_shoot_ratio"
## [37] "soc" "soil_drp"
## [39] "soil_n2o_flux" "soil_nh4"
## [41] "soil_nh4-n_supply_rate" "soil_no3_leaching"
## [43] "soil_no3-n_supply_rate" "soil_np"
## [45] "soil_p" "soil_p_supply_rate"
## [47] "soil_pbray" "soil_phkcl"
## [49] "soil_total_c" "soil_total_cn"
## [51] "soil_total_n" "som_cn"
## [53] "stem_biomass" "swc"
## [55] "total_biomass" "total_biomass_group"
## [57] "total_biomass_np" "total_biomass_production"Note: no NxP fertilisation experiments in MESI database include leaf-level physiology measurements beyond stomatal conductance (as with P fertilisation experiments). Leaf nutrient content is available, however.
Select variables
use_response_np <- c("anpp",
"agb",
"bgb",
"fine_root_biomass",
"root_production",
"fine_root_production",
"root_shoot_ratio",
"leaf_n_mass",
"leaf_p_mass",
"lai",
"gpp",
"soil_no3-n",
"soil_nh4-n",
"soil_nh4",
"soil_no3",
"soil_solution_nh4",
"soil_solution_no3",
"root_n_uptake",
"root_nh4_uptake",
"root_no3_uptake"
)
npfert_mesi_responses <- explore_mesi_npfert %>%
filter(response %in% use_response_np) %>%
mutate(myvar = response) %>%
mutate(myvar = ifelse(myvar %in% c("agb",
"agb_coarse"),
"agb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("bgb",
"fine_root_biomass",
"coarse_root_c_stock",
"bgb_coarse"),
"bgb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("bgb",
"fine_root_biomass"),
"bgb", myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("root_production",
"fine_root_production",
"coarse_root_production"),
"root_production",
myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("root_n_uptake",
"root_nh4_uptake",
"root_no3_uptake"),
"n_uptake",
myvar)) %>%
mutate(myvar = ifelse(myvar %in% c("soil_no3-n",
"soil_nh4-n",
"soil_nh4",
"soil_no3",
"soil_solution_nh4",
"soil_solution_no3"),
"n_inorg",
myvar))
use_vars_np <- unique(npfert_mesi_responses$myvar)Analysis
Calculate "ROM" - the log transformed ratio of means
(Hedges et al., 1999; Lajeunesse, 2011) for each observation pair
(ambient and elevated).
npfert_mesi_responses2 <- npfert_mesi_responses %>%
## keep only essential variables and drop rows containing missing values for
## essential variables
select(id, duplicate_id, exp, myvar, treatment, sampling_year,
x_t, x_c, sd_t, sd_c, rep_t, rep_c) %>%
## Get logarithm of response ratio and its variance
metafor::escalc(
measure = "ROM",
m1i = x_t, sd1i = sd_t, n1i = rep_t,
m2i = x_c, sd2i = sd_c, n2i = rep_c,
data = .,
append = TRUE,
var.names = c("logr", "logr_var")
) %>%
## to keep the output readable from the console output
as_tibble() %>%
## get standard error
mutate( logr_se = sqrt(logr_var) / sqrt(rep_t) )
head(npfert_mesi_responses2)## # A tibble: 6 × 15
## id duplicate_id exp myvar treat…¹ sampl…² x_t x_c sd_t sd_c rep_t
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 s171 s171 -15.… agb f 2 8.38 9.43 5.76 5.24 4
## 2 s3499 s3499 17.2… agb f 3 506. 106. 40.0 47.3 3
## 3 s3502 s3502 17.2… agb f 3 364. 200. 76.4 65.5 3
## 4 s3505 s3505 17.2… agb f 3 262. 80.1 76.4 29.1 3
## 5 s174 s174 -18.… agb f 2 9.96 4.19 6.29 6.29 4
## 6 s169 s169 -18.… agb f 1 163. 51.1 64.3 17.1 4
## # … with 4 more variables: rep_c <dbl>, logr <dbl>, logr_var <dbl>,
## # logr_se <dbl>, and abbreviated variable names ¹treatment, ²sampling_year# Aggregate all measurements (multiple years, sampling dates and plots) by experiment (and response variable - although here only one) for meta analysis.
npfert_mesi_responses3 <- npfert_mesi_responses2 %>%
# suggested addition by Kevin, email 02.10.2023 10:03
dplyr::distinct(duplicate_id, x_t, x_c, .keep_all = TRUE) |>
filter(!is.na(logr_var) & !is.na(logr)) %>%
# re-create ID (common ID per experiment and response variable)
select(-id) %>%
mutate( id = paste(exp, myvar, sep = "_XXX_")) %>%
MAd::agg(
id = id,
es = logr,
var = logr_var,
cor = 1.0,
method = "BHHR",
data = .
) %>%
## to keep the output readable from the console output
as_tibble() %>%
# separate ID again for ease of data use
mutate( id = str_split(id, "_XXX_") ) %>%
mutate( exp = purrr::map_chr(id, 1),
myvar = purrr::map_chr(id, 2) ) %>%
## rename again
select(exp, myvar, logr = es, logr_var = var) %>%
## add number of observations (sum of plots and repeated samplings)
left_join(
npfert_mesi_responses2 %>%
group_by(exp, myvar, treatment) %>%
summarise(n_c = sum(rep_c), n_t = sum(rep_t)),
by = c("exp", "myvar")
) %>%
## get standard error. Verify if number available observations are identical
## for ambient and elevated. Use N from control here (n_c).
mutate( logr_se = sqrt(logr_var) / sqrt(n_c) )## `summarise()` has grouped output by 'exp', 'myvar'. You can override using the
## `.groups` argument.head(npfert_mesi_responses3)## # A tibble: 6 × 8
## exp myvar logr logr_var treatment n_c n_t logr_se
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 -15.44_23.25_f agb -0.118 0.195 f 4 4 0.221
## 2 17.25_-88.77_f agb 1.11 0.0634 f 9 9 0.0839
## 3 -18.66_25.5_f agb 1.01 0.287 f 8 8 0.190
## 4 -2.98_-47.52_f agb 0.245 0.000560 f 4 4 0.0118
## 5 22.13_-159.63_f agb 0.571 0.0281 f 11 11 0.0506
## 6 -22.283_117.666_f agb 0.213 0.00928 f 6 6 0.0393Meta-analysis
Aggregate log-ratios across multiple experiments, taking into account their respective variance and using the experiment identity as a grouping factor for random intercepts.
source("R/analyse_meta.R")
out_np <- purrr::map(as.list(use_vars_np),
~analyse_meta(npfert_mesi_responses3 %>%
rename(var = myvar), nam_target = .))
names(out_np) <- use_vars_np
df_box_np <- purrr::map_dfr(out_np, "df_box") |>
left_join(
npfert_mesi_responses3 |>
group_by(myvar) |>
summarise(logr_min = min(logr), logr_max = max(logr)) |>
rename(var = myvar),
by = "var"
)
saveRDS(df_box_np, file = paste0(here::here(), "df_box_mesi_npfert_eperkowski.rds"))Final data size
Number of data points (plot-level measurements) per variable:
npfert_mesi_responses3 %>%
group_by(myvar) %>%
summarise(n_plots = sum(n_c, na.rm = TRUE), n_exp = n()) %>%
rename_("Variable"="myvar", "N observations"="n_plots", "N experiments"="n_exp") %>%
knitr::kable()| Variable | N observations | N experiments |
|---|---|---|
| agb | 1036 | 108 |
| anpp | 646 | 18 |
| bgb | 328 | 42 |
| gpp | 57 | 3 |
| leaf_n_mass | 8 | 4 |
| leaf_p_mass | 8 | 4 |
| n_inorg | 20 | 1 |
| n_uptake | 36 | 4 |
| root_production | 24 | 1 |
| root_shoot_ratio | 10 | 3 |
Number of data points (plot-level measurements) per experiment:
npfert_mesi_responses3 %>%
group_by(exp) %>%
summarise(n_plots = sum(n_c), n_exp = n()) %>%
rename_("Experiment"="exp", "N observations"="n_plots", "N experiments"="n_exp") %>%
knitr::kable()| Experiment | N observations | N experiments |
|---|---|---|
| -15.44_23.25_f | 4 | 1 |
| -18.66_25.5_f | 8 | 1 |
| -2.98_-47.52_f | 4 | 1 |
| -22.283_117.666_f | 6 | 1 |
| -22.41_21.71_f | 4 | 1 |
| -22.78_31.25_f | 16 | 2 |
| -23.75_31.43_f | 16 | 2 |
| -24.17_21.89_f | 4 | 1 |
| -24.4_31.75_f | 16 | 2 |
| -25.12_31.23_f | 16 | 2 |
| -25.29_31.91_f | 16 | 2 |
| -3.5_36_f | 6 | 1 |
| -3.95_-79.03_f | 12 | 2 |
| -34.58_-58.48_f | 9 | 1 |
| 0.28_37.88_f | 4 | 1 |
| 17.25_-88.77_f | 9 | 1 |
| 19.6_-155.33_fnp | 3 | 1 |
| 22.13_-159.63_f | 11 | 1 |
| 26.52_109.78_fn2p | 3 | 1 |
| 31.3_-81.28_f | 4 | 1 |
| 31.37_90.02_f | 48 | 1 |
| 31.42_-88.45_f | 8 | 2 |
| 31.55_-81.78_f | 4 | 1 |
| 32.54_-116.7_fnp | 12 | 1 |
| 34.92_102.88_f2np | 5 | 1 |
| 34.92_102.88_f3np | 5 | 1 |
| 34.92_102.88_fnp | 5 | 1 |
| 35.97_101.88_f | 40 | 2 |
| 37.25_-121.75_forb_fnp | 3 | 1 |
| 37.25_-121.75_grass_fnp | 3 | 1 |
| 37.48_101.2_fnp | 12 | 2 |
| 37.55_-122.3_f | 3 | 1 |
| 37.6_101.32_fnp | 16 | 2 |
| 37.87_-122.52_f | 14 | 2 |
| 39.15_-86.52_f | 3 | 1 |
| 39.25_-121.28_fnp | 20 | 1 |
| 39.51_-108.19_fnp | 4 | 1 |
| 41.35_36.25_fn2p | 4 | 1 |
| 41.35_36.25_fn2p2 | 4 | 1 |
| 41.35_36.25_fn3p | 4 | 1 |
| 41.35_36.25_fn3p2 | 4 | 1 |
| 41.35_36.25_fnp | 4 | 1 |
| 41.35_36.25_fnp2 | 4 | 1 |
| 41.62_-71.32_fnp | 8 | 2 |
| 42.58_122.21_fnp | 12 | 2 |
| 44.8_116.03_fnp | 24 | 2 |
| 51.85_5.62_fnp | 3 | 1 |
| 52.07_5.58_fnp | 3 | 1 |
| 52.37_5.1_fnp | 6 | 1 |
| 52.5_5.7_fnp | 6 | 1 |
| 54.63_8.83_fnp | 5 | 1 |
| 68.2_-149.6_f | 16 | 2 |
| 68.38_-104.54_f | 4 | 1 |
| 69.43_-133.02_fn2p2 | 10 | 1 |
| 69.43_-133.02_fnp | 10 | 1 |
| 9.6_-79.5_f | 3 | 1 |
| alpflix_fnp | 60 | 1 |
| andorra_mire_fnp | 5 | 1 |
| bcnm_fnp | 20 | 1 |
| bennekom_drained_fnp | 5 | 1 |
| bennekom_undrained_fnp | 5 | 1 |
| buitengoor_1992_fnp | 5 | 1 |
| damxung_2015_f | 476 | 1 |
| damxung_f | 120 | 2 |
| daqinggou_fnp | 12 | 2 |
| drentsche_aa_drained_fnp | 5 | 1 |
| drentsche_aa_wet_fnp | 5 | 1 |
| duolun15_fn2p | 4 | 1 |
| duolun15_fnp | 4 | 1 |
| escambia_county_f | 8 | 2 |
| ewenke_f_np | 6 | 1 |
| gdfr_f | 10 | 1 |
| gusewell_s1_fnp | 4 | 1 |
| gusewell_s2_fnp | 4 | 1 |
| gusewell_t1_fnp | 4 | 1 |
| gusewell_t2_fnp | 4 | 1 |
| gusewell_v1_fnp | 4 | 1 |
| gusewell_v2_fnp | 4 | 1 |
| gusewell_v3_fnp | 4 | 1 |
| gusewell_w1_fnp | 4 | 1 |
| gusewell_w2_fnp | 4 | 1 |
| gusewell_w3_fnp | 4 | 1 |
| gusewell_w4_fnp | 4 | 1 |
| haibei_fn1p | 33 | 1 |
| haibei_fn1pp | 33 | 1 |
| haibei_fn2p | 33 | 1 |
| haibei_fn2pp | 33 | 1 |
| haibei_fn3p | 33 | 1 |
| haibei_fn3pp | 33 | 1 |
| harbin_f | 12 | 1 |
| hawke_f | 20 | 1 |
| ibge_fnp | 28 | 2 |
| imgers_ng_2006np_fn2 | 6 | 1 |
| imgers_ng_2006np_fn3 | 6 | 1 |
| imgers_ng_2006np_fn4 | 6 | 1 |
| imgers_ng_2006np_fn5 | 6 | 1 |
| imgers_ng_2006np_fn6 | 6 | 1 |
| imgers_ng_2006np_fp2 | 6 | 1 |
| imgers_ng_2006np_fp3 | 6 | 1 |
| imgers_ng_2006np_fp4 | 6 | 1 |
| imgers_ng_2006np_fp5 | 6 | 1 |
| imgers_ng_2006np_fp6 | 6 | 1 |
| katelijne_2016_fnp | 10 | 2 |
| katelijne_2017_f | 5 | 1 |
| katelijne_2017_f2 | 5 | 1 |
| luneburg_field_2006_fnp | 10 | 1 |
| luneburg_field_2008_fnp | 10 | 1 |
| luneburg_gh_drought_molinia_fnp | NA | 3 |
| luneburg_gh_fert_molinia_fnp | 10 | 1 |
| michigan_underc_bog_fnp | 39 | 3 |
| michigan_underc_intermfen_fnp | 37 | 3 |
| michigan_underc_richfen_fnp | 40 | 3 |
| nashfield_pooled_fnp | 8 | 2 |
| newyork_tnp_fnp | 45 | 2 |
| niwot_ridge2_dm_f2np | 5 | 1 |
| niwot_ridge2_dm_fnp | 10 | 2 |
| niwot_ridge2_wm_fnp | 15 | 2 |
| pituffik_fnp | 12 | 1 |
| rengen_fnp | 5 | 1 |
| sanpedro_fnp | 4 | 2 |
| schiermonnikoog_old_fn2p | 12 | 1 |
| schiermonnikoog_old_fn2p2 | 12 | 1 |
| schiermonnikoog_old_fnp | 12 | 1 |
| schiermonnikoog_old_fnp2 | 12 | 1 |
| schiermonnikoog_young_fn2p | 12 | 1 |
| schiermonnikoog_young_fn2p2 | 12 | 1 |
| schiermonnikoog_young_fnp | 12 | 1 |
| schiermonnikoog_young_fnp2 | 12 | 1 |
| sedgwick_fn2p | 3 | 1 |
| sedgwick_fnp | 3 | 1 |
| shaaxi_330_bungeana_fnp | 3 | 1 |
| shaaxi_330_davarica_fnp | 3 | 1 |
| shaaxi_330_fnp | 3 | 1 |
| shaaxi_6_bungeana_fn3p | 3 | 1 |
| shaaxi_6_davarica_fn3p | 3 | 1 |
| shaaxi_6_fn3p | 3 | 1 |
| tambopata_fnp | 4 | 2 |
| teberda_f | 8 | 1 |
| tono_fnp | 4 | 2 |
| toolik_acidic_1981_f | 4 | 1 |
| toolik_acidic_f | 20 | 1 |
| toolik_nonacidic_fnp | 13 | 3 |
| toolik_shrub_f | 4 | 2 |
| wayqecha_fnp | 4 | 2 |
| yucatan_marsh_highsalinity_fnp | 30 | 1 |
| yucatan_marsh_lowsalinity_fnp | 30 | 1 |
| yucatan_marsh_mediumsalinity_fnp | 30 | 1 |
| zwarte_beek_drained_fnp | 5 | 1 |
| zwarte_beek_wet_fnp | 5 | 1 |
Some quick plots:
ggplot(data = npfert_mesi_responses3, aes(x = myvar, y = logr)) +
geom_jitter(color = rgb(0,0,0,0.3),
aes( size = 1/logr_se ),
position = position_jitter(w = 0.2, h = 0),
show.legend = FALSE) +
geom_crossbar( data = df_box_np %>% drop_na(var),
aes(x = var, y = middle, ymin = ymin, ymax = ymax),
fill = "royalblue",
color = "royalblue4",
alpha = 0.6,
width = 0.5 ) +
geom_hline( yintercept = 0.0, linewidth = 0.5, linetype = "dotted" ) +
#scale_x_discrete("", labels = mylabl) +
labs(x = "",
y = "Log response to N*P addition") +
coord_flip() +
theme_classic(base_size = 18)## Warning: Removed 3 rows containing missing values (`geom_point()`).