Using initial values in RStan

I am developing a hierarchical model that fits a non-linear regression to individuals nested within species. One of the parameters in the regression controls the convexity of the curve. If I estimate a single convexity value shared across all individuals and all species, the model runs and converges nicely. If I try to estimate convexity values for each individual (Since it’s just a tuning parameter, I don’t bother with nesting them within species), the model fails to initialize.

It’s not a model specification problem, because if I specify inits for this parameter (based loosely on the posterior mean and standard deviation of the parameter estimated from analysis of individual curves), the model runs and converges nicely. To generate the inits, though, I have to use a bit of a hack.

Specifically, I am working with RStan, and I set up an environment to hold the number of individuals:

## define new environment
init <- new.env()
## set n_indiv to default (illegal) value
n_indiv <- -999.0
## return to parent environment

Then just before the call to stan() when I know the number of individuals in the data set I set n_indiv:

## set n_indiv for use in init_theta
init$n_indiv <- nrow(unique_table)

That all works just fine, but it would be a lot cleaner if I could just specify

init = init_theta(n_indiv)

(or something similar) in the call to stan(). Am I missing something?

The init argument can be a function that returns a list of (partial) initial values

So, this could be something like

init_theta <- function(chain_id) {
  list(theta = ...)

It gets called in the environment of the data list (I think).

I don’t think I asked my question clearly. That would let me return a list of values from init_theta(). My problem is that I want to feed a data dependent value (n_indiv) into init_theta(). In outline

n_indiv <- 10
fit <- stan(model, data, pars, init=init_theta(n_indiv))

It doesn’t appear that I can do that. Am I missing something?