Passing output from other R functions to brms for a non-linear regression

Interfaces brms
,
1

Hi everyone

I have written an R package to help people fit unimodal thermal performance curves to their data using non-linear regression. This includes helper functions that return sensible start parameters, lower, and upper parameter bounds, and functions that contain the equation to fit.

However lots of users (me included) have data in an experimental format where it would make sense to fit non-linear mixed effects models in brms.

The general workflow for the nls() (or other non-linear least squares approach) is:

# load packages
library(rTPC) # remotes::install_github('padpadpadpad/rTPC')
library(brms)
library(tidyverse)

# current approach for using non linear regression 

# filter data
data("bacteria_tpc")
d <- filter(bacteria_tpc, phage == 'nophage')

# get start parameters
start_vals <- get_start_vals(d$temp, d$rate, 'gaussian_1987')

# get upper limits
upper_lims <- get_upper_lims(d$temp, d$rate, 'gaussian_1987')

# get lower limits
lower_lims <- get_lower_lims(d$temp, d$rate, 'gaussian_1987')

# fit model
fit <- nls(rate ~ gaussian_1987(temp = temp, rmax, topt, a),
           data = d,
           start = start_vals,
           lower = lower_lims,
           upper = upper_lims,
           algorithm = 'port')

# get predictions and plot
preds <- broom::augment(fit)

ggplot(preds, aes(temp, rate)) +
  geom_point(size = 3, shape = 21, fill = 'white') +
  geom_line(aes(temp, .fitted)) +
  theme_bw()

Currently to fit it in brms I have to view start_vals, lower_lims, and upper_lims in the console to hardcode priors, and then also view the code of the function in the nls() call to copy and paste it into bf().

# set relatively uninformative priors for each parameter
# using the output from get_start_vals() and upper_lims() and lower_lims()
start_vals
lower_lims
upper_lims

priors <- prior(normal(0.84, 5), nlpar = names(start_vals[1]), lb = 0.28, ub = 8.44) +
  prior(normal(25, 5), nlpar = "topt", lb = 15, ub = 37) +
  prior(normal(22, 5), nlpar = 'a', lb = 0, ub = 220)

# get model formula out of function call
gaussian_1987
# est <- rmax * exp(-0.5 * (abs(temp - topt)/a)^2)

# set formula using bf
brms_form <- bf(rate ~ rmax * exp(-0.5 * (abs(temp - topt)/a)^2),
                rmax + topt + a ~ 1,
                nl = TRUE)

# fit brms model
fit_brms <- brm(brms_form,
           data = d, 
           prior = priors)

Was just wondering if there was a smarter way of doing this and what would everyone’s approach be. I would love to include a method to the helper functions I have which would output code/something else that could be more easily integrated into using brms for these problems. I saw stanvar() but it still feels like there is a lot of hard coding for the user that I would prefer to avoid.