Load packages and data

## Libraries
library(dplyr)       # for working with ease
library(tidyr)       # for working with ease
library(stringr)
library(ggplot2)
library(readr)
library(metafor)     # for meta analysis
library(MAd)         # for meta analysis

## MESI data
df.mesi <- read_csv("/Users/eaperkowski/git/mesi-db/data/mesi_main.csv")

## NutNet data from Cleland et al. (2019). Calculate aboveground biomass from root biomass and root mass fraction, then calculate root:shoot ratio
df.nutnet <- read_csv("/Users/eaperkowski/Desktop/Cleland_2019_RMF.csv") %>%
  rename(bgb = Rootsgperm2, rmf = rootmassfraction) %>%
  mutate(agb = (bgb/rmf) - bgb,
         root_shoot_ratio = bgb / agb)

Cleland et al. (2019) data cleaning and analysis

Clean data and calculate log RR for root mass fraction, belowground biomass, aboveground biomass, and root:shoot ratio for each site in a single step

df.nutnet.summary <- df.nutnet %>%
  filter(treatment_name == "N" | treatment_name == "Control") %>%
  group_by(site_code, treatment_name) %>%
  summarize(
    
    aridity = mean(Aridity, na.rm = TRUE),
    light = mean(Light, na.rm = TRUE),
    
    x_rmf = mean(rmf, na.rm = TRUE),
    sd_rmf = sd(rmf, na.rm = TRUE),
    n_rmf = length(rmf),
    
    x_bgb = mean(bgb, na.rm = TRUE),
    sd_bgb = sd(bgb, na.rm = TRUE),
    n_bgb = length(bgb),
    
    x_agb = mean(agb, na.rm = TRUE),
    sd_agb = sd(agb, na.rm = TRUE),
    n_agb = length(agb),
            
    x_rootshoot = mean(root_shoot_ratio, na.rm = TRUE),
    sd_rootshoot = sd(root_shoot_ratio, na.rm = TRUE),
    n_rootshoot = length(root_shoot_ratio)
    
    ) %>%
  
  pivot_wider(names_from = treatment_name, values_from = x_rmf:n_rootshoot) %>%
  rename(
    
    x_c_rmf = x_rmf_Control, x_t_rmf = x_rmf_N,
    sd_c_rmf = sd_rmf_Control, sd_t_rmf = sd_rmf_N,
    rep_c_rmf = n_rmf_Control, rep_t_rmf = n_rmf_N,
    
    x_c_bgb = x_bgb_Control, x_t_bgb = x_bgb_N,
    sd_c_bgb = sd_bgb_Control, sd_t_bgb = sd_bgb_N,
    rep_c_bgb = n_bgb_Control, rep_t_bgb = n_bgb_N,
         
    x_c_agb = x_agb_Control, x_t_agb = x_agb_N,
    sd_c_agb = sd_agb_Control, sd_t_agb = sd_agb_N,
    rep_c_agb = n_agb_Control, rep_t_agb = n_agb_N,   
         
    x_c_rootshoot = x_rootshoot_Control, x_t_rootshoot = x_rootshoot_N,
    sd_c_rootshoot = sd_rootshoot_Control, sd_t_rootshoot = sd_rootshoot_N,
    rep_c_rootshoot = n_rootshoot_Control, rep_t_rootshoot = n_rootshoot_N
    )

Calculate ROM for RMF data:

df.nutnet.rmf <- df.nutnet.summary[1:9] %>%
  rename(x_c = x_c_rmf,
         x_t = x_t_rmf,
         sd_c = sd_c_rmf,
         sd_t = sd_t_rmf,
         rep_c = rep_c_rmf,
         rep_t = rep_t_rmf) %>%
  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")) %>%
  mutate(logr_se = sqrt(logr_var) / sqrt(rep_t),
         myvar = "rmf")

Calculate ROM for belowground biomass:

df.nutnet.bgb <- df.nutnet.summary[, c(1:3, 10:15)] %>%
  rename(x_c = x_c_bgb, x_t = x_t_bgb,
         sd_c = sd_c_bgb, sd_t = sd_t_bgb,
         rep_c = rep_c_bgb, rep_t = rep_t_bgb) %>%
  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")) %>%
  mutate(logr_se = sqrt(logr_var) / sqrt(rep_t),
         myvar = "bgb")

Calculate ROM for aboveground biomass data:

df.nutnet.agb <- df.nutnet.summary[c(1:3, 16:21)] %>%
  rename(x_c = x_c_agb,
         x_t = x_t_agb,
         sd_c = sd_c_agb,
         sd_t = sd_t_agb,
         rep_c = rep_c_agb,
         rep_t = rep_t_agb) %>%
  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")) %>%
  mutate(logr_se = sqrt(logr_var) / sqrt(rep_t),
         myvar = "agb")

Calculate ROM for root:shoot ratio data:

df.nutnet.rootshoot <- df.nutnet.summary[c(1:3, 22:27)] %>%
  rename(x_c = x_c_rootshoot,
         x_t = x_t_rootshoot,
         sd_c = sd_c_rootshoot,
         sd_t = sd_t_rootshoot,
         rep_c = rep_c_rootshoot,
         rep_t = rep_t_rootshoot) %>%
  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")) %>%
  mutate(logr_se = sqrt(logr_var) / sqrt(rep_t),
         myvar = "root_shoot_ratio")

Merge lnRR data for all belowground traits, format for easy merge into MESI dataset

nutnet.total <- df.nutnet.bgb %>% full_join(df.nutnet.rmf) %>%
  full_join(df.nutnet.agb) %>% full_join(df.nutnet.rootshoot) %>%
  dplyr::select(exp = site_code, myvar, logr, logr_var, n_c = rep_c, n_t = rep_t,
                logr_se) %>%
  mutate(exp = str_c("nutnet_", exp),
         treatment = "f")
head(nutnet.total)
## 
##                  exp myvar    logr logr_var n_c n_t    logr_se treatment 
## 1     nutnet_bldr.us   bgb  0.9717   0.9436   2   2 0.68687273         f 
## 2     nutnet_bnch.us   bgb -0.0911   0.0238   3   3 0.08910833         f 
## 3   nutnet_bogong.au   bgb  0.8837   0.0453   3   3 0.12288166         f 
## 4 nutnet_burrawan.au   bgb -1.6745   0.0701   3   3 0.15283252         f 
## 5     nutnet_cbgb.us   bgb -0.0423   0.0445   6   6 0.08611377         f 
## 6     nutnet_cdcr.us   bgb -0.0114   0.1620   3   3 0.23239268         f

MESI dataset and merge with Cleland et al. (2019)

Follows exact approach from Beni’s analysis from Rpubs. Link: https://rpubs.com/stineb/analysis_mesi

Data subsetting

df0 <- df.mesi %>%
  
  ## fertilisation experiments only
  filter(treatment == "f") %>% 
  
  ## only experiments conducted in the field
  filter(experiment_type == "field")

## Data subsetting
use_exp <- df0 %>% 
  filter(!is.na(sampling_year)) %>% 
  group_by(exp) %>% 
  summarise(nyears = max(sampling_year)) %>% 
  filter(nyears >= 3) %>% 
  pull(exp)
# df1 <- df0 %>%
#   filter(exp %in% use_exp)
## xxx try: do not subset to minimum experiment length
df1 <- df0

Select variables:

use_response <- c("anpp", "agb", "bgb", "fine_root_biomass", "root_production", "fine_root_production", "root_shoot_ratio", "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")

df2 <- df1 %>% 
  filter(response %in% use_response) %>% 
  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 <- unique(df2$myvar)

Analysis

Calculate "ROM" - the log transformed ratio of means (Hedges et al., 1999; Lajeunesse, 2011) for each observation pair (ambient and elevated).

df3 <- df2 %>%
  
  ## keep only essential variables and drop rows containing missing values for essential variables
  select(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) )


# Aggregate all measurements (multiple years, sampling dates and plots) by experiment (and response variable - although here only one) for meta analysis.

df4 <- df3 %>% 
  
  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(
    df3 %>%
      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) ) %>%
  
  ## Merge with nutnet dataset
  full_join(nutnet.total)
## `summarise()` has grouped output by 'exp', 'myvar'. You can override using the
## `.groups` argument.
## Joining, by = c("exp", "myvar", "logr", "logr_var", "treatment", "n_c", "n_t",
## "logr_se")

Add root mass fraction to var list (for plots)

use_vars2 <- c(use_vars, "rmf")

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("/Users/eaperkowski/git/lt_cn_review/R/analyse_meta.R")
out  <- purrr::map(as.list(use_vars2), ~analyse_meta(df4 %>% rename(var = myvar), 
                                                    nam_target = .))
names(out) <- use_vars2
df_box <- purrr::map_dfr(out, "df_box")
saveRDS(df_box, file = "/Users/eaperkowski/Desktop/df_box_mesi_nfert.rds")

Final data size

Number of data points (plot-level measurements) per variable:

df4 %>% 
  group_by(myvar) %>% 
  summarise(n_plots = sum(n_c), 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 6048 266
anpp NA 50
bgb 1665 138
n_inorg 1043 52
n_uptake 102 4
rmf 91 29
root_production 374 21
root_shoot_ratio 240 54

Number of data points (plot-level measurements) per experiment:

df4 %>% 
  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
-43.52_172.58_f 5 1
34.2_106.43_f 10 1
35.9_-79.03_f 3 1
37.25_-121.75_forb_fn 3 1
37.25_-121.75_forb_fnp 3 1
37.25_-121.75_forb_fp 3 1
37.25_-121.75_grass_fn 3 1
37.25_-121.75_grass_fnp 3 1
37.25_-121.75_grass_fp 3 1
38.53_-76.33_f 50 2
38.883_-107.033_f 24 1
40.4_-86.9_f 9 1
41.77_111.88_f 24 1
44.1_-96.41_f 8 1
47.26_8.31_f 4 1
47.26_8.31_f2 8 1
47.26_8.31_f3 8 1
68.38_-104.54_f 4 1
69.05_20.88_f 21 1
abag_degraded_f 36 1
abag_degraded_f2 36 1
abag_degraded_f3 36 1
abag_degraded_f4 36 1
abag_degraded_f5 36 1
abag_f 38 2
abag_f2 83 4
abag_f3 83 4
abag_f4 83 4
abag_f5 83 4
abag_mature_f 72 2
abag_mature_f2 72 2
abag_mature_f3 72 2
abag_mature_f4 72 2
abag_mature_f5 72 2
alaska_f 12 2
alp_weissenstein2_f 130 3
antu_f 12 1
antu_f2 12 1
arizona_f 16 1
bcnm_fk 20 1
bcnm_fn 20 1
bcnm_fnk 20 1
bcnm_fnp 20 1
bcnm_fnpk 20 1
bcnm_fp 20 1
bcnm_fpk 20 1
bennekom_drained_fn 25 1
bennekom_drained_fnp 5 1
bennekom_drained_fp 5 1
bennekom_undrained_fn 25 1
bennekom_undrained_fnp 5 1
bennekom_undrained_fp 5 1
biocon_f 238 6
bleckstugan_f 1 1
bonanza_creek_2005_f 108 1
bordeaux_f 24 1
bordeaux_f2 24 1
changbai_mountain_f 6 1
changbai_mountain_f2 18 2
changshan_f 200 2
chaux-des-breuleux_f 25 3
chiriqui_f 8 1
chubut_f 91 1
climaite_f 6 1
climaite_f2 6 1
cuiliugou_f 240 2
damxung_2013_f 24 2
damxung_2014_f 48 1
damxung_2014_f2 48 1
damxung_2014_f3 48 1
damxung_2015_f 476 1
damxung_b_f 46 2
damxung_b_f2 39 2
damxung_b_f3 20 2
damxung_f 120 2
daqinggou_fn 36 2
daqinggou_fnp 12 2
daqinggou_fp 12 2
dbr2_f 12 2
dbr2_f2 12 2
dbr2_f3 12 2
dbr2_f4 12 2
deqing_b_f 15 1
deqing_b_f3 27 1
deqing_d_f 21 1
deqing_d_f2 24 1
deqing_d_f3 21 1
deqing_f 5 1
deqing_f2 5 1
dukeface_f 4 1
duolun_2010b_f 55 5
duolun1_f 88 4
duolun12_f 10 1
duolun15_fn 4 1
duolun15_fn2 4 1
duolun15_fn2p 4 1
duolun15_fnp 4 1
duolun2_f 24 1
duolun2_f2 64 1
duolun2_f3 24 1
duolun2_f4 64 1
duolun2_f5 64 1
duolun2_f6 69 2
duolun2_f7 24 1
duolun3_f 185 1
duolun3_f2 175 1
duolun3_f3 185 1
duolun3_f4 185 1
duolun3_f5 185 1
duolun4_f 32 2
duolun7_f 161 4
duolun9_f 132 5
duolun9_w 6 1
eastern_japan_f 10 1
euroface_pa_f 27 3
euroface_pe_f 27 3
euroface_pn_f 27 3
euroface_pooled_f 54 2
flakaliden_f 18 2
flakaliden_f2 4 1
fruebuel2_f 190 3
gaoyao_f 4 1
gaoyao_f2 4 1
gaoyao_f3 8 1
haeggsjoeliden_f 1 1
hamr_f 8 3
hamr_f2 8 3
harbin_f 18 3
hawaii_fn 10 1
hf_f_mh_f 3 1
hf_f_mh_f2 3 1
hf_f_pr_f 3 1
hf_f_pr_f2 3 1
horqin_zuoyi_houqi_f 12 1
horqin_zuoyi_houqi_f2 12 1
horqin_zuoyi_houqi_f3 12 1
horqin_zuoyi_houqi_f4 12 1
horqin_zuoyi_houqi_f5 12 1
horsham_face10bc_f 4 1
horsham_face10bk_f 4 1
imgers_hg_2005_f 96 3
imgers_hg_2005_f2 96 3
imgers_hg_2008_f 32 3
imgers_hg_d_2008_f 2 1
imgers_lg_d_2008_f 2 1
imgers_mg_2008_f 36 3
imgers_mg_d_2007_f 2 1
imgers_ng_2006_f 6 1
imgers_ng_2006_f2 6 1
imgers_ng_2006_f3 6 1
imgers_ng_2006_f4 6 1
imgers_ng_2006_f5 6 1
imgers_ng_2007b_f 40 3
imgers_ng_d_2007_f 2 1
imgers_pooled_2005_f 16 2
imgers_pooled_2006_f 8 1
imgers_pooled_2006_f2 8 1
imgers_pooled_2008c_f 4 1
irvine_ranch_1_f 45 5
itatinga_f 15 1
itatinga_f2 15 1
jilin_f 16 2
jingtai_f 56 2
jingtai_f2 64 2
jingtai_f3 64 2
jrbp_face_f 178 5
kemijaervi_f 1 1
kisa_grahamiana_f 20 1
liaoning_f 6 1
liudaogou_f 72 2
liujiang_b_f 3 1
liujiang_b_f2 3 1
liujiang_b_f3 3 1
luangcheng_f 9 1
luneburg_field_2006_fn 10 1
luneburg_field_2006_fnp 10 1
luneburg_field_2006_fp 10 1
luneburg_field_2008_fn 10 1
luneburg_field_2008_fnp 10 1
luneburg_field_2008_fp 10 1
maoershan_fraxinus_f 3 1
maoershan_larix_f 33 3
maoxian_f 36 3
maoxian_pc_f 5 1
maoxian_pp_f 5 1
massachusetts_f 80 1
massachusetts_f2 80 1
mbs_pt2_f 10 1
menyuan_b_f 12 1
menyuan_b_f2 12 1
menyuan_f 14 2
menyuan_f2 14 2
menyuan_f3 14 2
menyuan_f4 14 2
menyuan_f5 14 2
menyuan_f6 14 2
menyuan_f62 14 2
menyuan_f63 14 2
menyuan_f64 19 2
michigan_c_f 3 1
michigan_e_f 36 1
michigan_pooled_f 3 1
michigan-kbs_grassland_f 120 1
michigana_f 15 2
michiganb_f 3 1
michiganc_f 18 3
niwot_ridge2_dm_f2np 5 1
niwot_ridge2_dm_fn 35 3
niwot_ridge2_dm_fnp 10 2
niwot_ridge2_dm_fp 10 2
niwot_ridge2_wm_fn 25 3
niwot_ridge2_wm_fnp 15 2
niwot_ridge2_wm_fp 10 2
norrliden_f 3 1
norrliden_f2 3 1
norrliden_f3 3 1
norrliden_f4 3 1
norrliden_f5 3 1
norrliden_f6 3 1
nutnet_bldr.us 8 4
nutnet_bnch.us 12 4
nutnet_bogong.au 12 4
nutnet_burrawan.au 12 4
nutnet_cbgb.us 24 4
nutnet_cdcr.us 12 4
nutnet_cdpt.us 12 4
nutnet_cowi.ca 12 4
nutnet_elliot.us 12 4
nutnet_frue.ch 12 4
nutnet_gilb.za 12 4
nutnet_hall.us 12 4
nutnet_hart.us 12 4
nutnet_konz.us 12 4
nutnet_lancaster.uk 8 4
nutnet_look.us 12 4
nutnet_mtca.au 12 4
nutnet_sage.us 12 4
nutnet_saline.us 12 4
nutnet_sgs.us 12 4
nutnet_shps.us 12 4
nutnet_sier.us 12 4
nutnet_smith.us 12 4
nutnet_spin.us 12 4
nutnet_summ.za 12 4
nutnet_trel.us 12 4
nutnet_ukul.za 24 4
nutnet_unc.us 12 4
nutnet_valm.ch 12 4
ornl_face_liqui2_f 19 3
ottawa_fn 3 1
ottawa_fn2 3 1
ottawa_fn3 3 1
pepeekeo_f NA 3
puerto_b_f 3 1
reckenholz_f 24 1
riceface_nianyufarm_2002_f 3 1
riceface_nianyufarm_2003_f 3 1
riceface_zhongcun_2008_f 9 1
riceface_zhongcun_hybrid_2006_f 18 1
rosinedal_f 26 3
rosinedal_f2 20 1
salmisuo_f 25 3
sanjiang_mire_maize_fn 6 2
sanjiang_mire_maize_fn2 6 2
sanjiang_mire_nfert2_f 3 1
sanjiang_mire_nfert2_f2 6 1
sanjiang_mire_nfert2_f3 6 1
sanjiang_mire_pfert_fp 9 1
sanjiang_mire_pfert_fp2 9 1
sanjiang_mire_pfert_fp3 9 1
santa_rosa_f 24 1
santa_rosa_f2 24 1
santa_rosa_f3 24 1
sarobetsu_f 27 1
sarobetsu_f2 27 1
sarobetsu_f3 27 1
sarobetsu_f4 27 1
sca_f 24 2
scbg_f 6 2
schiermonnikoog_old_fn 48 1
schiermonnikoog_old_fn2 48 1
schiermonnikoog_old_fn2p 12 1
schiermonnikoog_old_fn2p2 12 1
schiermonnikoog_old_fnp 12 1
schiermonnikoog_old_fnp2 12 1
schiermonnikoog_old_fp 12 1
schiermonnikoog_old_fp2 12 1
schiermonnikoog_young_fn 48 1
schiermonnikoog_young_fn2 48 1
schiermonnikoog_young_fn2p 12 1
schiermonnikoog_young_fn2p2 12 1
schiermonnikoog_young_fnp 12 1
schiermonnikoog_young_fnp2 12 1
schiermonnikoog_young_fp 12 1
schiermonnikoog_young_fp2 12 1
shaaxi_330_bungeana_fn 3 1
shaaxi_330_bungeana_fn3 6 1
shaaxi_330_bungeana_fnp 3 1
shaaxi_330_bungeana_fp 3 1
shaaxi_330_davarica_fn 3 1
shaaxi_330_davarica_fn3 6 1
shaaxi_330_davarica_fnp 3 1
shaaxi_330_davarica_fp 3 1
shaaxi_330_fn 6 1
shaaxi_330_fn2 3 1
shaaxi_330_fn3 9 1
shaaxi_330_fnp 3 1
shaaxi_330_fp 3 1
shaaxi_6_bungeana_fn 3 1
shaaxi_6_bungeana_fn2 3 1
shaaxi_6_bungeana_fn3 3 1
shaaxi_6_bungeana_fn3p 3 1
shaaxi_6_bungeana_fp 3 1
shaaxi_6_davarica_fn 3 1
shaaxi_6_davarica_fn2 3 1
shaaxi_6_davarica_fn3 3 1
shaaxi_6_davarica_fn3p 3 1
shaaxi_6_davarica_fp 3 1
shaaxi_6_fn 6 1
shaaxi_6_fn2 6 1
shaaxi_6_fn3 6 1
shaaxi_6_fn3p 3 1
shaaxi_6_fp 3 1
shaxian_f 9 2
shaxian_f2 9 2
shaxian_f3 9 2
silwood_f 80 1
sjaellirimsheden_f 1 1
skogaby_f2 4 1
slattatjakka_f 80 1
sodankylae_f 1 1
songen_f 96 2
strasan_f 2 1
strasan_f2 2 1
strasan_f3 2 1
swissface_lolium2_f 30 4
swissface_loliumtrifolium2_f 18 2
swissface_trifolium2_f 30 4
sydney_f 48 1
sydney_f2 48 1
teberda_f 32 1
tiao_xi_nan_lu_f 16 1
tiao_xi_nan_lu_f2 16 1
tiao_xi_nan_lu_f3 16 1
toolik_nonacidic_f 32 1
trebon_basin_biosphere_reserve_f 8 3
trebon_basin_biosphere_reserve_f2 8 3
tu_2011_f 6 1
tu_2011_f2 6 1
tu_2011_f3 6 1
tu_b_f 3 1
tu_b_f2 3 1
tu_b_f3 3 1
ulan_f 216 1
ulan_f2 216 1
ulan_f3 216 1
ulan_f4 216 1
ulan_f5 216 1
vikki_f 30 1
vikki_f2 30 1
vikki_f3 30 1
winnter_f 74 4
xiang_dao_f 12 1
xiaojin_b_f 15 1
xiaojin_b_f2 15 1
xiaojin_b_f3 12 1
xun_f 8 1
xun_f2 8 1
xun_f3 8 1
xun_f4 8 1
yakeshi_f 6 1
yakeshi_f2 6 1
zengcheng_f 8 1

Plots

df4 %>%
  drop_na(myvar) %>% 
  
  ## give it a nice order (for plotting)
  mutate(myvar = factor(myvar, levels = rev(c("anpp", "agb", "root_production", "bgb", "root_shoot_ratio", "rmf", "n_uptake", "n_inorg")))) %>% 
  ggplot( aes(x = myvar, y = logr)) +
  geom_jitter(color = rgb(0,0,0,0.5), 
                aes( size = 1/logr_se ), 
                position = position_jitter(w = 0.2, h = 0) 
  ) +
  geom_crossbar( data = df_box %>% drop_na(var), 
                 aes(x = var, y = middle, ymin = ymin, ymax = ymax),
                 fill = "tomato", color = "black", alpha = 0.6, width = 0.5 ) +
  geom_hline( yintercept = 0.0, size = 0.5 ) +
  scale_y_continuous(limits = c(-2, 2), breaks = seq(-2, 2, 1)) +
  labs(x = "Variable", y = "Log Response Ratio", size = expression(paste("Error"^{-1}))) +
  coord_flip() +
  labs(title = "MESI + NutNet data", subtitle = "Response to N-fertilisation")

Couplings

Make a wide table with ANPP and BNPP:

df4_agb_bnpp <- df4 %>% 
  filter(myvar %in% c("agb", "bgb")) %>% 
  select(exp, myvar, logr) %>%
  tidyr::spread( myvar, logr )

## add standard error 
df4_agb_bnpp <- df4 %>% 
  filter(myvar %in% c("agb", "bgb")) %>% 
  select(exp, myvar, logr_se) %>%
  tidyr::spread( myvar, logr_se ) %>% 
  rename(se_agb = agb, se_bgb = bgb) %>%
  right_join(df4_agb_bnpp, by = "exp") %>% 
  mutate(se = se_agb * se_bgb)

df4_agb_bnpp %>% 
  ggplot(aes(x = agb, y = bgb, label = exp)) +
  geom_point(aes(size =  1/se), color = "tomato") +
  xlim(0, 0.6) + ylim(0, 0.6) +
  geom_abline(linetype = "dotted") + 
  ggrepel::geom_text_repel(size = 3, point.padding = 0.5, segment.alpha = 0, color = "grey50") +
  labs(x = "Log RR of AGB", y = "Log RR of BGB", size=expression(paste("Error"^{-1}))) +
  theme_classic()