- brms Version: 2.10.0
Hi,
This is my first foray into non-normal modelling, so please bear with me if I’m asking obvious questions. I have searched the forum and looked at various documentations, but I remain unclear on a couple of things.
I have data which can be considered long reaction times (between about 20 and 100 seconds). The shape of each participant’s data follows more or a less a shifted lognormal distribution. The design of the experiment is such that there are two within-participant factors with two levels each, resulting in four conditions: A1B1, A2B1, A1B2, A2B2.
I would like to test the effect of condition on the shift parameter and mean of the shifted lognormal. I have defined the model as the following:
formula <- bf(RT ~ 0 + intercept + A + B + (1 | ID),
ndt ~ 0 + intercept + A + B + (1 | ID))
priors <- c(set_prior("normal(0,20)", class = "b", coef = "intercept"),
set_prior("normal(0,20)", class = "b", coef = "groupunfair"),
set_prior("normal(0,20)", class = "b", coef = "partnerunfair"),
set_prior("cauchy(0,1)", class = "sd"),
set_prior("normal(50,20)", class = "b", coef = "intercept", dpar = "ndt"),
set_prior("normal(0,50)", class = "b", coef = "groupunfair", dpar = "ndt"),
set_prior("normal(0,50)", class = "b", coef = "partnerunfair", dpar = "ndt")
)
model <- brm(formula, family = shifted_lognormal(),
prior = priors,
data = dat, warmup = 2000,
iter = 10000, chains = 3,
control = list(adapt_delta = 0.97),
cores = 4, inits = "0",
sample_prior = TRUE,
save_all_pars = TRUE)
The model seemed to converge, with well-mixed chains etc, and the posterior predictive check fits the data surprisingly well. However, I have some questions about the defining the priors and understanding the output.
-
As you can see, I defined wide priors on the RT scale. However, in this post, Paul seems to suggest that priors of the lognormal should be defined on the log scale. When I tried to do that, I received an error message stating that no samples could be drawn (annoyingly I seem to have deleted the code where I defined these priors - sorry, bad practice, I know).
-
Here is the output of the model (I’ve altered the predictor names to match the definitions above for ease):
Family: shifted_lognormal
Links: mu = identity; sigma = identity; ndt = log
Formula: LT ~ 0 + intercept + A + B + (1 | ID)
ndt ~ 0 + intercept + A + B + (1 | ID)
Data: dat (Number of observations: 1296)
Samples: 3 chains, each with iter = 10000; warmup = 2000; thin = 1;
total post-warmup samples = 24000
Group-Level Effects:
~ID (Number of levels: 25)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.67 0.10 0.50 0.90 1.00 4570 8184
sd(ndt_Intercept) 2.67 0.74 1.63 4.48 1.00 5916 8951
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
intercept 3.38 0.14 3.10 3.66 1.00 3057 6018
A2 -0.05 0.06 -0.16 0.06 1.00 11896 15970
B2 -0.16 0.05 -0.26 -0.07 1.00 14393 15952
ndt_intercept 0.13 0.74 -1.58 1.34 1.00 4574 7373
ndt_A2 -0.31 0.07 -0.44 -0.16 1.00 12563 13073
ndt_B2 0.47 0.11 0.25 0.66 1.00 10483 13871
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.63 0.02 0.60 0.67 1.00 14286 17161
My understanding was that the estimates should be on the scale of the outcome variable (seconds), yet clearly here this is not the case (the mean RT in the data is 34s). Am I right in assuming the estimates here are in log(RT)? If so, this makes sense for the estimate of the intercept (exp(3.38 = 29.37), but I’m confused about rest of the estimates. When I run a mixed model in lme4::lmer, the estimates are about a 3 and 9 second difference for the different parameters. I know this a different kettle of fish completely given the fact it’s a different model and approach, but these tiny effects are baffling me a bit.
I fully expect that it is to do with a complete non-understanding of shifted lognormal parameters or how to define the model/priors correctly, so would really appreciate a hand in understanding what’s going on.
Thanks in advance,
Anthony
