A quick question to hopefully address a bit of confusion regarding the specification of the init argument (in particular for Wiener models). My model is giving me trouble because of impossible values for the bias parameter.
The formula is:
formula <- bf(rt | dec(response) ~ condition,
bias ~ condition)
family <- wiener(link_bs = "identity",
link_ndt = "identity",
link_bias = "identity")
- I started by creating a function to generate a list of init values, naively set to their posterior variable names (as obtained with
names(as.data.frame(m))):
init_func <- function(chain_id=1) {
list(bs = 1, ndt = 0.1,
b_condition = 0, b_Intercept = 0,
b_Intercept_bias = 0.5, b_bias = 0)
}
Supplying this function to the init argument didn’t seem to improve (note: I have fairly narrow priors on the parameter so I don’t think that is the source of the issue).
- I then read more carefully the documentation of
init, that states that:
If specifying initial values using a list or a function then currently the parameter names must correspond to the names used in the generated Stan code (not the names usedi n R).
I ran stancode() on the model, and tried to locate there the initial values (note that I’m not familiar with the stan language). If I’m not mistaken (?), it could correspond to this paragraph:
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
real<lower=0> bs; // boundary separation parameter
real<lower=0,upper=min_Y> ndt; // non-decision time parameter
vector<lower=0,upper=1>[Kc_bias] b_bias; // population-level effects
real<lower=0,upper=1> Intercept_bias; // temporary intercept for centered predictors
}
As we can see, while the posterior variable is called for instance b_Intercept_bias, it is named Intercept_bias in the stan code. So I changed the correspond name in the init function, but it didn’t seem to help (note that no warning was displayed in all cases).
- I stumbled then onto this similar thread, that deals with the
ndtparameter. Centrally, Paul says to initialize inits likelist(temp_ndt_Intercept = -3).
Question: what are the correct names of the init values? Should we add temp_ before them even though it doesn’t appear in the stan code?
On a developer side, it would be really helpful to throw some warnings if the supplied init values do not correspond to any parameter, or are wrongly specified. Or to be able to easily access the init names and their default values.
Here’s all the stan code for reference in case I am looking at the wrong bit:
> stancode(m)
// generated with brms 2.17.4
functions {
/* Wiener diffusion log-PDF for a single response
* Args:
* y: reaction time data
* dec: decision data (0 or 1)
* alpha: boundary separation parameter > 0
* tau: non-decision time parameter > 0
* beta: initial bias parameter in [0, 1]
* delta: drift rate parameter
* Returns:
* a scalar to be added to the log posterior
*/
real wiener_diffusion_lpdf(real y, int dec, real alpha,
real tau, real beta, real delta) {
if (dec == 1) {
return wiener_lpdf(y | alpha, tau, beta, delta);
} else {
return wiener_lpdf(y | alpha, tau, 1 - beta, - delta);
}
}
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=0,upper=1> dec[N]; // decisions
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> K_bias; // number of population-level effects
matrix[N, K_bias] X_bias; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
real min_Y = min(Y);
int Kc = K - 1;
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
int Kc_bias = K_bias - 1;
matrix[N, Kc_bias] Xc_bias; // centered version of X_bias without an intercept
vector[Kc_bias] means_X_bias; // column means of X_bias before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
for (i in 2:K_bias) {
means_X_bias[i - 1] = mean(X_bias[, i]);
Xc_bias[, i - 1] = X_bias[, i] - means_X_bias[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
real<lower=0> bs; // boundary separation parameter
real<lower=0,upper=min_Y> ndt; // non-decision time parameter
vector<lower=0,upper=1>[Kc_bias] b_bias; // population-level effects
real<lower=0,upper=1> Intercept_bias; // temporary intercept for centered predictors
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += student_t_lpdf(b | 2, 0, 2);
lprior += student_t_lpdf(Intercept | 2, 0, 2);
lprior += gamma_lpdf(bs | 2, 1.5);
lprior += gamma_lpdf(ndt | 1.5, 3)
- 1 * gamma_lcdf(min_Y | 1.5, 3);
lprior += normal_lpdf(b_bias | 0.5, 0.05)
- 1 * log_diff_exp(normal_lcdf(1 | 0.5, 0.05), normal_lcdf(0 | 0.5, 0.05));
lprior += beta_lpdf(Intercept_bias | 2, 2);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
// initialize linear predictor term
vector[N] bias = Intercept_bias + Xc_bias * b_bias;
for (n in 1:N) {
target += wiener_diffusion_lpdf(Y[n] | dec[n], bs, ndt, bias[n], mu[n]);
}
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// actual population-level intercept
real b_bias_Intercept = Intercept_bias - dot_product(means_X_bias, b_bias);
}