#install.packages(c(readxl, brms, loo, forecast, ggplot2, ggpubr, reshape2, rstan))
rm(list=ls(all=T)) #Clean work environment
library(readxl)
library(rstan)
library(brms)
library(loo)
library(forecast)
library(ggplot2)
library(ggpubr)
library(reshape2)
library(dplyr)
#READ TRAINNING DATA: LABILITY EXPERIMENTS WITH NATURAL CONDITIONS
setwd("D:/ebals/Documents/Alumnos/Alumnos/Marcela/paper_DOM/Modelos/Modelo_def") #change to your file location
data <- read_xlsx("Data_base_Bastidas_Schenone.xlsx",
sheet = 'Experiments') #change to your file name
#Rename columns for analysis
data <- data %>%
rename(time = "exp. time", INI = "DOC0 (mg L-1)", TDP = "TDP (µmol L-1)", TDN = "TDN (µmol L-1)", DOC = "DOC (mg L-1)")
#BAYESIAN REGRESSION MODELS USING STAN
chains = 3 #number of Markovian chains to run
iter = 3000 #iterations per Markovian chain
warmup = 1500 #burn-in iterations before parameter sampling
cores = 6 #set a computer core for each chain if possible
seed = 7777 #make model run reproducible
control = list(adapt_delta = 0.99, max_treedepth = 15) # Control parameter for chain iterations
#Null model
m_null <- brm(formula = bf(DOC ~ lab * exp(-k * time) + res, #model equation
lab ~ 1, k ~ 1, res ~ 1, nl = TRUE), #parameter dependencies
data = data, family = gaussian(),
prior = c(set_prior("normal(0, 1)", nlpar = "lab", lb = 0), #prior constraints
set_prior("normal(0, 1)", nlpar = "k", lb = 0),
set_prior("normal(0, 1)", nlpar = "res", lb = 0)),
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, #save parameter sampling information
file = "m_null", #save model output
control = control) #control chain steps
print(m_null, digits = 5) #to see posterior summary results
#Priors for the alternative models
priors = c(set_prior("normal(0.5, .25)", nlpar = "plab", lb = 0),
set_prior("normal(0, 1)", nlpar = "k", lb = 0),
set_prior("normal(0.5, .25)", nlpar = "pres", lb = 0))
#Control model
m_cntr <- brm(formula = bf(DOC ~ (plab * exp(-k * time) + pres) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_cntr", control = control)
print(m_cntr, digits = 5)
#Model with TDP in k and TDP in pres
m_P_P <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_P", control = control)
print(m_P_P, digits = 5)
#Model with TDP in k and TDN in pres
m_P_N <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_N", control = control)
print(m_P_N, digits = 5)
#Model with TDP in k and Cprot in pres
m_P_C <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_C", control = control)
print(m_P_C, digits = 5)
#Model with TDP in k and Cprot * TDP in pres
m_P_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_PC", control = control)
print(m_P_PC, digits = 5)
#Model with TDP in k and Cprot * sqrt(TDP) in pres
m_P_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_sPC", control = control)
print(m_P_sPC, digits = 5)
#Model with TDP in k and Cprot * TDN in pres
m_P_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_NC", control = control)
print(m_P_NC, digits = 5)
#Model with TDP in k and Cprot * sqrt(TDN) in pres
m_P_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDP * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_P_sNC", control = control)
print(m_P_sNC, digits = 5)
#Model with TDN in k and TDP in pres
m_N_P <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_P", control = control)
print(m_N_P, digits = 5)
#Model with TDN in k and TDN in pres
m_N_N <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_N", control = control)
print(m_N_N, digits = 5)
#Model with TDN in k and Cprot in pres
m_N_C <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_C", control = control)
print(m_N_C, digits = 5)
#Model with TDN in k and Cprot * TDP in pres
m_N_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_PC", control = control)
print(m_N_PC, digits = 5)
#Model with TDN in k and Cprot * sqrt(TDP) in pres
m_N_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_sPC", control = control)
print(m_N_sPC, digits = 5)
#Model with TDN in k and Cprot * TDN in pres
m_N_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_NC", control = control)
print(m_N_NC, digits = 5)
#Model with TDN in k and Cprot * sqrt(TDN) in pres
m_N_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * TDN * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_N_sNC", control = control)
print(m_N_sNC, digits = 5)
#Model with Cprot in k and TDP in pres
m_C_P <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_P", control = control)
print(m_C_P, digits = 5)
#Model with Cprot in k and TDN in pres
m_C_N <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_N", control = control)
print(m_C_N, digits = 5)
#Model with Cprot in k and Cprot in pres
m_C_C <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_C", control = control)
print(m_C_C, digits = 5)
#Model with Cprot in k and Cprot * TDP in pres
m_C_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_PC", control = control)
print(m_C_PC, digits = 5)
#Model with Cprot in k and Cprot * sqrt(TDP) in pres
m_C_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_sPC", control = control)
print(m_C_sPC, digits = 5)
#Model with Cprot in k and Cprot * TDN in pres
m_C_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_NC", control = control)
print(m_C_NC, digits = 5)
#Model with Cprot in k and Cprot * sqrt(TDN) in pres
m_C_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_C_sNC", control = control)
print(m_C_sNC, digits = 5)
#Model with Cprot * TDP in k and TDP in pres
m_PC_P <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_P", control = control)
print(m_PC_P, digits = 5)
#Model with Cprot * TDP in k and TDN in pres
m_PC_N <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_N", control = control)
print(m_PC_N, digits = 5)
#Model with Cprot * TDP in k and Cprot in pres
m_PC_C <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_C", control = control)
print(m_PC_C, digits = 5)
#Model with Cprot * TDP in k and Cprot * TDP in pres
m_PC_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_PC", control = control)
print(m_PC_PC, digits = 5)
#Model with Cprot * TDP in k and Cprot * sqrt(TDP) in pres
m_PC_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_sPC", control = control)
print(m_PC_sPC, digits = 5)
#Model with Cprot * TDP in k and Cprot * TDN in pres
m_PC_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_NC", control = control)
print(m_PC_NC, digits = 5)
#Model with Cprot * TDP in k and Cprot * sqrt(TDN) in pres
m_PC_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDP * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_PC_sNC", control = control)
print(m_PC_sNC, digits = 5)
#Model with Cprot * sqrt(TDP) in k and TDP in pres
m_sPC_P <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_P", control = control)
print(m_sPC_P, digits = 5)
#Model with Cprot * sqrt(TDP) in k and TDN in pres
m_sPC_N <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_N", control = control)
print(m_sPC_N, digits = 5)
#Model with Cprot * sqrt(TDP) in k and Cprot in pres
m_sPC_C <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_C", control = control)
print(m_sPC_C, digits = 5)
#Model with Cprot * sqrt(TDP) in k and Cprot * TDP in pres
m_sPC_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_PC", control = control)
print(m_sPC_PC, digits = 5)
#Model with Cprot * sqrt(TDP) in k and Cprot * sqrt(TDP) in pres
m_sPC_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_sPC", control = control)
print(m_sPC_sPC, digits = 5)
#Model with Cprot * sqrt(TDP) in k and Cprot * TDN in pres
m_sPC_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_NC", control = control)
print(m_sPC_NC, digits = 5)
#Model with Cprot * sqrt(TDP) in k and Cprot * sqrt(TDN) in pres
m_sPC_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDP) * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sPC_sNC", control = control)
print(m_sPC_sNC, digits = 5)
#Model with Cprot * TDN in k and TDP in pres
m_NC_P <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_P", control = control)
print(m_NC_P, digits = 5)
#Model with Cprot * TDN in k and TDN in pres
m_NC_N <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_N", control = control)
print(m_NC_N, digits = 5)
#Model with Cprot * TDN in k and Cprot in pres
m_NC_C <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_C", control = control)
print(m_NC_C, digits = 5)
#Model with Cprot * TDN in k and Cprot * TDP in pres
m_NC_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_PC", control = control)
print(m_NC_PC, digits = 5)
#Model with Cprot * TDN in k and Cprot * sqrt(TDP) in pres
m_NC_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_sPC", control = control)
print(m_NC_sPC, digits = 5)
#Model with Cprot * TDN in k and Cprot * TDN in pres
m_NC_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_NC", control = control)
print(m_NC_NC, digits = 5)
#Model with Cprot * TDN in k and Cprot * sqrt(TDN) in pres
m_NC_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * TDN * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_NC_sNC", control = control)
print(m_NC_sNC, digits = 5)
#Model with Cprot * sqrt(TDN) in k and TDP in pres
m_sNC_P <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + TDP)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_P", control = control)
print(m_sNC_P, digits = 5)
#Model with Cprot * sqrt(TDN) in k and TDN in pres
m_sNC_N <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + TDN)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_N", control = control)
print(m_sNC_N, digits = 5)
#Model with Cprot * sqrt(TDN) in k and Cprot in pres
m_sNC_C <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_C", control = control)
print(m_sNC_C, digits = 5)
#Model with Cprot * sqrt(TDN) in k and Cprot * TDP in pres
m_sNC_PC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + TDP * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_PC", control = control)
print(m_sNC_PC, digits = 5)
#Model with Cprot * sqrt(TDN) in k and Cprot * sqrt(TDP) in pres
m_sNC_sPC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + sqrt(TDP) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_sPC", control = control)
print(m_sNC_sPC, digits = 5)
#Model with Cprot * sqrt(TDN) in k and Cprot * TDN in pres
m_sNC_NC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + TDN * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_NC", control = control)
print(m_sNC_NC, digits = 5)
#Model with Cprot * sqrt(TDN) in k and Cprot * sqrt(TDN) in pres
m_sNC_sNC <- brm(formula = bf(DOC ~ (plab * exp(-k * Cprot * sqrt(TDN) * time) + pres / (1 + sqrt(TDN) * Cprot)) * INI,
plab ~ 1, k ~ 1, pres ~ 1, nl = TRUE),
data = data, family = gaussian(), prior = priors,
chains = chains, iter = iter, warmup = warmup, cores = cores, seed = seed,
sample_prior = TRUE, file = "m_sNC_sNC", control = control)
print(m_sNC_sNC, digits = 5)
#FIT COMPARISON of all 51 models using Leave One Out Cross Validation (LOOCV)
loo <- loo(m_null, m_cntr,
m_P_P, m_P_N, m_P_C, m_P_PC, m_P_sPC, m_P_NC, m_P_sNC,
m_N_P, m_N_N, m_N_C, m_N_PC, m_N_sPC, m_N_NC, m_N_sNC,
m_C_P, m_C_N, m_C_C, m_C_PC, m_C_sPC, m_C_NC, m_C_sNC,
m_PC_P, m_PC_N, m_PC_C, m_PC_PC, m_PC_sPC, m_PC_NC, m_PC_sNC,
m_sPC_P, m_sPC_N, m_sPC_C, m_sPC_PC, m_sPC_sPC, m_sPC_NC, m_sPC_sNC,
m_NC_P, m_NC_N, m_NC_C, m_NC_PC, m_NC_sPC, m_NC_NC, m_NC_sNC,
m_sNC_P, m_sNC_N, m_sNC_C, m_sNC_PC, m_sNC_sPC, m_sNC_NC, m_sNC_sNC)
# reloo=TRUE)
#save the LOOIC (effective log posterior density (ELPD) * -2) and its error
LOO_Diffs <- data.frame(m_names=row.names(loo$diffs), loo$diffs[,1:2] * -2)
#READ TEST DATASET TO MAKE PREDICTIONS
indep <- read_xlsx("Data_base_Bastidas_Schenone.xlsx", sheet = 'Enrichments') #change to your file name
#Rename columns for analysis
indep <- indep %>%
rename(time = "exp. time", INI = "DOC0 (mg L-1)", TDP = "TDP (µmol L-1)", TDN = "TDN (µmol L-1)", DOC = "DOC (mg L-1)")
indep <- indep[order(indep$Experiment),]
#PREDICTIVE skills of each model using indep test dataset
p_null <- predict(m_null, newdata = indep, type = 'response')
p_cntr <- predict(m_cntr, newdata = indep, type = 'response')
p_P_P <- predict(m_P_P, newdata = indep, type = 'response')
p_P_N <- predict(m_P_N, newdata = indep, type = 'response')
p_P_C <- predict(m_P_C, newdata = indep, type = 'response')
p_P_PC <- predict(m_P_PC, newdata = indep, type = 'response')
p_P_sPC <- predict(m_P_sPC, newdata = indep, type = 'response')
p_P_NC <- predict(m_P_NC, newdata = indep, type = 'response')
p_P_sNC <- predict(m_P_sNC, newdata = indep, type = 'response')
p_N_P <- predict(m_N_P, newdata = indep, type = 'response')
p_N_N <- predict(m_N_N, newdata = indep, type = 'response')
p_N_C <- predict(m_N_C, newdata = indep, type = 'response')
p_N_PC <- predict(m_N_PC, newdata = indep, type = 'response')
p_N_sPC <- predict(m_N_sPC, newdata = indep, type = 'response')
p_N_NC <- predict(m_N_NC, newdata = indep, type = 'response')
p_N_sNC <- predict(m_N_sNC, newdata = indep, type = 'response')
p_C_P <- predict(m_C_P, newdata = indep, type = 'response')
p_C_N <- predict(m_C_N, newdata = indep, type = 'response')
p_C_C <- predict(m_C_C, newdata = indep, type = 'response')
p_C_PC <- predict(m_C_PC, newdata = indep, type = 'response')
p_C_sPC <- predict(m_C_sPC, newdata = indep, type = 'response')
p_C_NC <- predict(m_C_NC, newdata = indep, type = 'response')
p_C_sNC <- predict(m_C_sNC, newdata = indep, type = 'response')
p_PC_P <- predict(m_PC_P, newdata = indep, type = 'response')
p_PC_N <- predict(m_PC_N, newdata = indep, type = 'response')
p_PC_C <- predict(m_PC_C, newdata = indep, type = 'response')
p_PC_PC <- predict(m_PC_PC, newdata = indep, type = 'response')
p_PC_sPC <- predict(m_PC_sPC, newdata = indep, type = 'response')
p_PC_NC <- predict(m_PC_NC, newdata = indep, type = 'response')
p_PC_sNC <- predict(m_PC_sNC, newdata = indep, type = 'response')
p_sPC_P <- predict(m_sPC_P, newdata = indep, type = 'response')
p_sPC_N <- predict(m_sPC_N, newdata = indep, type = 'response')
p_sPC_C <- predict(m_sPC_C, newdata = indep, type = 'response')
p_sPC_PC <- predict(m_sPC_PC, newdata = indep, type = 'response')
p_sPC_sPC<- predict(m_sPC_sPC, newdata = indep, type = 'response')
p_sPC_NC <- predict(m_sPC_NC, newdata = indep, type = 'response')
p_sPC_sNC<- predict(m_sPC_sNC, newdata = indep, type = 'response')
p_NC_P <- predict(m_NC_P, newdata = indep, type = 'response')
p_NC_N <- predict(m_NC_N, newdata = indep, type = 'response')
p_NC_C <- predict(m_NC_C, newdata = indep, type = 'response')
p_NC_PC <- predict(m_NC_PC, newdata = indep, type = 'response')
p_NC_sPC <- predict(m_NC_sPC, newdata = indep, type = 'response')
p_NC_NC <- predict(m_NC_NC, newdata = indep, type = 'response')
p_NC_sNC <- predict(m_NC_sNC, newdata = indep, type = 'response')
p_sNC_P <- predict(m_sNC_P, newdata = indep, type = 'response')
p_sNC_N <- predict(m_sNC_N, newdata = indep, type = 'response')
p_sNC_C <- predict(m_sNC_C, newdata = indep, type = 'response')
p_sNC_PC <- predict(m_sNC_PC, newdata = indep, type = 'response')
p_sNC_sPC<- predict(m_sNC_sPC, newdata = indep, type = 'response')
p_sNC_NC <- predict(m_sNC_NC, newdata = indep, type = 'response')
p_sNC_sNC<- predict(m_sNC_sNC, newdata = indep, type = 'response')
#MODEL ACCURACY USING Mean Absolute Error (MAE)
MAE_indep <- c(accuracy(indep$DOC, p_null[,1])[1,3], accuracy(indep$DOC, p_cntr[,1])[1,3],
accuracy(indep$DOC, p_P_P[,1])[1,3], accuracy(indep$DOC, p_P_N[,1])[1,3],
accuracy(indep$DOC, p_P_C[,1])[1,3], accuracy(indep$DOC, p_P_PC[,1])[1,3],
accuracy(indep$DOC, p_P_sPC[,1])[1,3], accuracy(indep$DOC, p_P_NC[,1])[1,3],
accuracy(indep$DOC, p_P_sNC[,1])[1,3],
accuracy(indep$DOC, p_N_P[,1])[1,3], accuracy(indep$DOC, p_N_N[,1])[1,3],
accuracy(indep$DOC, p_N_C[,1])[1,3], accuracy(indep$DOC, p_N_PC[,1])[1,3],
accuracy(indep$DOC, p_N_sPC[,1])[1,3], accuracy(indep$DOC, p_N_NC[,1])[1,3],
accuracy(indep$DOC, p_N_sNC[,1])[1,3],
accuracy(indep$DOC, p_C_P[,1])[1,3], accuracy(indep$DOC, p_C_N[,1])[1,3],
accuracy(indep$DOC, p_C_C[,1])[1,3], accuracy(indep$DOC, p_C_PC[,1])[1,3],
accuracy(indep$DOC, p_C_sPC[,1])[1,3], accuracy(indep$DOC, p_C_NC[,1])[1,3],
accuracy(indep$DOC, p_C_sNC[,1])[1,3],
accuracy(indep$DOC, p_PC_P[,1])[1,3], accuracy(indep$DOC, p_PC_N[,1])[1,3],
accuracy(indep$DOC, p_PC_C[,1])[1,3], accuracy(indep$DOC, p_PC_PC[,1])[1,3],
accuracy(indep$DOC, p_PC_sPC[,1])[1,3], accuracy(indep$DOC, p_PC_NC[,1])[1,3],
accuracy(indep$DOC, p_PC_sNC[,1])[1,3],
accuracy(indep$DOC, p_sPC_P[,1])[1,3], accuracy(indep$DOC, p_sPC_N[,1])[1,3],
accuracy(indep$DOC, p_sPC_C[,1])[1,3], accuracy(indep$DOC, p_sPC_PC[,1])[1,3],
accuracy(indep$DOC, p_sPC_sPC[,1])[1,3], accuracy(indep$DOC, p_sPC_NC[,1])[1,3],
accuracy(indep$DOC, p_sPC_sNC[,1])[1,3],
accuracy(indep$DOC, p_NC_P[,1])[1,3], accuracy(indep$DOC, p_NC_N[,1])[1,3],
accuracy(indep$DOC, p_NC_C[,1])[1,3], accuracy(indep$DOC, p_NC_PC[,1])[1,3],
accuracy(indep$DOC, p_NC_sPC[,1])[1,3], accuracy(indep$DOC, p_NC_NC[,1])[1,3],
accuracy(indep$DOC, p_NC_sNC[,1])[1,3],
accuracy(indep$DOC, p_sNC_P[,1])[1,3], accuracy(indep$DOC, p_sNC_N[,1])[1,3],
accuracy(indep$DOC, p_sNC_C[,1])[1,3], accuracy(indep$DOC, p_sNC_PC[,1])[1,3],
accuracy(indep$DOC, p_sNC_sPC[,1])[1,3], accuracy(indep$DOC, p_sNC_NC[,1])[1,3],
accuracy(indep$DOC, p_sNC_sNC[,1])[1,3])
m_names <- c('m_null', 'm_cntr',
'm_P_P', 'm_P_N', 'm_P_C', 'm_P_PC', 'm_P_sPC', 'm_P_NC', 'm_P_sNC',
'm_N_P', 'm_N_N', 'm_N_C', 'm_N_PC', 'm_N_sPC', 'm_N_NC', 'm_N_sNC',
'm_C_P', 'm_C_N', 'm_C_C', 'm_C_PC', 'm_C_sPC', 'm_C_NC', 'm_C_sNC',
'm_PC_P', 'm_PC_N', 'm_PC_C', 'm_PC_PC', 'm_PC_sPC', 'm_PC_NC', 'm_PC_sNC',
'm_sPC_P', 'm_sPC_N', 'm_sPC_C', 'm_sPC_PC', 'm_sPC_sPC', 'm_sPC_NC', 'm_sPC_sNC',
'm_NC_P', 'm_NC_N', 'm_NC_C', 'm_NC_PC', 'm_NC_sPC', 'm_NC_NC', 'm_NC_sNC',
'm_sNC_P', 'm_sNC_N', 'm_sNC_C', 'm_sNC_PC', 'm_sNC_sPC', 'm_sNC_NC', 'm_sNC_sNC')
MAE_Diffs <- data.frame(m_names, MAE_indep)
#FIT AND PREDICTION RESULTS
table <- merge(LOO_Diffs, MAE_Diffs, by = "m_names")
write.csv(table, 'results.csv', row.names = FALSE)
#Read brms models already estimated for plots
m_null <- readRDS('m_null.rds')
m_cntr <- readRDS('m_cntr.rds')
m_C_PC <- readRDS('m_C_PC.rds')
m_sPC_C <- readRDS('m_sPC_C.rds')
#Set plot preferences
theme_set(theme_bw(18))
theme_update(axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line(),
legend.title = element_blank())
#FIT PLOTS of best models using train dataset
f_null <- fitted(m_null, type = 'response')
f_cntr <- fitted(m_cntr, type = 'response')
f_C_PC <- fitted(m_C_PC, type = 'response')
f_sPC_C <- fitted(m_sPC_C, type = 'response')
fit <- data.frame(data, f_null, f_cntr, f_C_PC, f_sPC_C)
g_null <- ggplot(data = fit, aes(x = DOC, y = Estimate, fill = 'black')) +
geom_point() + geom_abline(intercept = 0, slope = 1) + theme_update(legend.position = 'none') +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_cntr <- ggplot(data = fit, aes(x = DOC, y = Estimate.1, fill = 'black')) +
geom_point() + geom_abline(intercept = 0, slope = 1) + theme_update(legend.position = 'none') +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_C_PC <- ggplot(data = fit, aes(x = DOC, y = Estimate.2, fill = 'Black')) +
geom_point() + geom_abline(intercept = 0, slope = 1) + theme_update(legend.position = 'none') +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_sPC_C <- ggplot(data = fit, aes(x = DOC, y = Estimate.3, fill = 'Black')) +
geom_point() + geom_abline(intercept = 0, slope = 1) + theme_update(legend.position = 'none') +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_fit <- ggarrange(g_null + rremove('x.text') + rremove('xlab') + rremove('ylab') + rremove('legend'),
g_cntr + rremove('x.text') + rremove('xlab') + rremove('y.text') + rremove('ylab'),
g_C_PC + rremove('xlab') + rremove('ylab'),
g_sPC_C + rremove('xlab') + rremove('y.text') + rremove('ylab'),
ncol = 2, nrow = 2, align = 'hv', hjust = c(-4.4, -4, -4.4, -4),
labels = c('a', 'b', 'c', 'd'))
setEPS()
postscript("Figure 5.eps")
annotate_figure(g_fit, bottom = text_grob('Observed DOC', size = 18), left = text_grob('Predicted DOC', size = 18, rot = 90))
dev.off()
#PREDICTION PLOTS of best models using indep dataset
p_null <- predict(m_null, newdata = indep, type = 'response')
p_cntr <- predict(m_cntr, newdata = indep, type = 'response')
p_C_PC <- predict(m_C_PC, newdata = indep, type = 'response')
p_sPC_C <- predict(m_sPC_C, newdata = indep, type = 'response')
p_indep <- data.frame(indep, p_null, p_cntr, p_C_PC, p_sPC_C)
g_null <- ggplot(data = p_indep, aes(x = DOC, y = Estimate, color = Experiment)) +
geom_point() + geom_abline(intercept = 0, slope = 1) +
scale_color_manual(values=c("Grey", "firebrick1", "firebrick4",
"seagreen1", "seagreen")) +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_cntr <- ggplot(data = p_indep, aes(x = DOC, y = Estimate.1, color = Experiment)) +
geom_point() + geom_abline(intercept = 0, slope = 1) +
scale_color_manual(values=c("Grey", "firebrick1", "firebrick4",
"seagreen1", "seagreen")) +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_C_PC <- ggplot(data = p_indep, aes(x = DOC, y = Estimate.2, color = Experiment)) +
geom_point() + geom_abline(intercept = 0, slope = 1) +
scale_color_manual(values=c("Grey", "firebrick1", "firebrick4",
"seagreen1", "seagreen")) +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_sPC_C <- ggplot(data = p_indep, aes(x = DOC, y = Estimate.3, color = Experiment)) +
geom_point() + geom_abline(intercept = 0, slope = 1) +
scale_color_manual(values=c("Grey", "firebrick1", "firebrick4",
"seagreen1", "seagreen")) +
scale_y_continuous(limits = c(0, 15)) + scale_x_continuous(limits = c(0, 15))
g_indep <- ggarrange(g_null + rremove('x.text') + rremove('xlab') + rremove('ylab'),
g_cntr + rremove('x.text') + rremove('xlab') + rremove('y.text') + rremove('ylab'),
g_C_PC + rremove('xlab') + rremove('ylab'),
g_sPC_C + rremove('xlab') + rremove('y.text') + rremove('ylab'),
ncol = 2, nrow = 2, align = 'hv', hjust = c(-4.4, -4, -4.4, -4),
labels = c('a', 'b', 'c', 'd'), common.legend = TRUE)
setEPS()
postscript("Figure 6.eps")
annotate_figure(g_indep, bottom = text_grob('Observed DOC', size = 18), left = text_grob('Predicted DOC', size = 18, rot = 90))
dev.off()
#Environmental Null model for Table 2
m_envi <- brm(formula = bf(DOC ~ nlab * exp(-nk * time) + nres, #model equation
nlab ~ 0 + Environment, nk ~ 0 + Environment, nres ~ 0 + Environment, nl = TRUE), #parameter dependencies
data = data, family = gaussian(),
prior = c(set_prior("normal(0, 1)", nlpar = "nlab", lb = 0), #prior constraints
set_prior("normal(0, 1)", nlpar = "nk", lb = 0),
set_prior("normal(0, 1)", nlpar = "nres", lb = 0)),
chains = chains, iter = 6000, warmup = 3000, cores = cores, seed = seed,
sample_prior = TRUE, #save parameter sampling information
file = "m_envi", #save model output
control = control) #control chain steps
print(m_envi, digits = 5) #to see posterior summary results
#########