---
title: Statistical analysis
subtitle: Litter decomposition
author: <a href="https://marce10.github.io/">Marcelo Araya-Salas</a> & Andrea Vincent
date: "`r Sys.Date()`"
toc: true
toc-depth: 2
toc-location: left
number-sections: true
highlight-style: pygments
format:
html:
self-contained: true
df-print: kable
code-fold: true
code-tools: true
css: qmd.css
editor_options:
chunk_output_type: console
---
```{r set root directory, echo = FALSE}
# set working directory
knitr::opts_knit$set(root.dir = "..")
```
```{r add link to github repo, echo = FALSE, results='asis'}
# print link to github repo if any
if (file.exists("./.git/config")){
config <- readLines("./.git/config")
url <- grep("url", config, value = TRUE)
url <- gsub("\\turl = |.git$", "", url)
cat("\nSource code and data found at [", url, "](", url, ")", sep = "")
}
```
```{r setup style, echo = FALSE, message = FALSE, warning=FALSE}
# options to customize chunk outputs
knitr::opts_chunk$set(
class.source = "numberLines lineAnchors", # for code line numbers
tidy.opts = list(width.cutoff = 65),
tidy = TRUE,
message = FALSE,
warning = FALSE
)
```
<!-- skyblue box -->
<div class="alert alert-info">
# Purpose {.unnumbered .unlisted}
- Evaluate role of nutrient availability on litter decomposition
</div>
```{r load packages, echo = FALSE, message = FALSE, warning=FALSE, eval = TRUE}
# remotes::install_github("maRce10/sketchy")
remotes::install_local("~/Dropbox/R_package_testing/sketchy/")
# github packages must include user name ("user/package")
pkgs <- c("kableExtra", "knitr", "readxl", "ggplot2", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "viridis", "brms", "ggdist", "brmsish")
# install/ load packages
sketchy::load_packages(pkgs, quite = TRUE, upgrade.deps = TRUE)
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
```
```{r functions and global parameters, echo = TRUE, message = FALSE, warning=FALSE}
# custom function based on gghalves::geom_half_point
geom_half_point <- function(
mapping = NULL,
data = NULL,
side = "l",
range_scale = 0.4,
alpha = 1,
pch = 16,
size = 1.5,
transformation = NULL, # <- absorb but ignore
...
) {
ggplot2::geom_point(
mapping = mapping,
data = data,
alpha = alpha,
shape = pch,
size = size,
position = position_half_jitter(
width = 0.1,
height = 0,
side = side,
range_scale = range_scale
),
...
)
}
position_half_jitter <- function(width = 0.1,
height = 0,
side = "l",
range_scale = 0.4) {
ggplot2::ggproto(
NULL, ggplot2::Position,
compute_panel = function(self, data, params, scales) {
# apply regular jitter first
data$x <- data$x + stats::runif(
nrow(data),
-width / 2,
width / 2
)
# shift to one side
shift <- if (side %in% c("l", "left")) {
-range_scale / 2
} else {
range_scale / 2
}
data$x <- data$x + shift
data
}
)
}
cols <- viridis(10, alpha = 0.7)
fill_color <- viridis(10)[7]
# brms models
chains <- 4
iters <- 4000
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
prior(student_t(3, 0, 20), "sd"))
# set ggplot2 them
ggplot2::theme_set(theme_classic(base_size = 20))
# standard error
se <- function(x) sd(x) / sqrt(length(x))
```
```{r, echo=FALSE}
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
```
# Site characteristics
```{r}
soil <- read.csv("./data/raw/soils_effex.csv")
```
```{r, eval = TRUE}
preds <- names(soil)[7:17]
soil$pool.n <- factor(soil$pool.n)
levels(soil$pool.n) <- c("no.n", "n")
for (i in preds){
form <- as.formula(paste(i, "~ pool.n + pool.p"))
print(i)
mod <-
brm(
form,
data = soil,
chains = chains,
family = gaussian(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = c(
prior(normal(0, 10), "b"),
prior(normal(0, 50), "Intercept")),
file = paste0("./data/processed/", i, "_pooled_model"),
file_refit = "on_change"
)
}
```
```{r, results='asis', eval = TRUE}
preds <- names(soil)[7:17]
for (i in preds){
extended_summary(read.file = paste0("./data/processed/", i, "_pooled_model.rds"), gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
}
```
## Microbial CNP
```{r}
soil_master <- read.csv("./data/raw/strimastersoils.csv")
```
# Nutrient change
## Nitrogen
Content
```{r}
nutr <- read.csv("./data/raw/litter-nutrients-mas.csv")
nutr$plot.f <- as.factor(nutr$plot)
nutr$days.sc <- scale(nutr$days)
nutr$litter.n.content.prop.initial <- nutr$litter.n.content.perc.initial / 100
nutr$litter.n.content.prop.initial <- ifelse(nutr$litter.n.content.prop.initial >= 1, 0.99999, nutr$litter.n.content.prop.initial)
# remove plot 9
sub.nutr <- nutr[nutr$plot != 9, ]
agg_n <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, mean)
agg_n$sd <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, sd)$litter.n.content.perc.initial
agg_n$se <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, se)$litter.n.content.perc.initial
agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_n, aes(x = days, y = litter.n.content.perc.initial, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = litter.n.content.perc.initial + se, ymin = litter.n.content.perc.initial - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_n$days),
labels = unique(agg_n$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_N_content.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_N_content.tiff)
```{r, eval = TRUE}
fit.n <- brm(
litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f),
data = nutr,
chains = chains,
cores = chains,
family = Beta(link = "logit"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
prior = prior,
file = "./data/processed/litter.n.content.rem_model",
file_refit = "on_change"
)
```
```{r, results='asis', warning=FALSE}
extended_summary(read.file = "./data/processed/litter.n.content.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Concentration
```{r}
agg_conc_n <- aggregate(perc.n ~ colecta + treat + days, nutr, mean)
agg_conc_n$sd <- aggregate(perc.n ~ colecta + treat + days, nutr, sd)$perc.n
agg_conc_n$se <- aggregate(perc.n ~ colecta + treat + days, nutr, se)$perc.n
agg_conc_n$treat <- factor(agg_conc_n$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_conc_n, aes(x = days, y = perc.n, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.n + se, ymin = perc.n - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_n$days),
labels = unique(agg_conc_n$days)) + theme(legend.position = c(0.3, 0.8))
ggsave("./output/Litter_N_concentration.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_N_concentration.tiff)
```{r, eval = TRUE}
nutr$prop.n <- nutr$perc.n / 100
fit.perc.n <- brm(
prop.n ~ treat * days.sc + (1 | plot.f),
data = nutr,
chains = chains,
core = chains,
family = Beta(link = "logit"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
prior = prior,
file = "./data/processed/litter.n.perc_model",
file_refit = "on_change"
)
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/litter.n.perc_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Proportion N
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$prop.n <- sub.nutr$perc.n / 100
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit_np_treat_by_rem_leave <- brm(
prop.n ~ pool.n * days.sc + pool.p * days.sc + (1 | plot),
data = sub.nutr,
chains = chains,
iter = iters, backend = "cmdstanr",
family = Beta(link = "logit"),
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/n_perc_by_pooled_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/n_perc_by_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Content N
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$prop.n <- sub.nutr$perc.n / 100
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit_np_treat_by_rem_leave <- brm(
litter.n.content.prop.initial ~ pool.n * days.sc + pool.p * days.sc + (1 | plot),
data = sub.nutr,
chains = chains,
family = Beta(link = "logit"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/n_content_by_pooled_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/n_content_by_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Phosphorus
Content
```{r}
nutr$litter.p.content.prop.initial <- nutr$litter.p.content.perc.initial / 100
# excluding plot 9
agg_p <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, mean)
agg_p$sd <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, sd)$litter.p.content.perc.initial
agg_p$se <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, se)$litter.p.content.perc.initial
agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_p, aes(x = days, y = litter.p.content.perc.initial, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = litter.p.content.perc.initial + se, ymin = litter.p.content.perc.initial - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter P content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_p$days),
labels = unique(agg_p$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_P_content.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_P_content.tiff)
```{r, eval = TRUE}
fit.p <- brm(
litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f),
data = nutr[nutr$litter.p.content.prop.initial < 1, ],
chains = chains,
prior = prior,
cores = chains,
iter = iters, backend = "cmdstanr",
family = Beta(link = "logit"),
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/litter.p.content.rem_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.content.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$prop.n <- sub.nutr$perc.n / 100
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
sub.nutr$litter.p.content.prop.initial <- sub.nutr$litter.p.content.perc.initial / 100
fit.p <- brm(
litter.p.content.prop.initial ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f),
data = sub.nutr[sub.nutr$litter.p.content.prop.initial < 1, ],
family = Beta(link = "logit"),
chains = chains,
prior = prior,
cores = chains,
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/litter.p.content.perc.initial_pooled_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.content.perc.initial_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Concentration
```{r}
agg_conc_p <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, mean)
agg_conc_p$sd <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, sd)$perc.p
agg_conc_p$se <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, se)$perc.p
agg_conc_p$treat <- factor(agg_conc_p$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_conc_p, aes(x = days, y = perc.p, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.p + se, ymin = perc.p - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter P concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_p$days),
labels = unique(agg_conc_p$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_P_concentration.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_P_concentration.tiff)
```{r, eval = TRUE}
# convert to proportions to use beta distribution
sub.nutr$prop.p <- sub.nutr$perc.p / 100
fit.perc.p <- brm(
prop.p ~ treat * days.sc + (1 |
plot.f),
data = sub.nutr,
chains = chains,
cores = chains,
family = Beta(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
prior = prior,
file = "./data/processed/litter.p.perc_model",
file_refit = "on_change"
)
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.perc_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Pooled P
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit.perc.p <- brm(prop.p ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, cores = chains, family = Beta(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), prior = prior, file = "./data/processed/litter.p.perc_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.perc_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
<div class="alert alert-success">
# Takeaways {.unnumbered .unlisted}
- Litter P content is significantly lower in plus N treatment plots than in control plot, after accounting for variation explained by time
</div>
## Remaining litter
```{r, eval = TRUE, out.width = "100%", echo = FALSE, fig.align= "center"}
dat <- read_excel("./data/raw/litter_data.xlsx")
dat$plot.f <- as.factor(dat$plot)
dat$prop.litter.rem <- dat$perc.litter.rem / 100
dat$days.sc <- scale(dat$days)
agg_rem <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, mean)
agg_rem$sd <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, sd)$perc.litter.rem
agg_rem$se <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, se)$perc.litter.rem
agg_rem <- agg_rem[order(agg_rem$trat), ]
pd <- position_dodge(15)
ggplot(agg_rem, aes(x = days, y = perc.litter.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.litter.rem + se, ymin = perc.litter.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter mass remaining (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem$days),
labels = unique(agg_rem$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_mass_remaining.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_mass_remaining.tiff)
```{r, eval = TRUE, out.width = "100%", echo = FALSE, fig.align= "center"}
agg_rem_n <- aggregate(perc.litter.rem ~ colecta + pool.n + days, dat, mean)
agg_rem_n$sd <- aggregate(perc.litter.rem ~ colecta + pool.n + days, dat, sd)$perc.litter.rem
agg_rem_n$se <- aggregate(perc.litter.rem ~ colecta + pool.n + days, dat, se)$perc.litter.rem
agg_rem_p <- aggregate(perc.litter.rem ~ colecta + pool.p + days, dat, mean)
agg_rem_p$sd <- aggregate(perc.litter.rem ~ colecta + pool.p + days, dat, sd)$perc.litter.rem
agg_rem_p$se <- aggregate(perc.litter.rem ~ colecta + pool.p + days, dat, se)$perc.litter.rem
names(agg_rem_n)[2] <- names(agg_rem_p)[2] <- "pool"
agg_rem_p$treatment <- "p"
agg_rem_n$treatment <- "n"
agg_rem_pooled <- rbind(agg_rem_p, agg_rem_n)
agg_rem_pooled <- agg_rem_pooled[order(agg_rem_pooled$pool), ]
pd <- position_dodge(15)
ggplot(agg_rem_pooled, aes(x = days, y = perc.litter.rem, color = pool)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.litter.rem + se, ymin = perc.litter.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("+N", "-N", "+P", "-P")) +
labs(x= "Time (days)", y = "Remaining mass (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem_pooled$days),
labels = unique(agg_rem_pooled$days)) +
facet_wrap(~ treatment, nrow = 2)
ggsave("./output/Litter_mass_remaining_pooled_treatments.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_mass_remaining_pooled_treatments.tiff)
```{r, eval = TRUE}
fit <- brm(prop.litter.rem ~ trat * days.sc + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model", file_refit = "on_change")
# dat$days.fc <- as.numeric(as.factor(dat$days))
#
# # monotonic effect of time
# fit_mo <- brm(prop.litter.rem ~ trat * mo(days.fc) + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model_monotonic", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.litter.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r}
dat$pool.n <- factor(dat$pool.n)
levels(dat$pool.n) <- c("no.n", "n")
fit <- brm(prop.litter.rem ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_pooled_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.litter.rem_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Remaining wood
```{r, eval = TRUE, out.width = "100%", echo = FALSE, fig.align= "center"}
agg_rem_w <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, mean)
agg_rem_w$sd <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, sd)$perc.wood.rem
agg_rem_w$se <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, se)$perc.wood.rem
pd <- position_dodge(15)
ggplot(agg_rem_w, aes(x = days, y = perc.wood.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.wood.rem + se, ymin = perc.wood.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Wood mass remaining (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem_w$days),
labels = unique(agg_rem_w$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/wood_mass_remaining.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/wood_mass_remaining.tiff)
```{r, eval = TRUE}
dat$prop.wood.rem <- dat$perc.wood.rem / 100
fit2 <- brm(
prop.wood.rem ~ trat * days.sc + (1 | plot.f),
data = dat,
chains = chains,
family = Beta(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/prop.wood.rem_model",
file_refit = "on_change"
)
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.wood.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
dat$prop.wood.rem <- dat$perc.wood.rem / 100
dat$pool.n <- factor(dat$pool.n)
levels(dat$pool.n) <- c("no.n", "n")
fit2 <- brm(prop.wood.rem ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.wood.rem_pooled_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.wood.rem_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# K
## Litter
Correlation between Silvia's and Andrea's K
```{r, eval = TRUE}
k_vals <- read.csv("./data/raw/k-values-corr.csv")
names(k_vals)[names(k_vals) %in% c("k.sil.litt.corr", "k.sil.wood.corr")] <- c("k.sil.litt", "k.sil.wood")
k_vals <- k_vals[complete.cases(k_vals), ]
```
Comparing all treatments vs control
```{r}
agg_kvals <- aggregate(k.sil.litt ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.litt ~ Treatment, k_vals, se)$k.sil.litt
agg_kvals$sd <- aggregate(k.sil.litt ~ Treatment, k_vals, sd)$k.sil.litt
agg_kvals$n <- aggregate(k.sil.litt ~ Treatment, k_vals, length)$k.sil.litt
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
# composed box plot
ggplot(k_vals, aes(x = Treatment, y = k.sil.litt)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
## adjust height
width = .6,
.width = 0,
## move geom to the right
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)") +
# ylim(c(-0.39, 0.145)) +
geom_text(data = agg_kvals, aes(y = rep(0.4, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
"NP" = "+NP", "P" = "+P"))
ggsave("./output/litter_decomposition_k_by_treatment.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/litter_decomposition_k_by_treatment.tiff)
```{r, eval = TRUE}
fit_k.litter <- brm(k.sil.litt ~ Treatment + (1 | idplot), data = k_vals, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_model.rds")
```
pooled N and P
```{r, eval = TRUE}
k_vals$pool.p <- ifelse(k_vals$Treatment %in% c("P", "NP"), "p", "no.p")
k_vals$pool.n <- ifelse(k_vals$Treatment %in% c("N", "NP"), "n", "no.n")
k_vals$pool.n <- factor(k_vals$pool.n)
levels(k_vals$pool.n) <- c("no.n", "n")
fit_k.litter <- brm(k.sil.litt ~ pool.p + pool.n + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_pooled_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_pooled_model.rds")
```
Phosphorus vs no-phosphorus
```{r, eval = TRUE}
fit_k.litter.p.np <- brm(k.sil.litt ~ p.treat + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", , control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_p_nop_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_litter_p_nop_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r}
agg_kvals <- aggregate(k.sil.litt ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.litt ~ p.treat, k_vals, se)$k.sil.litt
agg_kvals$sd <- aggregate(k.sil.litt ~ p.treat, k_vals, sd)$k.sil.litt
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(k.sil.litt ~ p.treat, k_vals, length)$k.sil.litt
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$p.treat <- ifelse(k_vals$p.treat == "no.p", "No P", "P")
k_vals$p.treat <- factor(k_vals$p.treat)
# composed box plot
rn_p_litter <- ggplot(k_vals, aes(x = p.treat, y = k.sil.litt)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_p_litter
ggsave("./output/litter_phosphorus_decomposition_k.tiff", dpi = 300)
```
Nitrogen vs no-nitrogen
```{r, eval = TRUE}
fit_k.litter.n.nn <- brm(k.sil.litt ~ n.treat + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_n_no_n_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_litter_n_no_n_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r}
agg_kvals <- aggregate(k.sil.litt ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.litt ~ n.treat, k_vals, se)$k.sil.litt
agg_kvals$sd <- aggregate(k.sil.litt ~ n.treat, k_vals, sd)$k.sil.litt
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(k.sil.litt ~ n.treat, k_vals, length)$k.sil.litt
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))
# composed box plot
rn_n_litter <- ggplot(k_vals, aes(x = n.treat, y = k.sil.litt)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_n_litter
ggsave("./output/litter_nitrogen_decomposition_k.tiff", dpi = 300)
```
## Wood
Comparing all treatments vs control
```{r}
agg_kvals <- aggregate(k.sil.wood ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.wood ~ Treatment, k_vals, se)$k.sil.wood
agg_kvals$sd <- aggregate(k.sil.wood ~ Treatment, k_vals, sd)$k.sil.wood
agg_kvals$n <- aggregate(k.sil.wood ~ Treatment, k_vals, length)$k.sil.wood
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
# composed box plot
ggplot(k_vals, aes(x = Treatment, y = k.sil.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)") +
# ylim(c(-0.39, 0.145)) +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
"NP" = "+NP", "P" = "+P"))
ggsave("./output/wood_decomposition_k_by_treatment.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/wood_decomposition_k_by_treatment.tiff)
```{r, eval = TRUE}
fit_k.wood.p.treat <- brm(k.sil.wood ~ Treatment + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters * 2, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_treatment_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_treatment_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
pooled N and P
```{r, eval = TRUE}
k_vals$pool.p <- ifelse(k_vals$Treatment %in% c("P", "NP"), "p", "no.p")
k_vals$pool.n <- ifelse(k_vals$Treatment %in% c("N", "NP"), "n", "no.n")
k_vals$pool.n <- factor(k_vals$pool.n)
levels(k_vals$pool.n) <- c("no.n", "n")
fit_k.wood.p.pooled <- brm(k.sil.wood ~ pool.p + pool.n + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_pooled_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_wood_pooled_model.rds")
```
## Phosphorus vs no-phosphorus
```{r}
agg_kvals <- aggregate(k.sil.wood ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.wood ~ p.treat, k_vals, se)$k.sil.wood
agg_kvals$sd <- aggregate(k.sil.wood ~ p.treat, k_vals, sd)$k.sil.wood
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(k.sil.wood ~ p.treat, k_vals, length)$k.sil.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
# k_vals$p.treat <- ifelse(k_vals$p.treat == "no.p", "No P", "P")
# k_vals$p.treat <- factor(k_vals$p.treat)
# composed box plot
rn_p_wood <- ggplot(k_vals, aes(x = p.treat, y = k.sil.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_p_wood
ggsave("./output/phosphorus_decomposition_k.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/phosphorus_decomposition_k.tiff)
```{r, eval = TRUE}
fit_k.wood.p.no.p <- brm(k.sil.wood ~ p.treat + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_p_no_p_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_p_no_p_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Nitrogen vs no-nitrogen
```{r}
agg_kvals <- aggregate(k.sil.wood ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(k.sil.wood ~ n.treat, k_vals, se)$k.sil.wood
agg_kvals$sd <- aggregate(k.sil.wood ~ n.treat, k_vals, sd)$k.sil.wood
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(k.sil.wood ~ n.treat, k_vals, length)$k.sil.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))
# composed box plot
rn_n_wood <- ggplot(k_vals, aes(x = n.treat, y = k.sil.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_n_wood
ggsave("./output/nitrogen_decomposition_k.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nitrogen_decomposition_k.tiff)
```{r, eval = TRUE}
fit_k.wood.n.no.n <- brm(k.sil.wood ~ n.treat + (1 | idplot), data = k_vals, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_n_no_n_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_n_no_n_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# Nutrient content by remaining litter mass
## Nitrogen
```{r}
nutr$prom.hoja.reman.sc <- scale(nutr$prom.hoja.reman)
ggplot(data = nutr, aes(x = prom.hoja.reman, y = perc.n, color = treat)) +
geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter N concentration (%)", color = "Treatment") +
scale_color_viridis_d(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
scale_x_reverse()
ggsave("./output/nitrogen_by_litter_mass.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nitrogen_by_litter_mass.tiff)
```{r, eval = TRUE}
nutr$prop.n <- nutr$perc.n / 100
fit_pern_by_rem_leave <- brm(prop.n ~ prom.hoja.reman.sc * treat + (1 | plot), data = nutr, chains = chains, family = Beta(link = "logit"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
fit_pern_by_rem_leave <- brm(prop.n ~ prom.hoja.reman.sc * treat + (1 | plot), data = nutr, chains = chains, family = Beta(link = "logit"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Phosphorus
```{r}
ggplot(data = sub.nutr, aes(x = prom.hoja.reman, y = perc.p, color = treat)) +
geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter P concentration (%)", color = "Treatment") +
scale_color_viridis_d(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
scale_x_reverse()
ggsave("./output/phosphorus_by_litter_mass.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/phosphorus_by_litter_mass.tiff)
```{r, eval = TRUE}
sub.nutr$prom.hoja.reman.sc <- scale(sub.nutr$prom.hoja.reman)
sub.nutr$prop.p <- sub.nutr$perc.p / 100
fit_perp_by_rem_leave <- brm(prop.p ~ prom.hoja.reman.sc * treat + (1 | plot), data = sub.nutr, chains = chains, family = Beta(link = "logit"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/p_per_by_leave_rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/p_per_by_leave_rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
<!-- light brown box -->
<div class="alert alert-success">
# Takeaways {.unnumbered .unlisted}
- Nitrogen percentage, but no Phosphorus percentage, increases along with remaining litter mass
</div>
# Nutrient ratios
## C:N
```{r}
# excluding plot 9
agg_cn <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, mean)
agg_cn$sd <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, sd)$cn.mol.kg
agg_cn$se <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, se)$cn.mol.kg
agg_cn$treat <- factor(agg_cn$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_cn, aes(x = days, y = cn.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = cn.mol.kg + se, ymin = cn.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter C:N ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_cn$days),
labels = unique(agg_cn$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/cn_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/cn_ratio_through_time.tiff)
```{r, eval = TRUE}
fit.cn <- brm(cn.mol.kg ~ treat * days.sc + (1 | plot.f), data = nutr, chains = chains, family = Gamma(link = "log"), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/cn_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cn_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit.cn <- brm(cn.mol.kg ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/cn_pooled_model2", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cn_pooled_model2.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## N:P
```{r}
# excluding plot 9
agg_np <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, mean)
agg_np$sd <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, sd)$np.mol.kg
agg_np$se <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, se)$np.mol.kg
agg_np$treat <- factor(agg_np$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_np, aes(x = days, y = np.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = np.mol.kg + se, ymin = np.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N:P ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_np$days),
labels = unique(agg_np$days)) + theme(legend.position = c(0.2, 0.8))
ggsave("./output/np_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/np_ratio_through_time.tiff)
```{r, eval = TRUE}
fit.np <- brm(np.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/np_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/np_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit.cn <- brm(np.mol.kg ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/np_pooled_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/np_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## C:P
```{r}
# excluding plot 9
agg_cp <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, mean)
agg_cp$sd <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, sd)$cp.mol.kg
agg_cp$se <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, se)$cp.mol.kg
agg_cp$treat <- factor(agg_cp$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_cp, aes(x = days, y = cp.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = cp.mol.kg + se, ymin = cp.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter C:P ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_cp$days),
labels = unique(agg_cp$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/cp_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/cp_ratio_through_time.tiff)
```{r, eval = TRUE}
fit.cp <- brm(cp.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", cores = chains, prior = prior, file = "./data/processed/cp_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cp_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r, eval = TRUE}
sub.nutr$pool.p <- ifelse(sub.nutr$treat %in% c("P", "NP"), "p", "no.p")
sub.nutr$pool.n <- ifelse(sub.nutr$treat %in% c("N", "NP"), "n", "no.n")
sub.nutr$pool.n <- factor(sub.nutr$pool.n)
levels(sub.nutr$pool.n) <- c("no.n", "n")
fit.cn <- brm(cp.mol.kg ~ pool.n * days.sc + pool.p * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, backend = "cmdstanr", control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/cp_pooled_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cp_pooled_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
---
# Combined plots
## Remaining mass
```{r, out.width="100%", out.height="140%"}
agg_rem_w$substrate <- "Wood"
agg_rem$substrate <- "Litter"
names(agg_rem) <- names(agg_rem_w) <- c("colecta", "trat", "days", "perc.rem", "sd", "se", "substrate")
agg_rem_pooled <- rbind(agg_rem, agg_rem_w)
ggplot(agg_rem_pooled, aes(x = days, y = perc.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.rem + se, ymin = perc.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Remaining mass (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem_pooled$days),
labels = unique(agg_rem_pooled$days)) +
facet_wrap(~ substrate, nrow = 2)
ggsave("./output/remaining_mass_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/remaining_mass_combined.tiff)
## Litter content
```{r, out.width="100%", out.height="140%"}
agg_p$nutrient <- "P"
agg_n$nutrient <- "N"
names(agg_n) <- names(agg_p) <- c("colecta", "treat", "days", "perc", "sd", "se", "nutrient")
agg_np2 <- rbind(agg_n, agg_p)
ggplot(agg_np2, aes(x = days, y = perc, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc + se, ymin = perc - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_p$days),
labels = unique(agg_p$days)) +
facet_wrap(~ nutrient, nrow = 2)
ggsave("./output/litter_content_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/litter_content_combined.tiff)
## Concentration
```{r, out.width="100%", out.height="140%"}
agg_conc_p$nutrient <- "P"
agg_conc_n$nutrient <- "N"
names(agg_conc_n) <- names(agg_conc_p) <- c("colecta", "treat", "days", "perc", "sd", "se", "nutrient")
agg_conc_np <- rbind(agg_conc_n, agg_conc_p)
ggplot(agg_conc_np, aes(x = days, y = perc, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc + se, ymin = perc - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_p$days),
labels = unique(agg_conc_p$days)) +
facet_wrap(~ nutrient, nrow = 2, scales = "free_y")
ggsave("./output/concentration_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/concentration_combined.tiff)
## Nutrient ratios
```{r, out.width="100%", fig.height=12}
agg_cn$nutrients <- "C:N"
agg_cp$nutrients <- "C:P"
agg_np$nutrients <- "N:P"
names(agg_cn) <- names(agg_cp) <- names(agg_np) <- c("colecta", "treat", "days", "mol.kg", "sd", "se", "nutrients")
agg_ratios <- rbind(agg_cn[, names(agg_cn)], agg_cp[, names(agg_cn)], agg_np[, names(agg_cn)])
ggplot(agg_ratios, aes(x = days, y = mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = mol.kg + se, ymin = mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient ratios", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_ratios$days),
labels = unique(agg_ratios$days)) +
facet_wrap(~ nutrients, nrow = 3, scales = "free_y")
ggsave("./output/nutrient_ratios_combined.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nutrient_ratios_combined.tiff)
## Decomposition constant
```{r, out.width="100%", fig.height=12}
rn_n_wood_dat <- rn_n_wood$data[, c("n.treat", "k.sil.wood")]
rn_n_litter_dat <- rn_n_litter$data[, c("n.treat", "k.sil.litt")]
rn_p_wood_dat <- rn_p_wood$data[, c("p.treat", "k.sil.wood")]
rn_p_litter_dat <- rn_p_litter$data[, c("p.treat", "k.sil.litt")]
names(rn_p_wood_dat) <- names(rn_n_wood_dat) <- names(rn_p_litter_dat) <- names(rn_n_litter_dat) <- c("treatment", "k")
rn_p_wood_dat$substrate <- rn_n_wood_dat$substrate <- "wood"
rn_p_litter_dat$substrate <- rn_n_litter_dat$substrate <- "litter"
rn_n_wood_dat$nutrient <- rn_n_litter_dat$nutrient <- "N"
rn_p_litter_dat$nutrient <- rn_p_wood_dat$nutrient <- "P"
rn_k_dat <- rbind(rn_n_wood_dat, rn_p_wood_dat, rn_n_litter_dat, rn_p_litter_dat)
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
# geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
facet_grid(~ substrate)
ggsave("./output/decomposition_k_combined_v1.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/decomposition_k_combined_v1.tiff)
```{r}
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
# geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
facet_grid(substrate ~ nutrient, scales = "free_x")
ggsave("./output/decomposition_k_combined_v12.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/decomposition_k_combined_v2.tiff)
# Carbon inputs/outputs
## Litter
### N
```{r}
# read data
c_input_litter <- read.csv("./data/raw/annual_litter_C_input.csv")
# same with ggplot2
ggplot(c_input_litter, aes(x = n.pooled, y = cr.total)) +
geom_boxplot() +
labs(y = "Carbon input (Mg C ha-1 yr-1)", x = "N treatment") +
theme_classic(base_size = 18) +
scale_x_discrete(labels=c("no.n" = "-N", "with.n" = "+N"))
# mod1 <- lm(cr.total ~ n.pooled, data = c_input_litter[complete.cases(c_input_litter), ])
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"))
litter_input_n <- brm(
cr.total ~ n.pooled,
data = c_input_litter[complete.cases(c_input_litter), c("cr.total", "n.pooled")],
chains = chains,
family = gaussian(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/litter_input_n_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(fit = litter_input_n, highlight = TRUE, remove.intercepts = TRUE)
```
### P
```{r}
ggplot(c_input_litter, aes(x = p.pooled, y = cr.total)) +
geom_boxplot() +
labs(y = "Carbon input (Mg C ha-1 yr-1)", x = "P treatment") +
theme_classic(base_size = 18) +
scale_x_discrete(labels=c("no.p" = "-P", "with.p" = "+P"))
litter_input_p <- brm(
cr.total ~ p.pooled,
data = c_input_litter[complete.cases(c_input_litter), c("cr.total", "p.pooled")],
chains = chains,
family = gaussian(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/litter_input_p_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(fit = litter_input_p, highlight = TRUE, remove.intercepts = TRUE)
```
## CO2 output
### N
```{r}
c_avg_co2 <- read.csv("./data/raw/co2_plot_averages.csv")
ggplot(c_avg_co2, aes(x = pool.n, y = annual.flux)) +
geom_boxplot() +
labs(y = "Carbon output (Mg C ha-1 yr-1)", x = "N treatment") +
theme_classic(base_size = 18)
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"))
co2_output_n <- brm(
annual.flux ~ pool.n,
data = c_avg_co2[complete.cases(c_avg_co2), ],
chains = chains,
family = gaussian(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/co2_output_n_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(fit = co2_output_n, highlight = TRUE, remove.intercepts = TRUE)
```
### P
```{r}
ggplot(c_avg_co2, aes(x = pool.p, y = annual.flux)) +
geom_boxplot() +
labs(y = "Carbon output (Mg C ha-1 yr-1)", x = "P treatment") +
theme_classic(base_size = 18)
co2_output_p <- brm(
annual.flux ~ pool.p,
data = c_avg_co2[complete.cases(c_avg_co2), ],
chains = chains,
family = gaussian(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/co2_output_p_model",
file_refit = "on_change"
)
```
```{r, eval = TRUE, results='asis'}
extended_summary(fit = co2_output_p, highlight = TRUE, remove.intercepts = TRUE)
```
# Microbial biomass
## N
```{r}
soil_mic <- read.csv("./data/raw/strimastersoils.csv")
# names(soil_mic)
soil_mic$BS <- soil_mic$BS / 100
soil_mic$Al.sat <- soil_mic$Al.sat / 100
variables <- c("mic.c", "mic.n", "mic.p", "mic.cn", "mic.cp", "mic.np",
"AG", "BG", "XYL", "CEL", "NAG", "LAP", "MUP", "BIS", "S",
"BG.NAG", "BG.MUP", "NAG.MUP", "BG.S", "NAG.S", "MUP.S",
"ph.h2o", "ph.cacl2", "tc.g.kg.18", "tn.g.kg.18", "tp.g.kg.18",
"CN", "CP", "NP", "k2so4.C", "nitrate.2017", "resin.p",
"ox.Al", "ox.Fe", "Al", "Ca", "Fe", "K", "Mg", "Mn", "Na",
"TEB", "ECEC", "BS", "Al.sat")
variables[!variables %in% names(soil_mic)]
for (i in variables){
soil_mic$var <- soil_mic[,i]
gg <- ggplot(soil_mic, aes(x = treat.pool.n, y = var)) +
geom_boxplot() +
geom_jitter() +
labs(y = i, x = "N treatment") +
theme_classic(base_size = 18)
print(gg)
}
soil_mic$var <- NULL
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"))
for (i in variables){
soil_mic$var <- soil_mic[,i]
if(i %in% c("ph.cacl2", "ph.h2o"))
fam <- gaussian() else
fam <- Gamma(link = "log")
if(i %in% c("BS", "Al.sat"))
fam <- Beta(link = "logit") else
if (any(soil_mic$var < 0))
soil_mic$var <- soil_mic$var + min(abs(soil_mic$var)) + 0.0001
co2_output_n <- brm(
var ~ treat.pool.n,
data = soil_mic,
chains = chains,
family = fam,
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = paste("./data/processed/soil_", i, "_model", sep = ""),
file_refit = "on_change"
)
}
```
```{r, eval = TRUE, results='asis'}
for (i in variables) {
mod <- paste("./data/processed/soil_", i, "_model.rds", sep = "")
# print(paste("###", i))
# print("<br>")
extended_summary(read.file = mod,
highlight = TRUE,
remove.intercepts = TRUE)
}
```
## P
```{r}
for (i in variables){
soil_mic$var <- soil_mic[,i]
gg <- ggplot(soil_mic, aes(x = treat.pool.p, y = var)) +
geom_boxplot() +
geom_jitter() +
labs(y = i, x = "N treatment") +
theme_classic(base_size = 18)
print(gg)
}
soil_mic$var <- NULL
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"))
for (i in variables){
soil_mic$var <- soil_mic[,i]
fam <- Gamma(link = "log")
if(i %in% c("BS", "Al.sat"))
fam <- Beta(link = "logit") else
soil_mic$var <- soil_mic$var + min(abs(soil_mic$var)) + 0.0001
co2_output_n <- brm(
var ~ treat.pool.p,
data = soil_mic,
chains = chains,
family = fam,
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = paste("./data/processed/soil_", i, "_p_model", sep = ""),
file_refit = "on_change"
)
}
```
```{r, eval = TRUE, results='asis'}
for (i in variables) {
# print(paste("###", i))
#
# print("<br>")
mod <- paste("./data/processed/soil_", i, "_p_model.rds", sep = "")
extended_summary(read.file = mod,
highlight = TRUE,
remove.intercepts = TRUE)
}
```
## BIS by Nitrate rate
```{r, results='asis'}
soil_mic$Nrat <- soil_mic$nitrate.2017 / soil_mic$resin.p
ggplot(soil_mic, aes(x = Nrat, y = BIS)) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic(base_size = 18)
BIS_Nitrate_mod <- brm(
BIS ~ Nrat,
data = soil_mic,
chains = chains,
family = Gamma(link = "log"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/soil_BIS_by_Nitrate_rate_model",
file_refit = "on_change"
)
extended_summary(fit = BIS_Nitrate_mod,
highlight = TRUE,
remove.intercepts = TRUE)
```
# CO2 data
## All sampling events
```{r, results='asis'}
co2mean <- read.csv("./data/raw/CO2_dat_mean.csv")
co2mean$date4 <- as.Date(co2mean$date4)
co2mean$day <- co2mean$date4 - min(co2mean$date4)
exp_flux_pool_mod <- brm(
exp.flux ~ pool.n * day + (1|id.plot/collar.generic),
data = co2mean,
chains = chains,
family = student(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/exp_flux_pool_n_model",
file_refit = "on_change"
)
extended_summary(fit = exp_flux_pool_mod,
highlight = TRUE,
remove.intercepts = TRUE)
exp_flux_pool_p_mod <- brm(
exp.flux ~ pool.p * day + (1|id.plot/collar.generic),
data = co2mean,
chains = chains,
family = student(),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/exp_flux_pool_p_model",
file_refit = "on_change"
)
extended_summary(fit = exp_flux_pool_p_mod,
highlight = TRUE,
remove.intercepts = TRUE)
```
## Only first and last sampling events
```{r, results='asis'}
co2mean <- read.csv("./data/raw/CO2_dat_mean.csv")
co2mean$date4 <- as.Date(co2mean$date4)
co2mean$day <- co2mean$date4 - min(co2mean$date4)
sub_co2mean <- co2mean[co2mean$day %in% c(0, max(co2mean$day)), ]
sub_co2mean$day_f <- ifelse(sub_co2mean$day == 0, "first", "last")
exp_flux_pool_fl_mod <- brm(
exp.flux ~ pool.n * day_f + (1|id.plot/collar.generic),
data = sub_co2mean,
chains = chains,
family = Gamma(link = "log"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/exp_flux_pool_n_first_last_model",
file_refit = "on_change"
)
extended_summary(fit = exp_flux_pool_fl_mod,
highlight = TRUE,
remove.intercepts = TRUE)
exp_flux_pool_p_fl_mod <- brm(
exp.flux ~ pool.p * day_f + (1|id.plot/collar.generic),
data = sub_co2mean,
chains = chains,
family = Gamma(link = "log"),
iter = iters, backend = "cmdstanr",
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = chains,
prior = prior,
file = "./data/processed/exp_flux_pool_p_first_last_model",
file_refit = "on_change"
)
extended_summary(fit = exp_flux_pool_p_fl_mod,
highlight = TRUE,
remove.intercepts = TRUE)
```
# Combined model diagnostics
```{r, eval = FALSE}
check_rds_fits(path = "./data/processed", html = TRUE)
```
<!-- add packages used, system details and versions -->
# Session information {.unnumbered .unlisted}
```{r session info, echo=F}
sessionInfo()
```
## Phosphorus vs no-phosphorus
