Simple non-linear multilevel model gone awry


#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.