Non-linear model fitting problem

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:

  1. Are there some guidelines on how to approach the choice of priors in this case?
  2. I guess this is an identifiability issue? So, should I just try to reparameterize the model in a different way?
  3. 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.

I’m hardly expert on this issue given that I’m stuck with my own non-linear model, but I see two improvements that will get yours to fit.

1) Strengthen the prior on "lambda"
\mathcal{N}(-1,10) is pretty weak, and as I understand it non-linear models need extra attention to the prior specification in order to converge. How you go about this might be a bit more challenging. It’s clear how to proceed if, say, lambda = 10 is biologically unreasonable in your study system. But if those extreme values are plausible, others would have to weigh in.

2) Rescale model parameters
From another post, @paul.buerkner recommends:

Putting the two together, I come to:

fit_lactin <- brm(
  bf(growthrate ~ exp(rho*temperature) - exp(rho*TmaxSC*10 - (TmaxSC*10-temperature)/delta) + lambdaSC/10,
     rho ~ 1,
     TmaxSC ~ 1,
     delta ~ 1,
     lambdaSC ~ 1,
     nl = TRUE),
  data = ryandata, family = gaussian(),
  prior = c(
    prior(normal(0.1,10), nlpar = "rho"),
    prior(normal(4.2,1.0), nlpar = "TmaxSC", lb = 0.1, ub = 5.0),
    prior(normal(7,10), nlpar = "delta"),
    prior(normal(-10,2), nlpar = "lambdaSC")
  ),
  warmup = 5000, iter = 10000, chains = 2, cores = 2,
  control = list(adapt_delta = 0.95)
)

You can see that I’ve just multiplied Tmax and divided lambda—both in the model formula and in the priors—to bring them onto a similar order of magnitude as the other parameters. I append “SC” as a suffix to remind myself that the posteriors for these parameters need to be back-transformed. Although it’s a bit harder to see because of this scaling, I’ve also dramatically narrowed the prior on lambda. This fits reasonably quickly and gives the following:

Population-Level Effects: 
                   Estimate Est.Error l-95% CI
rho_Intercept          0.08      0.01     0.07
TmaxSC_Intercept       4.38      0.01     4.35
delta_Intercept        7.07      0.56     5.94
lambdaSC_Intercept    -9.52      0.62   -10.52
                   u-95% CI Eff.Sample Rhat
rho_Intercept          0.09       1902 1.00
TmaxSC_Intercept       4.40       3831 1.00
delta_Intercept        8.03       2010 1.00
lambdaSC_Intercept    -8.01       2593 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI
sigma     0.17      0.02     0.14     0.22
      Eff.Sample Rhat
sigma       4260 1.00

These correspond (after backtransforming the “SC” parameters) quite nicely to the estimates from nls:

Parameters:
        Estimate Std. Error t value Pr(>|t|)    
rho     0.073298   0.006612  11.086 2.57e-13 ***
Tmax   43.789369   0.148651 294.579  < 2e-16 ***
delta   6.928933   0.858966   8.067 1.13e-09 ***
lambda -0.980015   0.044927 -21.813  < 2e-16 ***

I hope this helps in your case, but so far these “tricks” haven’t gotten me out of my issue. I would also be keen to learn about the strategies for re-parameterization. I also haven’t been able to find specific, practical suggestions on this front. Perhaps these strategies could help you make use of weaker priors that you may have to use when fitting other data to this model.

2 Likes

Hi, thanks a lot for the advice. Rescaling seems to solve the problem in this case. Good to know.