Please also provide the following information in addition to your question:
- Operating System: Red Hat Enterprise Linux
- brms Version: version 2.8.0 (package)
Hi,
I’m starting to use brms for non-linear modeling. I’m trying to fit a non-linear model, modeling the growth rate of a fungus in response to temperature. Example data from Ryan et al. 1943. Am. J. Bot. 30: 784-799
ryandata <- data.frame(temperature = c(42.30, 41.0, 39.0, 37.0, 35.7, 34.80, 34.35, 33.50, 32.80, 31.30, 29.80, 29.10, 28.50, 27.20, 26.10, 25.30, 24.80, 23.40, 22.70, 21.30, 20.0, 19.70, 18.0, 16.50, 15.0, 14.50, 13.10, 11.60, 11.40, 10.50, 10.20, 9.20, 8.30, 7.50, 6.70, 6.10, 5.80, 5.50, 5.10, 4.50, 3.80), growthrate = c(1.30, 2.50, 4.0, 4.90, 5.20, 5.20, 5.20, 5.0, 5.20, 4.60, 4.40, 4.30, 4.30, 3.80, 3.70, 3.70, 3.50, 3.30, 3.10, 2.90, 2.70, 2.70, 2.30, 2.10, 1.90, 1.70, 1.60, 1.30, 1.20, 1.10, 1.0, 0.86, 0.72, 0.54, 0.37, 0.31, 0.36, 0.30, 0.18, 0.07, 0.06))
I’m fitting a lactin model to the data, which is motivated by some theory and biological considerations:
$$ y = \exp{\rho T} - \exp{\rho T_{max} - \frac{T_{max} - T}{\delta}} + \lambda $$
I can use standard nls in R to fit this model:
start.val <- c(rho = 0.05, Tmax = 43, delta = 8, lambda = -1)
model.lactin <- nls(growthrate ~ exp(rho*temperature) - exp(rho*Tmax - (Tmax-temperature)/delta) + lambda, start = start.val, trace = T, data = ryandata)
summary(model.lactin)
Formula: growthrate ~ exp(rho * temperature) - exp(rho * Tmax - (Tmax -
temperature)/delta) + lambda
Parameters:
Estimate Std. Error t value Pr(>|t|)
rho 0.073298 0.006612 11.086 2.57e-13 ***
Tmax 43.789369 0.148651 294.578 < 2e-16 ***
delta 6.928934 0.858967 8.067 1.13e-09 ***
lambda -0.980015 0.044927 -21.813 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1657 on 37 degrees of freedom
Number of iterations to convergence: 18
Achieved convergence tolerance: 7.335e-06
Plotting reveals that while the model is not ideal, it is not completely unreasonable:
plot(ryandata$temperature, ryandata$growthrate)
lines(ryandata$temperature, predict(model.lactin))
I then try to use brms to fit the model:
fit_lactin <- brm(
bf(growthrate ~ exp(rho*temperature) - exp(rho*Tmax - (Tmax-temperature)/delta) + lambda,
rho ~ 1, Tmax ~ 1 , delta ~ 1 , lambda ~ 1 , nl = TRUE),
data = ryandata, family = gaussian(),
prior = c(
prior(normal(0.1,10), nlpar = "rho"),
prior(normal(42,10), nlpar = "Tmax", lb = 1, ub = 50),
prior(normal(7,10), nlpar = "delta"),
prior(normal(-1,10), nlpar = "lambda")
),
warmup = 5000, iter = 10000, chains = 2, cores = 2, control = list(adapt_delta = 0.95)
)
However, so far I’ve not managed to get the model to convergence. The chains seem to get stuck in different values. I’ve read brms non-linear vignette that states that non-linear model are difficult to fit in general. I’ve tried different priors but no luck so far. So my questions are:
- Are there some guidelines on how to approach the choice of priors in this case?
- I guess this is an identifiability issue? So, should I just try to reparameterize the model in a different way?
- Any advice for reparameterization, on how to go about doing this?
So, if anyone has any suggestions or could point me towards useful resources that would be great.