Working with brms version 2.16.3 on windows 10.
I am trying to specify the following model in brms:
Decision ~ Normal(u,sd)
u = baseline + social information
baseline = 1 + Score
social information = (Behaviour * (1+Success)) * Social
Behaviour is a variable with 3 levels (Increase, Neutral or Decrease) and Success is a dummy variable 0/1. Success being 0 represents the Majority condition. Social is a dummy variable where 1 indicates that participants saw social information. IE. I need the effect of social information to be 0 if they did not see any social information. I use non linear syntax 0 + Behaviour because I want to index the behaviour variable, rather than dummy code it.
I have attempted to build the model using the following definition.
Model <- brm(data = d,
family = gaussian,
bf(Decision ~ a + (b * (1 + Success) * Social),
a ~ 1 + Score,
b ~ 0 + Behaviour,
nl = TRUE),
prior = c(prior(normal(0,4), nlpar = a),
prior(normal(0,4), nlpar = b),
prior(exponential(1), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 222)
Where stancode generates this:
// generated with brms 2.16.3
functions {
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=1> K_a; // number of population-level effects
matrix[N, K_a] X_a; // population-level design matrix
int<lower=1> K_b; // number of population-level effects
matrix[N, K_b] X_b; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
}
parameters {
vector[K_a] b_a; // population-level effects
vector[K_b] b_b; // population-level effects
real<lower=0> sigma; // dispersion parameter
}
transformed parameters {
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] nlp_a = X_a * b_a;
// initialize linear predictor term
vector[N] nlp_b = X_b * b_b;
// initialize non-linear predictor term
vector[N] mu;
for (n in 1:N) {
// compute non-linear predictor values
mu[n] = nlp_a[n] + nlp_b[n];
}
target += normal_lpdf(Y | mu, sigma);
}
// priors including constants
target += normal_lpdf(b_a | 0, 4);
target += normal_lpdf(b_b | 0, 4);
target += exponential_lpdf(sigma | 1);
}
generated quantities {
}
The model runs without any issues. However, the model predictions change depending which of Success or Majority I include in the first line. ie.
bf(Decision ~ a + (b * (1 + Success) * Social),
Produces different results to:
bf(Decision ~ a + (b * (1 + Majority) * Social),
Which shouldn’t happen, and makes me think I am making some mistake in my specification of the model.
My questions are:
- Is there an obvious mistake in the way I am attempting to build this model?
- Is there an alternative way I can build a model that is equivalent to the definition I provide above?
I hope this is enough information for someone to offer some advice as I am unfortunately not able to just share the data in this case.
I really appreciate any help anyone can offer.