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), 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.