Init not using my initial values and seems to be defaulting to 0

Interfaces brms
,
1

I am having difficulties setting up init for fitting a nonlinear dataset. Here is a synthetic dataset.

library(tidyverse)
library(brms)
library(ggmcmc)

set.seed(1)

# synthetic Roche data
Time <- c(0, 2.1, 4.5, 6.5)

#replicates at each Time level
rep = 1:50

# set up first order curve

#curve parameters: Response ~ top - ((top-bottom)*exp(-rate * Time)
# "top" is max response, "bottom" is initial response and "rate" is the 
# first order rate constant in min-1

top = 14
bottom = 8
rate = 0.22

SynthData <- expand_grid(Time, rep) %>%
  mutate(Response = top - (top-bottom)*exp(-rate*Time)) %>%
  rowwise() %>%
  mutate(Response = Response + rnorm(1,0,0.5)) %>%
  ungroup()

Visualization so you can see this is a shallow first order curve with just a few timepoints but many reps.

ggplot(SynthData, aes(x=Time, y= Response))+
  geom_point()+
  theme_bw()

Now for the default run in brm, this does 4 chains:

# nonlinear brm fitting default settings
prior1 <- prior(normal(14, 0.1), nlpar = "top") + 
  prior(normal(8.0, 0.1), nlpar = "bottom")+         
  prior(normal(0.2, 0.05), nlpar = "rate")           

# default fitting
fit1 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Time)), top + bottom + rate ~ 1, nl = TRUE),
            data = SynthData, 
            prior = prior1,
)

summary(fit1)

All looks good, surprisingly good, but, hey I am just a new refugee from frequentist land.
Inspecting the movement of the arguments through the iterations just for the argument Top for the first order curve:

ggs(fit1, burnin = TRUE) %>%
  filter(Parameter == "b_top_Intercept") %>% 
  mutate(chain = factor(Chain),
         intercept = value) %>% 
  
  ggplot(aes(x = Iteration, y = intercept, color = chain)) +
  geom_line()+
  ggtitle("default settings")+
  theme_bw()

Like I said, in this case the default setting (apparently 0) found the global minimum, but I have been seeing similar datasets having maybe one of the chains stuck on a local minimum so it would be nice to specify initial values for the nonlinear model arguments. This is redoing the brm run for just one chain and using “init” to specify the starting values.

### just one chain and supplying initial arg values via init

prior1 <- prior(normal(14, 0.1), nlpar = "top") + #note estimates have to be quite close
  prior(normal(8.0, 0.1), nlpar = "bottom")+         #to get an accurate minimization
  prior(normal(0.2, 0.05), nlpar = "rate")          #sd had to be drasticaly reduced 


fit2 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Time)), top + bottom + rate ~ 1, nl = TRUE),
            data = SynthData, 
            prior = prior1,
            chains = 1,
            cores = 4,
            seed = 1,
            init = list(list ( top  =  14 ,
                               bottom  =  8 ,
                               rate  =  0.22))
)

# visualization
ggs(fit2, burnin = TRUE) %>%
  filter(Parameter == "b_top_Intercept") %>% 
  mutate(chain = factor(Chain),
         intercept = value) %>% 
  
  ggplot(aes(x = Iteration, y = intercept, color = chain)) +
  geom_line()+
  ggtitle("supplying top = 14")+
  theme_bw()

# changing initial conditions to top=5

fit3 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Time)), top + bottom + rate ~ 1, nl = TRUE),
            data = SynthData, 
            prior = prior1,
            chains = 1,
            cores = 4,
            seed = 1,
            init = list(list ( top  =  5 ,
                               bottom  =  8 ,
                               rate  =  0.22))
)


# visualization
ggs(fit3, burnin = TRUE) %>%
  filter(Parameter == "b_top_Intercept") %>% 
  mutate(chain = factor(Chain),
         intercept = value) %>% 
  
  ggplot(aes(x = Iteration, y = intercept, color = chain)) +
  geom_line()+
  ggtitle("supplying initial value top = 5")+
  theme_bw()

it seems that my initial values specified in the init setting are being ignored. I gathered from online vignettes (eg Solomon’s excellent article Solomon) that this should be specified as a list of lists. Do I have something wrong in the syntax? or do I have a more fundamental misunderstanding?

Thanks in advance,
Thomas

  • Operating System: “Windows” “10 x64” “build 26100” “US-LP-10681” “x86-64”
  • R version: 4.4.2 (2024-10-31 ucrt)
  • brms Version: 2.22.0
2

OK I found the problem. From this thread it seems that the underlying parameter names in the stan code must be used in the R script when setting up init parameters. The command stancode() can be used to inspect the underlying stan code.

In my case this can be observed in a parameter section:

parameters {
  vector[K_top] b_top;  // regression coefficients
  vector[K_bottom] b_bottom;  // regression coefficients
  vector[K_rate] b_rate;  // regression coefficients
  real<lower=0> sigma;  // dispersion parameter
}

indicating that init needs to be coded like this?:

fit2 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Time)), top + bottom + rate ~ 1, nl = TRUE),
            data = SynthData, 
            prior = prior1,
            chains = 1,
            cores = 4,
            seed = 1,
            init = list(list ( b_top  =  14 ,
                               b_bottom  =  8 ,
                               b_rate  =  0.22))
)

unfortunately now this happens:
SAMPLING FOR MODEL ‘anon_model’ NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b_top; dims declared=(1); dims found=() (in ‘string’, line 20, column 2 to column 22).

scary.

3

More progress. Here a similar problem is found with dimension errors and the recommendation was to wrap parameter values in an array.

fit2 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Time)), top + bottom + rate ~ 1, nl = TRUE),
            data = SynthData, 
            prior = prior1,
            chains = 1,
            cores = 4,
            seed = 1,
            init = function(){list(b_top=array(14),
                                   b_bottom=array(8), 
                                   b_rate=array(0.22))}
)

And this worked! The parameters are starting out at 14, 8, 0.22. Problem is that I dont really understand why this is working now. Anyone have any simple explanations and general rules I could apply to prevent these sorts of errors in the future?

Thanks
Thomas

4

Welcome to the Stan forum and sorry for the hassle. If I recall correctly the array issue comes up only when dealing with scalar values. Depending on how the underlying Stan code for this brms model is written Stan may be expecting an object with 1 dimension whereas in R a single number is dimensionless and wrapping it in array gives it a single dimension (compare dim(1) to dim(array(1)) in R).

In the future if you see something like dims declared=(1); dims found=() this should solve the problem.