Simple non-linear multilevel model gone awry

Interfaces brms
1

Hi,

Below is a reproducible example.

Issue: fitting non-linear, non-multilevel model yields acceptable results; however, fitting it on a multi-level model either doesn’t converge or converge but Rhat >>> 1.

Looking forward to insights as to why the above happens.

  library(brms)
  data <- data.frame(
    well = c(
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
      2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4,
      4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5
    ),
    t = c(
      0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8,
      9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6,
      7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
    ),
    q = c(
      600, 494.9847987, 414.0056588, 350.4265778, 299.720482,
      258.7206902, 225.1629665, 197.3966393, 174.1973846,
      154.6429526, 138.0288875, 1454.062171, 801.5207162, 585.898834,
      472.9998825, 401.9515008, 352.4708526, 315.7108008, 287.1471172,
      264.2064972, 245.3087319, 229.4255148, 1273.315867, 785.8348879,
      609.819227, 512.8691917, 449.616207, 404.3074031, 369.8600656,
      342.5634412, 320.2637214, 301.6150225, 285.7279857, 1283.14929,
      805.6565995, 630.1897515, 532.7537012, 468.8465601, 422.8895521,
      387.8406599, 359.9957651, 337.1981136, 318.0964517, 301.7958765,
      1394.391596, 855.1369184, 667.835811, 565.2573645, 498.3358037,
      450.3221429, 413.7394178, 384.6832399, 360.8906989, 340.9483724,
      323.9221504
    )
  )

  #non-multilevel model
  model <- bf(q ~ qi * (1 + D * b * t)^(-1 / b),
    qi ~ 1,
    D ~ 1,
    b ~ 1,
    nl = TRUE
  )

  priors <- c(
    prior(normal(0, 2), nlpar = "D", lb = 0),
    prior(normal(1, 1), nlpar = "b", lb = 0.0001),
    prior(normal(800, 500), nlpar = "qi", lb = 0),
    prior(student_t(3, 0, 20), class = "sigma")
  )

  fit <- brm(model,
    data = data, family = gaussian(), prior = priors,
    chains = 2, iter = 1000,
    warmup = 500
  )
  
  (prior <- brms::get_prior(model_multilevel,
                            data = data, family = gaussian()))
  
  summary(fit)
  
  #multilevel model
  model_multilevel <- bf(q ~ qi * (1 + D * b * t)^(-1 / b),
    qi ~ 1 + (1 | well),
    D ~ 1 + (1 | well),
    b ~ 1 + (1 | well),
    nl = TRUE
  )

  fit_multilevel <- brm(model_multilevel,
    data = data, family = gaussian(), prior = priors,
    chains = 4, iter = 10000,
    warmup = 500,
    control = list(adapt_delta = 0.95, max_treedepth = 25)
  )
  
  summary(fit_multilevel)

Please also provide the following information in addition to your question:

  • R version 3.5.1 (2018-07-02)
    Platform: x86_64-w64-mingw32/x64 (64-bit)
    Running under: Windows 10 x64 (build 17134)
  • brms Version:2.6.0
2

The boundaries on the priors you set will only affect the overall intercept not the varying intercepts per well. Thus, you need to use other ways to range restrict certain parameters. For instance, if you want D to be positive, drop the lb argument from prior and replace D with exp(D) in the non-linear formula.

3

Thanks for the feedback @paul.buerkner.

I’ll try your suggestion and update this post.

If I understand you correctly, brms currently doesn’t support specifying priors to group-level, nor does the population-level priors trickle down.

4

It’s not does currently not support it’s because of the generality of the non-linear framework there is no general way to impose boundaries via priors if the linear predictor contains more than one parameter.

5

@paul.buerkner: tried your suggestion; no luck, still very low Eff.Sample and high Rhat.

Also, wouldn’t substituting D for exp(D) in the non-linear formula completely change the formula definition?

I tried using lognormal() as family, using the same formula as per original post, and model converged.

6

Changing D to exp(D) changes the scale on which D is defined but makes sure that whatever the result is, it will be positive. This is likely one necessary step to solve your problem.

In general, non-linear models are hard and require carefully chosen priors and a lot of fine-tuning to get running reliably.

7

Noted; thanks for the advice @paul.buerkner.