Experiments data analysis with MESI

library(tidyverse)   # for working with ease

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.0      ✔ stringr 1.4.1 
## ✔ readr   2.1.2      ✔ forcats 0.5.1 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(metafor)     # for meta analysis

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loading required package: metadat
## 
## Loading the 'metafor' package (version 3.8-1). For an
## introduction to the package please type: help(metafor)

library(MAd)         # for meta analysis

Read data and libraries

Load and interpret MESI data files.

## second line provides units
col_names <- names(read_csv("~/mesi-db/data/mesi_main.csv", n_max = 0))
units <- names(read_csv("~/mesi-db/data/mesi_main.csv", n_max = 0, skip = 1))
df_units <- tibble(column = col_names,
                   units = units
                   )
df <- read_csv("~/mesi-db/data/mesi_main.csv", col_names = col_names, skip = 2)

## Warning: One or more parsing issues, see `problems()` for details

CO2 experiments

Select experiments from which we will use data.

df0 <- df %>%
  
  ## CO2-only manipulation (no other factors manipulated, no crossed manipulation)
  filter(treatment == "c") %>% 
  
  ## only experiments conducted in the field
  filter(experiment_type == "field") %>% 
  
  ## use only free air CO2 enrichment experiments and open-top chambers
  filter(fumigation_type %in% c("face", "otc"))

Data overview

Issues

• Some experiments don’t have a name, referred to by coordinates. Can they be identified to get an experiments name? (-> STUDENT)

Data subsetting

Let’s use just data from CO2-only experiments that lasted at least three years.

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

Issues:

• No data available for Narea (-> KEVIN)

## create own variable names (myvar), grouping available variables
use_response <- c("asat", "vcmax", "jmax", "gpp", "gs", "leaf_n", "leaf_cn", "anpp", "agb", "leaf_biomass", "lai", "lai_max", "leaf_area", "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) %>% 

  ## variables re-grouped by myself
  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("lai", "lai_max"), "lai", 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

Test analysis and plot of ANPP data. Calculate the response ratio of ANPP (mean and variance) for each experiment. To get that, we first need to calcuate the means and standard deviation for the ambient and elevated levels, pooling multiple measurements (years, sampling dates), each given with mean \(\mu_i\), number \(N_i\) (replicates/plots), and standard deviation \(\sigma_i\) or standard error. For the function metafor::escalc(), we need standard deviations (\(SD\)). Calculate them for those rows where only standard errors \(SE\) are given as: \[ SD = SE \sqrt{N} \]

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) )

## Warning: Some 'yi' and/or 'vi' values equal to +-Inf. Recoded to NAs.

# ## data reduction here!
# tmp3 <- purrr::map(as.list(use_vars), ~agg_by_var(., df3))
# names(tmp3) <- use_vars
# tmp3$root_shoot_ratio

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) )

## `summarise()` has grouped output by 'exp', 'myvar'. You can override using the
## `.groups` argument.

# tmp <- purrr::map(as.list(use_vars), ~agg_by_var(., df4 %>% rename(rep_c = n_c)))
# names(tmp) <- use_vars
# tmp$root_shoot_ratio

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("~/lt_cn_review/R/analyse_meta.R")
out  <- purrr::map(as.list(use_vars), ~analyse_meta(df4 %>% rename(var = myvar), nam_target = .))
names(out) <- use_vars
df_box <- purrr::map_dfr(out, "df_box")

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 plots"="n_plots", "N experiments"="n_exp") %>% 
  knitr::kable()

## Warning: `rename_()` was deprecated in dplyr 0.7.0.
## Please use `rename()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Variable	N plots	N experiments
agb	1715	61
anpp	336	22
asat	553	19
bgb	1237	50
gpp	128	7
gs	NA	17
jmax	205	7
lai	NA	17
leaf_area	78	6
leaf_biomass	426	16
leaf_cn	105	5
leaf_n	1801	17
n_inorg	NA	5
n_uptake	52	5
root_production	361	11
root_shoot_ratio	172	16
vcmax	249	9

Available experiments for root:shoot ratio

ivar <- "root_shoot_ratio"
df4 %>%
  filter(myvar == ivar) %>%
  select(exp, n_c) %>%
  knitr::kable()

exp	n_c
duolun7_c	12
climaite_c	6
46.77_9.87_c	3
swissface_forest_c	24
jrbp_face_c	8
euroface_pooled_c	3
swissface_lolium_c	6
sca_c	3
eucface_c	3
euroface_pa_c	6
euroface_pn_c	6
euroface_pe_c	6
maricopaface_cotton91_c	68
ornl_face_liqui_c	2
swissface_chalk_c	12
maricopaface_wheat94_c	4

Available experiments for inorganic N availability

ivar <- "n_inorg"
df4 %>%
  filter(myvar == ivar) %>%
  select(exp, n_c) %>%
  knitr::kable()

exp	n_c
climaite_c	54
ornl_face_liqui2_c	NA
dukeface_c	32
euroface_pooled_c	18
riceface_nianyufarm_triticum_2012_c	30

Available experiments for N uptake

ivar <- "n_uptake"
df4 %>%
  filter(myvar == ivar) %>%
  select(exp, n_c) %>%
  knitr::kable()

exp	n_c
climaite_c	36
ornl_face_liqui_c	2
euroface_pooled_c	6
aspenface_c	6
ornl_face_liqui2_c	2

Visualisations

Plot dots and my box.

df4 %>%
  ## give it a nice order (for plotting)
  mutate(myvar = factor(myvar, levels = rev(c("asat", "gpp", "vcmax", "jmax", "gs", "leaf_n", "leaf_cn", "lai", "leaf_area", "leaf_biomass", "anpp", "agb", "root_production", "bgb", "root_shoot_ratio", "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, 
                 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 ) +
  labs(x = "Variable", y = "Log Response Ratio", size = expression(paste("Error"^{-1}))) +
  coord_flip() +
  # ylim(-1, 1) +
  labs(title = "MESI data", subtitle = "Response to eCO2")

## Warning: Removed 6 rows containing missing values (geom_point).

Findings:

• Strong positive response of leaf-level assimilation under light saturation (asat) and GPP
• Small (non-significant) reduction in Vcmax, no change in Jmax
• Clear decline in leaf N (I think it’s Nmass, but to be confirmed -> KEVIN)
• Large range in plant-level (correct: plant-level? -> KEVIN) total leaf area
• Positive response in leaf biomass, even stronger than in LAI, indicating increase in (ecosystem-mean) LAI (-> KEVIN: to confirm: leaf_biomass and LAI are both ecosystem-level, while leaf_area is plant-level?)
• Strong positive response in root production - strongest response among all variables and stronger than response in ANPP (even considering negative outlier)
• Stronger response in belowground biomass (bgb) than in aboveground biomass (agb) (but outlier in agb), Consistent with the positive response in root:shoot ratio (but large variation between experiments).
• Positive response in N uptake, but large variation between experiments.
• Decline in mineral N availability (but not significant)
• Strong (in experiment extraordinarily strong) reduction in stomatal conductance (gs).

Standard forest plots

for (ivar in use_vars){
  print(ivar)
  try(metafor::forest(out[[ivar]]$modl))
}

## [1] "bgb"

## [1] "leaf_n"

## [1] "gs"

## [1] "asat"

## [1] "agb"

## [1] "gpp"

## [1] "root_shoot_ratio"

## [1] "n_inorg"

## [1] "leaf_biomass"

## [1] "root_production"

## [1] "n_uptake"

## [1] "anpp"

## [1] "vcmax"

## [1] "leaf_cn"

## [1] "jmax"

## [1] "leaf_area"

## [1] "lai"

# purrr::map(as.list(use_vars), ~try(metafor::forest(out[[.]]$modl)))

Issues:

• Find out in which experiments, belowground biomass (bgb) was actually measured and not estimated based on aboveground biomass (or other aboveground observations) and scaling assumptions (-> KEVIN, STUDENT).
• Negative outliers for agb and root_production and positive outliers for anpp (-> KEVIN, remove them?)
• Euroface and Aspenface are provided as separate experiments. And several (versions of?) “Aspenface” are treated as separate experiments. Should they be? How do they differ? (-> KEVIN)
• Very little data on several variables (N uptake, root:shoot ratio, inorganic N availability, BNPP, …). Wasn’t there more in GCME (-> Beni)? Can we get more from published experiments in general (-> STUDENT)? Show data that was available in GCME for these variables:

## file generated by cnreview/dataanalysis.Rmd
df_gcme <- readRDS("~/lt_cn_review/data/df_c_agg.rds")
write_csv(df_gcme, file = "~/lt_cn_review/data/df_co2_gcme_usedbybeni.csv")

N uptake

## GCME
df_gcme %>% 
  filter(my_varnam == "my_nup") %>% 
  select(exp_nam, my_varnam, n_plots)

## # A tibble: 11 × 3
##    exp_nam                             my_varnam n_plots
##    <chr>                               <chr>       <dbl>
##  1 Brandbjerg                          my_nup         24
##  2 RiceFACE_Japan_A_2003_39,38_140,57  my_nup         36
##  3 RiceFACE_Japan_Ki_2004_39,38_140,57 my_nup          8
##  4 POPFACE_pa                          my_nup          6
##  5 POPFACE_pe                          my_nup          6
##  6 POPFACE_pn                          my_nup          6
##  7 Durham_DukeFACE                     my_nup          6
##  8 ORNERP_liqui                        my_nup          2
##  9 Rhine-aspenFACE                     my_nup          6
## 10 ORNERP_liqui2                       my_nup          2
## 11 SwissFACE_lolium2                   my_nup          6

## MESI
df4 %>% 
  filter(myvar == "n_uptake") %>% 
  select(exp, myvar, n_c)

## # A tibble: 5 × 3
##   exp                myvar      n_c
##   <chr>              <chr>    <dbl>
## 1 climaite_c         n_uptake    36
## 2 ornl_face_liqui_c  n_uptake     2
## 3 euroface_pooled_c  n_uptake     6
## 4 aspenface_c        n_uptake     6
## 5 ornl_face_liqui2_c n_uptake     2

Root:shoot ratio

## GCME
df_gcme %>% 
  filter(my_varnam == "my_rootshootratio") %>% 
  select(exp_nam, my_varnam, n_plots)

## # A tibble: 18 × 3
##    exp_nam             my_varnam         n_plots
##    <chr>               <chr>               <dbl>
##  1 DRI                 my_rootshootratio       9
##  2 JRBP_OTCsand        my_rootshootratio       4
##  3 JRBP_OTCser         my_rootshootratio       4
##  4 MBS_pe              my_rootshootratio       5
##  5 MBS_pt              my_rootshootratio       4
##  6 ORNERP_liriodendron my_rootshootratio       2
##  7 ORNERP_quercus      my_rootshootratio       2
##  8 PHACE               my_rootshootratio      15
##  9 POPFACE_pa          my_rootshootratio       9
## 10 POPFACE_pe          my_rootshootratio       9
## 11 POPFACE_pn          my_rootshootratio       9
## 12 SCBG                my_rootshootratio       3
## 13 SwissFACE_lolium2   my_rootshootratio       6
## 14 UA_OTC              my_rootshootratio       8
## 15 USDA_citrus         my_rootshootratio       2
## 16 EUROFACE4_pa        my_rootshootratio       3
## 17 EUROFACE4_pe        my_rootshootratio       3
## 18 EUROFACE4_pn        my_rootshootratio       3

## MESI
df4 %>% 
  filter(myvar == "n_uptake") %>% 
  select(exp, myvar, n_c)

## # A tibble: 5 × 3
##   exp                myvar      n_c
##   <chr>              <chr>    <dbl>
## 1 climaite_c         n_uptake    36
## 2 ornl_face_liqui_c  n_uptake     2
## 3 euroface_pooled_c  n_uptake     6
## 4 aspenface_c        n_uptake     6
## 5 ornl_face_liqui2_c n_uptake     2

Inorganic N

## GCME
df_gcme %>% 
  filter(my_varnam == "my_navl") %>% 
  select(exp_nam, my_varnam, n_plots)

## # A tibble: 20 × 3
##    exp_nam                      my_varnam n_plots
##    <chr>                        <chr>       <dbl>
##  1 TasFACE                      my_navl       189
##  2 MI                           my_navl         8
##  3 ORNERP_liqui2                my_navl        38
##  4 Headley_fe                   my_navl         4
##  5 Headley_ps                   my_navl         4
##  6 Headley_qp                   my_navl         6
##  7 POPFACE_pa                   my_navl         6
##  8 POPFACE_pe                   my_navl         6
##  9 POPFACE_pn                   my_navl         6
## 10 DUKE                         my_navl        24
## 11 SwissFACE_lolium2            my_navl       126
## 12 EucFACE                      my_navl       228
## 13 POPFACE_pooled               my_navl        36
## 14 EUROFACE7_pooled             my_navl       144
## 15 Brandbjerg                   my_navl       252
## 16 BioCON                       my_navl        69
## 17 GiFACE                       my_navl       441
## 18 DUKE2                        my_navl        56
## 19 RiceFACE_China_32N_120E_Tr_1 my_navl        90
## 20 DRI                          my_navl        30

## MESI
df4 %>% 
  filter(myvar == "n_inorg") %>% 
  select(exp, myvar, n_c)

## # A tibble: 5 × 3
##   exp                                 myvar     n_c
##   <chr>                               <chr>   <dbl>
## 1 climaite_c                          n_inorg    54
## 2 ornl_face_liqui2_c                  n_inorg    NA
## 3 dukeface_c                          n_inorg    32
## 4 euroface_pooled_c                   n_inorg    18
## 5 riceface_nianyufarm_triticum_2012_c n_inorg    30

Data in GCME from:

df_navl <- read_rds("~/lt_cn_review/data/df_navl.rds")
df_navl %>% select(prev_name, Data_type, Fumigation_type, Source_Reference)

## # A tibble: 60 × 4
##    prev_name           Data_type      Fumigation_type Source_Reference        
##    <chr>               <chr>          <chr>           <chr>                   
##  1 EucFACE_c           soil_mineral_N FACE            Hasegawa et al 2016     
##  2 EucFACE_c           soil_mineral_N FACE            Ochoa-Hueso et al 2017  
##  3 EUROFACE7_pooled_c  soil_mineral_N FACE            Moscatelli et al 2005b  
##  4 POPFACE_pa_c        soil_mineral_N FACE            Calfapietera et al 2003a
##  5 POPFACE_pe_c        soil_mineral_N FACE            Calfapietera et al 2003a
##  6 POPFACE_pn_c        soil_mineral_N FACE            Calfapietera et al 2003a
##  7 POPFACE_pooled_c    soil_mineral_N FACE            Moscatelli et al 2005a  
##  8 POPFACE_pooled_c    soil_mineral_N FACE            Moscatelli et al 2005b  
##  9 SwissFACE_lolium2_c soil_mineral_N FACE            Fromin et al 2005       
## 10 TasFACE_c           soil_mineral_N FACE            Hayden et al. 2012      
## # … with 50 more rows

Couplings

Plot effect sizes for two different variables against each other, given that data is available for both variables for a given experiment.

ANPP-root production

## make a wide table with ANPP and BNPP
df4_anpp_bnpp <- df4 %>% 
  filter(myvar %in% c("anpp", "root_production")) %>% 
  select(exp, myvar, logr) %>%
  tidyr::spread( myvar, logr )

## add standard error 
df4_anpp_bnpp <- df4 %>% 
  filter(myvar %in% c("anpp", "root_production")) %>% 
  select(exp, myvar, logr_se) %>%
  tidyr::spread( myvar, logr_se ) %>% 
  rename(se_anpp = anpp, se_root_production = root_production) %>%
  right_join(df4_anpp_bnpp, by = "exp") %>% 
  mutate(se = se_anpp * se_root_production)

df4_anpp_bnpp %>% 
  ggplot(aes(x = anpp, y = root_production, label = exp)) +
  geom_point(aes(size =  1/se), color = "tomato") +
  xlim(0, 0.8) + ylim(0, 0.8) +
  geom_abline(linetype = "dotted") + 
  ggrepel::geom_text_repel(size = 3, point.padding = 0.5, segment.alpha = 0, color = "grey50") +
  labs(x = "Log RR of ANPP", y = "Log RR of root production", size = expression(paste("Error"^{-1}))) +
  theme_classic()

## Warning: Removed 19 rows containing missing values (geom_point).

## Warning: Removed 19 rows containing missing values (geom_text_repel).

AGB - BGB

## 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()

## Warning: Removed 52 rows containing missing values (geom_point).

## Warning: Removed 52 rows containing missing values (geom_text_repel).

Issues:

• SwissFACE Lolium: I seem to remember that there was a stronger stimulation of belowground than aboveground biomass. Why does this here suggest a much higher response in AGB than in BGB? In the SwissFACE experiment, ABG was determined from multiple cuts per year. Within each year, harvested (cut) biomass should be summed up, not averaged! Please check and correct: where sums are relevant, take sum over multiple samplings per year. (-> KEVIN, STUDENT)

N fertilisation experiments

Select experiments from which we will use data.

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

Data subsetting

Same as above

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

## select variables only variables as used for CO2 analysis and not provided by Liang et al.
## create own variable names (myvar), grouping available 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) %>% 

  ## variables re-grouped by myself
  ## variables re-grouped by myself
  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("lai", "lai_max"), "lai", 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) )

## Warning in log(m1i/m2i): NaNs produced

## Warning: Some 'yi' and/or 'vi' values equal to +-Inf. Recoded to NAs.

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) )

## `summarise()` has grouped output by 'exp', 'myvar'. You can override using the
## `.groups` argument.

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("~/lt_cn_review/R/analyse_meta.R")
out  <- purrr::map(as.list(use_vars), ~analyse_meta(df4 %>% rename(var = myvar), nam_target = .))
names(out) <- use_vars
df_box <- purrr::map_dfr(out, "df_box")

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()

Variable	N observations	N experiments
agb	5957	237
anpp	NA	50
bgb	1580	109
n_inorg	1043	52
n_uptake	102	4
root_production	374	21
root_shoot_ratio	149	25

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	184	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
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

Visualisations

Plot dots and my box.

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", "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 ) +
  labs(x = "Variable", y = "Log Response Ratio", size = expression(paste("Error"^{-1}))) +
  coord_flip() +
  labs(title = "MESI data", subtitle = "Response to N-fertilisation")

## Warning: Removed 1 rows containing missing values (geom_point).

Findings:

• Stronger effect on aboveground biomass than on belowground biomass - this comes out very clearly thanks to lots of data. One outlier of a very negative agb response (-> STUDENT?)
• Reduction in root:shoot ratio - opposite to CO2 response and consistent with expectations from a functional balance.
• Unsurprisingly, inorganic soil N increases in some cases very strongly, but large variation between experiments.
• No significant change in root production. Is this reliable? Or just to difficult to measure?

Issues:

• Get more data for root:shoot ratio from N-fertilisation experiments (-> Evan).
• Get more data for N-inorg. data from N-fertilisation experiments (-> Evan).

for (ivar in use_vars){
  print(ivar)
  try(metafor::forest(out[[ivar]]$modl))
}

## [1] "anpp"

## [1] "root_production"

## [1] "agb"

## [1] "bgb"

## [1] "root_shoot_ratio"

## [1] "n_inorg"

## [1] "n_uptake"

# purrr::map(as.list(use_vars), ~try(metafor::forest(out[[.]]$modl)))

Couplings

Plot effect sizes for two different variables against each other, given that data is available for both variables for a given experiment.

ANPP-root production

## make a wide table with ANPP and BNPP
df4_anpp_bnpp <- df4 %>% 
  filter(myvar %in% c("anpp", "root_production")) %>% 
  select(exp, myvar, logr) %>%
  tidyr::spread( myvar, logr )

## add standard error 
df4_anpp_bnpp <- df4 %>% 
  filter(myvar %in% c("anpp", "root_production")) %>% 
  select(exp, myvar, logr_se) %>%
  tidyr::spread( myvar, logr_se ) %>% 
  rename(se_anpp = anpp, se_root_production = root_production) %>%
  right_join(df4_anpp_bnpp, by = "exp") %>% 
  mutate(se = se_anpp * se_root_production)

df4_anpp_bnpp %>% 
  ggplot(aes(x = anpp, y = root_production, label = exp)) +
  geom_point(aes(size =  1/se), color = "tomato") +
  xlim(0, 0.8) + ylim(0, 0.8) +
  geom_abline(linetype = "dotted") + 
  ggrepel::geom_text_repel(size = 3, point.padding = 0.5, segment.alpha = 0, color = "grey50") +
  labs(x = "Log RR of ANPP", y = "Log RR of root production", size = expression(paste("Error"^{-1}))) +
  theme_classic()

## Warning: Removed 54 rows containing missing values (geom_point).

## Warning: Removed 54 rows containing missing values (geom_text_repel).

Issues:

• Not enough data. Drop this analysis.

AGB - BGB

## 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()

## Warning: Removed 282 rows containing missing values (geom_point).

## Warning: Removed 282 rows containing missing values (geom_text_repel).

Issues:

• Very few data on belowground biomass response. No estimates available from other experiments, e.g., from NutNet (-> Evan)