(sorry in advance for my extensive message… i’m a beginner in stan and BRMS)
Here’s my problem.
I have a dependent variable, denoted REP, binary (1 or 0). I wish to predict the probability of answering 1 with a non-linear model such as,
REP ~ inv_logit ((X/(1 + k*D)) - X2). with X and X2 two numerical variables , D a numerical variable and k : the free parameter to estimate.
I would like to estimate the parameter “k” according to the different experimental conditions I have realized.
Independent variable A: “DOMAIN” with two modalities (1: LOSS, 2: GAIN). Each subject is exposed to these 2 modalities of this independent variable.
Independent variable B: “STATUS” with two modalities (1: Control, 2: Threat). A participant is either “control” or “threat” but never both.
my data frame :
I would like to estimate the parameter “k” for each condition :
k1_Control_GAIN < k3_Threat_GAIN
k2_Control_LOSS < k4_Threat_LOSS
I have several ways to code this I think, the way I tried is this: (I’ve simplified the model here to make it more understandable)
main_model = bf( rep ~ inv_logit( ((X / (1 + exp(k1 * dummy_ctrl * dummy_gain + k2 * dummy_gain * dummy_threat + k3 * dummy_loss * dummy_ctrl + k4 * dummy_loss * dummy_threat)*D ) - X2)) , k1 ~ ( 1| ID ), k2 ~ ( 1| ID ) , k3 ~ ( 1| ID ), k4 ~ ( 1| ID ), nl = TRUE)
dummy_ctrl = 1 when the subject is in the “control” group and 0 otherwise
dummy_threat = 1 when the subject is in the “threat” group and 0 otherwise
dummy_GAIN= 1 when the subject is exposed to a trial of the “GAIN” condition and 0 otherwise
with dummy_LOSS= 1 when the subject is exposed to a trial of the “LOSS” condition and 0 otherwise
main_fit = brm(main_model, data = DF, family = bernoulli(link='identity'), prior = c(prior(normal(-5,3), nlpar = 'k1'), prior(normal(-5,3), nlpar = 'k2'), prior(normal(-5,3), nlpar = 'k3'), prior(normal(-5,3), nlpar = 'k4')), cores = parallel::detectCores(), seed = "12345", inits = "random", chains = 4, iter = 4000, warmup = 2000, control = list(adapt_delta=0.99))
I don’t think my way is right.
I think it’s better to do it like this?
main_model = bf( rep ~ inv_logit( ( (X / (1 + exp(k)*D ) - X2)) , k ~ DOMAIN*STATUS + (1 + DOMAIN | ID), nl = TRUE)
Thank you in advance for your responses.