Experiments data analysis with Van Groenigen et al., 2014 data

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

Load data files.

df <- read_csv("~/data/gcme/data_received_190325/soilc_decomp_vangroenigen14.csv")

## New names:
## Rows: 53 Columns: 23
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (4): exp_nam_vangroenigen, prev_name, treatment, Ecosystem dbl (19): I_ambient,
## k_ambient, C_ambient, tau...7, I_elevated, k_elevated, ...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `tau` -> `tau...7`
## • `tau` -> `tau...11`

Make data (more) conform with MESI.

df1 <- df %>% 
  ## Back-calculate SD of k from column wVk = 1/(SD_ambient^2/k_ambient^2 + SD_elevated^2/k_elevated^2)
  ## assuming that ambient_Sd = elevated_Sd
  rowwise() %>% 
  rename( rep_c = n_ambient, rep_t = n_elevated, x_c = k_ambient, x_t = k_elevated,
           exp = exp_nam_vangroenigen ) %>% 
  mutate( sd_c = x_c * x_t / sqrt(wVk * (x_c^2 + x_t^2)) ) %>%
  mutate( sd_t = sd_c,
          response = "kdecay_soil" )

Data overview

Issues

• …

Data subsetting

Let’s use just data from CO2-only experiments.

df2 <- df1 %>% 
  
  ## Use CO2-only treatment
  filter(treatment == "c")

Analysis

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

df3 <- df2 %>%
  mutate(id = 1:n()) %>% 
  filter(sd_c > 0, sd_t > 0) %>% 
  
  ## keep only essential variables and drop rows containing missing values for essential variables
  select(id, exp, response, x_t, x_c, sd_t, sd_c, rep_t, rep_c) %>% 
  drop_na() %>% 
  
  ## 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 (species, 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)
  mutate( id = exp ) %>% 
  
  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() %>% 
  
  ## rename again
  mutate(response = "kdecay_soil") %>% 
  select(exp = id, response, logr = es, logr_var = var) %>%

  ## add number of observations (sum of plots and repeated samplings)
  left_join(
    df3 %>%
      group_by(exp, response) %>%
      summarise(n_c = sum(rep_c), n_t = sum(rep_t)),
    by = c("exp", "response")
  ) %>% 
  
  ## 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'. 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  <- analyse_meta(df4 %>% rename(var = response), nam_target = "kdecay_soil")
df_box <- out$df_box

Visualisations

Plot dots and my box.

df4 %>%
  ggplot( aes(x = response, y = logr)) +
  geom_jitter( color = rgb(0,0,0,0.3), 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, linetype = "dotted") +
  labs(x = "Variable", y = "Log Response Ratio", size = expression(paste("Error"^{-1}))) +
  coord_flip() +
  ylim(-1, 1)  +
  labs(title = "Van Groenigen et al. 2014 data", subtitle = "Response to CO2")

Findings:

• Acceleration of soil decomposition rates (this refers to the decay rate parameter, not an absolute flux!)