RStudio__brms_play.R

# testing non-linear regression methods

# devtools::install_github('padpadpadpad/nlsLoop')

library(nlsLoop)
library(dplyr)
library(tidyr)
library(ggplot2)
library(brms)
library(ggridges)

data("Chlorella_TRC")

# make column for unique things apart from flux
Chlorella_TRC <- unite(Chlorella_TRC, id2, c(process, growth.temp, flux, curve_id), sep = ':', remove = F)

# filter for P and R
P <- filter(Chlorella_TRC, flux == 'photosynthesis')
R <- filter(Chlorella_TRC, flux == 'respiration')

# run nlsLoop
fits <- nlsLoop::nlsLoop(ln.rate ~ schoolfield.high(ln.c, Ea, Eh, Th, temp = K, Tc = 20),
data = Chlorella_TRC,
tries = 500,
id_col = 'id2',
param_bds = c(-30, 30, 0, 10, 0, 30, 245, 385),
r2 = 'Y',
supp.errors = 'Y',
AICc = 'Y',
na.action = na.omit,
lower = c(ln.c=-10, Ea=0, Eh=0, Th=0))

# run nlsLoop2
fits_2 <- nlsLoop::nlsLoop2(ln.rate ~ schoolfield.high(ln.c, Ea, Eh, Th, temp = K, Tc = 20),
data = Chlorella_TRC,
tries = 1000,
id_col = 'id2',
param_bds = c(-10, 10, 0.1, 2, 0.5, 5, 285, 330),
alg = 'plinear-random',
na.action = na.omit,
lower = c(ln.c=-10, Ea=0, Eh=0, Th=0))

# try brms - bayesian non-linear model fitting ####
Chlorella_TRC <- mutate(Chlorella_TRC, lnrate = ln.rate,
curve_id2 = ifelse(flux == 'photosynthesis', curve_id - 15, curve_id),
flux2 = ifelse(flux == 'respiration', 'aresp', flux))

ggplot(Chlorella_TRC, aes(K, lnrate, col = flux)) +
geom_point() +
facet_wrap(~curve_id2)

# test that functionin br() call creates the correct curve
K = 290 : 325; lnc = 1; Ea = 0.6; Eh = 4; Th = 315
rate <- lnc + log(exp(Ea/(8.62*10^-5) * (1/293.15 - 1/K))) + log(1/(1 + exp(Eh/(8.62*10^-5) * (1/Th - 1/K))))
plot(rate ~ K)

# make schoolfield high model into a non-linear brms model
nlform <- bf(lnrate ~ lnc + log(exp(Ea/(8.62*10^-5) * (1/293.15 - 1/K))) + log(1/(1 + exp(Eh/(8.62*10^-5) * (1/Th - 1/K)))),
Ea ~ 0 + flux + (1|curve_id),
lnc ~ 0 + flux + (1|curve_id),
Eh ~ 0 + flux + (1|curve_id),
Th ~ 0 + flux + (1|curve_id),
nl = TRUE)

# set up priors
nlprior <- c(prior(normal(0, 10), nlpar = 'lnc'),
prior(normal(0.6, 10), nlpar = 'Ea', lb = 0),
prior(normal(4, 10), nlpar = 'Eh', lb = 0),
prior(normal(300, 10), nlpar = 'Th', lb = 273.15))

# run model
fit_bayes1 <- brm(formula = nlform,
data = Chlorella_TRC,
family = gaussian(),
prior = nlprior,
control = list(adapt_delta = 0.9),
chains = 3,
iter = 2e3)

# get summary
summary(fit_bayes1)

# plot traceplots
plot(fit_bayes1)

# create dataframe for predictions over
conditions <- data.frame(curve_id = unique(Chlorella_TRC$curve_id))
row.names(conditions) <- unique(Chlorella_TRC$curve_id)

# calculate marginal effects
me_fit <- marginal_effects(fit_bayes1, conditions = conditions, re_formula = NULL, method = 'predict')

# plot marginal effects
plot(me_fit, ncol = 10, points = TRUE)
marginal_effects(fit_bayes1)

# extract posteriors and wrangle
posteriors <- posterior_samples(fit_bayes1, add_chain = TRUE) %>%
select(., starts_with('b'), chain) %>%
gather(., key = 'param', value = 'estimate', starts_with('b')) %>%
separate(., 'param', c('mean', 'param', 'flux'), sep = '_') %>%
mutate(., flux = ifelse(grepl('resp', flux), 'respiration', 'photosynthesis'))

# plot posteriors
ggplot(posteriors, aes(x = estimate)) +
geom_line(aes(group = interaction(chain, flux), col = flux), stat = 'density') +
facet_wrap(~ param, scale = 'free') +
scale_color_manual(values = c('green4', 'black')) +
scale_x_continuous(expand=c(0, 0))