Beta distribution: "Rejecting initial value" with proper initial values

Modeling
1

Hello,

I am fitting a non-linear model of the form

y(x) = p_2 + (p_4 - p_2) * e^{-e^{p_1} * (x - p_3) * g((x - p_3) * 10)}

where g is the inverse logit function. This is being fitted to a beta-distributed data (proportion):

data_str <- "y,x
0.838313611,0.09162907
0.879521184,0.09162907
0.850763463,0.09162907
0.865867666,0.09162907
0.881867882,0.09162907
0.849611448,0.09162907
0.999000000,0.91629073
0.920652405,0.91629073
0.919116384,0.91629073
0.947916772,1.60943791
0.948684782,1.60943791
0.893691199,1.60943791
0.890315996,2.30258509
0.926028477,2.30258509
0.916428348,2.30258509
0.662600928,2.70805020
0.628424467,2.70805020
0.590023950,2.70805020
0.446406015,3.21887582
0.366148933,3.21887582
0.335064725,3.21887582
0.201794719,3.55534806
0.287427873,3.55534806
0.298685288,3.55534806
0.076329756,3.91202301
0.114032082,3.91202301
0.104257405,3.91202301
0.007632976,4.60517019
0.007632976,4.60517019
0.007632976,4.60517019"

data <- read.table(text = data_str, sep = ",", header = TRUE)

I am trying to fit it in brms using my own initial values, but I am getting lots of the classic rejection “Rejecting initial value” sampling statement:

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: beta_lpdf: First shape parameter[1] is -0.071096, but must be > 0!  (in 'model4936186f04f5_844d16a56a3a8ac20fc7d8c9ff39d375' at line 47)

Line 47 in the generated stan code is:

target += beta_lpdf(Y | mu * phi, (1 - mu) * phi);

this means that mu < 1. I have checked for that though, and I am 100% positive that my initial values yield mu values that conform with the beta distribution expectation of shape parameters > 0 and < 1.

my_inits <- list(
  list(b_p1 = -0.07, b_p2 = 0.11, b_p3 = 4.52, b_p4 = 0.96, phi = 0.67),
  list(b_p1 = 0.19, b_p2 = 0.09, b_p3 = 3.36, b_p4 = 0.97, phi = 0.07),
  list(b_p1 = 1.50, b_p2 = 0.07, b_p3 = 3.58, b_p4 = 0.93, phi = 0.30),
  list(b_p1 = 0.49, b_p2 = 0.30, b_p3 = 3.83, b_p4 = 0.97, phi = 0.64)
)

mu_formula <- function(p1, p2, p3, p4, x) {
  p2 + (p4 - p2) * exp(-exp(p1) * (x - p3) *
    plogis((x - p3) * 10))
}

for (i in seq_along(my_inits)) {
  print(range(mu_formula(my_inits[[i]]$b_p1[1], my_inits[[i]]$b_p2[1],
                         my_inits[[i]]$b_p3[1], my_inits[[i]]$b_p4[1],
                         data$x)))
}
[1] 0.9139801 0.9611010
[1] 0.2852354 0.9999382
[1] 0.07869372 0.97270772
[1] 0.4891421 0.9883545

Yet, brms cannot fit the model. I do not quite understand what could be going on. Does anyone have an idea? Shouldn’t those values ensure that the very first calculation of mu is valid? Or are the initial values intended to do different thing?

Full reprex:

library(brms)

priors <- prior("normal(0, 1", nlpar = "p1") +
          prior("beta(1, 5)", nlpar = "p2") +
          prior("gamma(5, 1)", nlpar = "p3") +
          prior("beta(5, 1)", nlpar = "p4") +
          prior("gamma(0.1, 0.1)", class = "phi")

bform <- bf(y ~ p2 + (p4 - p2) * exp(-exp(p1) *
              (x - p3) * inv_logit((x - p3) * 10)),
            p2 + p4 + p1 + p3 ~ 1,
            nl = TRUE)

test <- brm(bform, data = data, family = Beta(link = "identity"),
            prior = priors, iter = 2e4, warmup = 1.6e4,
            chains = 4, cores = 4, inits = my_inits)
2

I never used brms, but is there a fixed_param option you can use to make sure your model is really doing what it’s supposed to? I also don’t use R so although the model is pretty simple I may be missing something obvious in your check (and and probably am not the best person to reply here, but there were none so far), like is Stan inv_logit the same as R plogis?.

The inits are not too many, so you can try not giving it any and see if it works, but it at this point you know the model is doing something differently from what you expect, so maybe you need some Rubber duck debugging.

1 Like
3

When running your code, I get the warning message

1: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior 
b_p2 ~ beta(1, 5)
b_p3 ~ gamma(5, 1)
b_p4 ~ beta(5, 1)
 
2: It appears as if you have specified an upper bounded prior on a parameter that has no natural upper bound.
If this is really what you want, please specify argument 'ub' of 'set_prior' appropriately.
Warning occurred for prior 
b_p2 ~ beta(1, 5)
b_p4 ~ beta(5, 1)
 
3: In .local(object, ...) :
  some chains had errors; consider specifying chains = 1 to debug

When I change the priors to

priors <- prior("normal(0, 1", nlpar = "p1") +
  prior("beta(1, 5)", nlpar = "p2",lb = 0, ub = 1) +
  prior("gamma(5, 1)", nlpar = "p3", lb = 0) +
  prior("beta(5, 1)", nlpar = "p4", lb = 0, ub = 1) +
  prior("gamma(0.1, 0.1)", class = "phi")

I get the following error (which I don’t understand unfortunately)

Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b_p2; dims declared=(1); dims found=()
[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                              
[2] "  mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b_p2; dims declared=(1); dims found=()"

However, the model samples normally if I don’t specify any inits.
What was the reason you wanted to specify inits? Are they important for your model or did you try to solve a problem that is now fixed with the upper and lower bounds?

Additionally, I seem to get different errors than you, so what brms version do you use?
I use brms 2.14.4.

2 Likes
4

I played around with your code some more, and I think I found the problem.
When assigning your my_inits, leave out the b_ in front of the parameter names.
Unfortunatly, I also had to learn the hard way that the program does not alert you when you set initial values for parameters the model doesn’t recognise. With these inits, the model samples normally for me (regardless of whether there are bounds on the priors set with ub and lb):

my_inits <- list(
  list(p1 = -0.07, p2 = 0.11, p3 = 4.52, p4 = 0.96, phi = 0.67),
  list(p1 = 0.19, p2 = 0.09, p3 = 3.36, p4 = 0.97, phi = 0.07),
  list(p1 = 1.50, p2 = 0.07, p3 = 3.58, p4 = 0.93, phi = 0.30),
  list(p1 = 0.49, p2 = 0.30, p3 = 3.83, p4 = 0.97, phi = 0.64)
)
1 Like
5

Hi @caesoma and @BrittaL,

Thanks for replying to my post, and apologies for taking so long to respond, I was away on holidays.
I realised a couple of things that were wrong with my initial code, which only became clearer after Britta’s investigation.

The b_ prefix should in fact be part of the named list; your case worked but brms chose its own initial values instead with the exception of the parameter phi which had the correct name specification (i.e. without b_).
e.g. dropping b_ from the list names makes brms spit out the following list of initial values (it fitted the model as it did in your case, but also returned some rejection statements):

> test$fit@inits
[[1]]
[[1]]$b_p2
[1] 0.2203274

[[1]]$b_p4
[1] 0.5453935

[[1]]$b_p1
[1] 1.446889

[[1]]$b_p3
[1] 0.2394709

[[1]]$phi
[1] 0.67


[[2]]
[[2]]$b_p2
[1] 0.5225934

[[2]]$b_p4
[1] 0.5426361

[[2]]$b_p1
[1] 0.1323563

[[2]]$b_p3
[1] 0.3963176

[[2]]$phi
[1] 0.07


[[3]]
[[3]]$b_p2
[1] 0.5020893

[[3]]$b_p4
[1] 0.01993598

[[3]]$b_p1
[1] -1.023359

[[3]]$b_p3
[1] 0.9610738

[[3]]$phi
[1] 0.3


[[4]]
[[4]]$b_p2
[1] 0.6035035

[[4]]$b_p4
[1] 0.2988594

[[4]]$b_p1
[1] 0.9129861

[[4]]$b_p3
[1] 1.538402

[[4]]$phi
[1] 0.64

One of the error messages you got made me realise that I needed to specify dimensions to my initial values as well.

my_inits <- list(
  list(b_p1 = -0.07, b_p2 = 0.11, b_p3 = 4.52, b_p4 = 0.96, phi = 0.67),
  list(b_p1 = 0.19, b_p2 = 0.09, b_p3 = 3.36, b_p4 = 0.97, phi = 0.07),
  list(b_p1 = 1.50, b_p2 = 0.07, b_p3 = 3.58, b_p4 = 0.93, phi = 0.30),
  list(b_p1 = 0.49, b_p2 = 0.30, b_p3 = 3.83, b_p4 = 0.97, phi = 0.64)
)

for (i in seq_along(my_inits)) {
  for (j in 1:4) {
    dim(my_inits[[i]][[j]]) <- 1
    print(dim(my_inits[[i]][[j]]))
  }
}

Note that I did not specify a dimension for phi because that is specified as a real in the stan code:

real<lower=0> phi;  // precision parameter

But the other parameters are defined as vectors of dimension 1:

  vector[K_p2] b_p2;  // population-level effects
  vector[K_p4] b_p4;  // population-level effects
  vector[K_p1] b_p1;  // population-level effects
  vector[K_p3] b_p3;  // population-level effects

I also added the bounds to the priors to avoid the warning messages at the end (and note that my original code had a typo on the p1 prior as well):

priors <- prior("normal(0, 1)", nlpar = "p1") +
          prior("beta(1, 5)", nlpar = "p2", lb = 0, ub = 1) +
          prior("gamma(5, 1)", nlpar = "p3", lb = 0) +
          prior("beta(5, 1)", nlpar = "p4", lb = 0, ub = 1) +
          prior("gamma(0.1, 0.1)", class = "phi")

Implementing the changes below eliminated all sampling rejection statements (meaning all initial values worked as intended, I think), and the returned list of initial values was identical to my_inits:

> test$fit@inits
[[1]]
[[1]]$b_p2
[1] 0.11

[[1]]$b_p4
[1] 0.96

[[1]]$b_p1
[1] -0.07

[[1]]$b_p3
[1] 4.52

[[1]]$phi
[1] 0.67


[[2]]
[[2]]$b_p2
[1] 0.09

[[2]]$b_p4
[1] 0.97

[[2]]$b_p1
[1] 0.19

[[2]]$b_p3
[1] 3.36

[[2]]$phi
[1] 0.07


[[3]]
[[3]]$b_p2
[1] 0.07

[[3]]$b_p4
[1] 0.93

[[3]]$b_p1
[1] 1.5

[[3]]$b_p3
[1] 3.58

[[3]]$phi
[1] 0.3


[[4]]
[[4]]$b_p2
[1] 0.3

[[4]]$b_p4
[1] 0.97

[[4]]$b_p1
[1] 0.49

[[4]]$b_p3
[1] 3.83

[[4]]$phi
[1] 0.64

Thanks for investigating this, the error messages you posted here provided clues to the solution!