N2O experiments cross-check

Read data

df_temp <- read_csv("~/data/n2o_Yunke/final_obs_dataset/obs_warming_dataset.csv") |> 
  
  # give it a unique ID - missing in the dataset. Also missing: an experiment name.
  mutate(id = 1:n()) |> 
  
  mutate(logrr = log(n2o_elv / n2o_amb)) |>
  mutate(logrr_norm = logrr / dT) |> 

  # log-transform - XXX not necessary, the errors should be normally distributed, not necessarily the predictors.
  mutate(PPFD_total_a = log(PPFD_total),
         vpd_a = log(vpd),
         orgc_a = log(ORGC),
         ndep_a = log(ndep)
         )

## Rows: 83 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (5): original_ref, location, pft, species, other_treatment
## dbl (18): lon, lat, n2o_amb, n2o_elv, dT, logr, Nfer_kgha, days_of_duration,...
## 
## ℹ 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.

Explore data

Normalising the log response reduced the variation in the log response.

df_temp |> 
  pivot_longer(cols = c(logrr, logrr_norm), names_to = "var", values_to = "logrr") |> 
  ggplot(aes(x = var, y = logrr)) +
  geom_boxplot()

Issue:

• Most data is from a crossed experimental treatment, with simultaneous treatments in CO2, N inputs, with/without plants, rain exclusion, etc. How do we go about that? It’s not very meaningful to regress the response to variables when the same variable may be experimentally manipulated (often the case for N inputs in experiments vs. N deposition as environmental predictor, precipitation manipulation vs. precipitation as environmental predictor, …).

df_temp |> 
  select(id, Nfer_kgha, other_treatment) |> 
  knitr::kable()

enrichment, Armidale C | | 81| 0.0|CO 2 enrichment, Armidale R | | 82| 0.0|Armidale R | | 83| 0.0|NA |

We could be conservative, using only data from experiments with no other simultaneous treatment (assuming that this is the case when other_treatment == NA). This yields 36 experiments (out of 83 originally).

df_temp_sub <- df_temp |> 
  filter(is.na(other_treatment))

df_temp |> dim()

## [1] 83 30

df_temp_sub |> dim()

## [1] 36 30

This further reduces the variation.

ggplot() +
  geom_boxplot(aes(x = 1, y = logrr), data = df_temp) +
  geom_boxplot(aes(x = 2, y = logrr_norm), data = df_temp) +
  geom_boxplot(aes(x = 3, y = logrr_norm), data = df_temp_sub)

Issue:

• Some data seem to be from the same experiment, but from multiple parallel treatments. Yet, environmental factors considered for modelling are extracted from global files (I guess) and are the same for all treatments within a given experiment. Hence, a (probably large) part of the uncertainty remains unexplained by these environmental factors.

Stepwise regression

# following example: https://www.statology.org/stepwise-regression-r/ 
# subset, using un-transformed predictors
df <- df_temp |> 
  select(logrr_norm, min_fapar, max_fapar, PPFD_total, vpd, ORGC, ndep) |>   # xxx Nfer doesn't exist in that data
  drop_na()

# define intercept-only model
intercept_only <- lm(logrr_norm ~ 1, data = df)

#define model with all predictors
all <- lm(logrr_norm ~ ., data = df)

#perform forward stepwise regression
forward <- step(intercept_only, direction = 'forward', scope = formula(all), trace = 0)

#view results of forward stepwise regression
forward$anova

##   Step Df Deviance Resid. Df Resid. Dev       AIC
## 1      NA       NA        64   20.01275 -74.57116

#view final model
forward$coefficients

## (Intercept) 
## -0.04563139

# subset, using un-transformed predictors
df <- df_temp |> 
  select(logrr_norm, min_fapar, max_fapar, PPFD_total_a, vpd_a, orgc_a, ndep_a) |>  # xxx Nfer doesn't exist in that data
  drop_na()

# define intercept-only model
intercept_only <- lm(logrr_norm ~ 1, data = df)

#define model with all predictors
all <- lm(logrr_norm ~ ., data = df)

#perform forward stepwise regression
forward <- step(intercept_only, direction = 'forward', scope = formula(all), trace = 0)

#view results of forward stepwise regression
forward$anova

##       Step Df  Deviance Resid. Df Resid. Dev       AIC
## 1          NA        NA        64   20.01275 -74.57116
## 2 + orgc_a -1 0.8807562        63   19.13199 -75.49665

#view final model
forward$coefficients

## (Intercept)      orgc_a 
##  -0.2875957   0.1161154

Predict

For spatial upscaling, use a map of SOC and use predict(). Here, the logarithm of the organic matter content (variable ORG is, say, 1.5)

newdata <- data.frame(orgc_a = 1.5)

logrr_predicted <- predict(forward, newdata = newdata)

In this case, the log-response ratio is -0.11. Convert this to a response ratio:

exp(logrr_predicted)

##         1 
## 0.8927733

Yes, that’s a decline in the N2O emissions for this particular choice SOC content.

For spatial uscaling, you apply the predict() function as above to the SOC content for each pixel.