Hi there,
up until version 2.1.0, the main parameter modelled with family = exgaussian was the mean of the Gaussian component. Since version 2.3.0, this has changed to model the mean of the distribution. Is there a way of modelling the mean on the Gaussian component in version 2.3.0 or newer?
Many thanks!
You can implement it via the custom_family feature. See vignette("brms_customfamilies").
And I recommend you update brms to version 2.4.0 from CRAN.
If something doesnāt work out, please feel free to ask for advice!
Dear Paul,
I am a bit confused about how the exGaussian model works in brms v2.7.0.
In the vignette " Parameterization of Response Distributions in brms" you say the following about the parameters:
āĪ² is the scale (inverse rate) of the exponential component, Ī¾ is the mean of the Gaussian component, Ļ is the standard deviation of the Gaussian component, and erfc is the complementary error function. We parameterize Ī¼=Ī¾+Ī² so that the main predictor term equals the mean of the distribution.ā
From this and the post here I understand that the point estimates in the model reflect the sum of the mean of the Gaussian component and the scale of the exponential component. Thus, the estimates are estimates of Ī¼ and not of Ī¾? -> Is this an accurate interpretation?
If the estimates represent Ī¾+Ī², then Iām a bit confused on how to interpret the way the formula is written in the case of the exGaussian model.
If I write the formula: āRT ~ predictorā
ā> does this mean that I state that the mean of the Gaussian component (Ī¾) depends on the predictor or am I stating that the sum of the Gaussian mean and scale of the exponential component (Ī¼) depends on the predictor? Since I get 1 point estimate of the scale of the exponential component (Ī²) in my summary, I assume that the formula means that only the Gaussian mean (Ī¾) depends on the predictor?
If only the mean of the Gaussian component (Ī¾) is modeled as a function of the predictor when the formula is written as āRT ~ predictorā and not as ābf(RT~predictor, beta~predictor)ā, why canāt the estimates of the mean of the Gaussian component (Ī¾) be derived from subtracting the predicted value of Ī² from the point estimate (Ī¼)?
Iām sorry if this is a stupid question, but itās all very new to meā¦
Thanks!
Hanne
Yes, in recent version of brms, the main formula predicts mu and not xi anymore. This used to be different but I changed it for consistency with most other families where we predict the mean via the main formula. You may of course derive the predictions for xi by computing mu - beta.
Thank you for the quick response!
Given this ānewā parameterization, am I correct in my understanding that the model structure for Ī¼ and Ī² imply a structure for Ī¾?
That is, if I have RT ~ a, beta ~ 1, this implies a structure such that mu ~ a (i.e. we express that the mean can vary by a but that the exponential component explicitly does not, implying that the gaussian mean component must be permitted to vary by a).
Similarly, if we have RT ~ 1, beta ~ a then by necessity Ī¾ will be perfectly anti-correlated with Ī²? This latter strikes me as troublesome for sampling. Arguably this is only a problem in the case where the user allows Ī² to vary by a given predictor but does not allow the overall mean to also vary by that predictor, but I can imagine scenarios where users building up increasingly complex models might not catch on to this implied structure. Possibly a warning is warranted?
You are right with your conclusion that mu and beta imply a strucutre for xi as is the case when we predict xi and beta which then imply a structure for mu.
I am not sure I can follow why a warning should be warranted though. I would need to see some (simulated or real) examples that show this pathodological behavior.
As a follow up on this post, I have been trying to set up a custom_family function that would provide the effects on the mean of the Gaussian component xi (instead of the default mean of the ex-Gaussian component). Hereās my attempt:
exgaussian_mg <- custom_family(
"exgaussian_mg",
dpars = c("mu", "sigma", "beta"),
links = c("identity", "log", "log"),
lb = c(NA, 0, 0),
type = "real"
)
stan_funs <- "
real exgaussian_mg_lpdf(real y, real mu, real sigma, real beta) {
// This is the log-probability density function (lpdf)
real lambda = 1.0 / beta;
real val = lambda * (y - mu);
return log(lambda) + // this is the rate of the function
normal_lpdf(y | mu, sigma) + // this is the gaussian lpdf (built-in)
log_diff_exp(0, normal_lcdf(0 | val, lambda * sigma)); // exponential component
}"
stanvars <- stanvar(scode = stan_funs, block = "functions")
This does seem to work in the following call:
formulas <- bf(RT ~ condition*primetype*experiment + (1||target) + (1||subj),
beta ~ condition*primetype*experiment + (1||target) + (1||subj))
# modelling
bay_fit <- brm(
formula = formulas,
data = data_sample,
family = exgaussian_mg,
prior = c(
set_prior('normal(0, 100)', class = "b") # for fixed effects only
),
stanvar = stanvars,
cores = getOption("mc.cores", 4),
threads = threading(2), #within-chain parallelization
chains = 4, # Number of MCMC chains to run at the same time
iter = 2000, # Total iterations per chain (including warmup)
warmup = 1000, # Number of warmup iterations (discarded for posterior inference)
seed = 42, # Random seed for reproducibility of MCMC sampling
init = 0, # set the initial values for the sampler; the inits_list won't work for some reason
control = list(adapt_delta = 0.99, max_treedepth = 15), # Increased adapt_delta and max_treedepth to prevent divergences
)
However, the chains seem to loop indefinitely even after reach the number of iterations:
Compiling Stan program...
Start sampling
starting worker pid=19488 on localhost:11205 at 16:57:39.837
starting worker pid=19501 on localhost:11205 at 16:57:40.097
starting worker pid=19514 on localhost:11205 at 16:57:40.350
starting worker pid=19527 on localhost:11205 at 16:57:40.591
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.020153 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 201.53 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.020182 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 201.82 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.014891 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 148.91 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4:
Chain 4: Gradient evaluation took 0.018561 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 185.61 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 2602.1 seconds (Warm-up)
Chain 1: 1371.14 seconds (Sampling)
Chain 1: 3973.24 seconds (Total)
Chain 1:
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 5014.59 seconds (Warm-up)
Chain 2: 1023.65 seconds (Sampling)
Chain 2: 6038.24 seconds (Total)
Chain 2:
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 9248.64 seconds (Warm-up)
Chain 4: 3010.83 seconds (Sampling)
Chain 4: 12259.5 seconds (Total)
Chain 4:
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 15129 seconds (Warm-up)
Chain 2: 1304.98 seconds (Sampling)
Chain 2: 16434 seconds (Total)
Chain 2:
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 117649 seconds (Warm-up)
Chain 4: 1081.54 seconds (Sampling)
Chain 4: 118731 seconds (Total)
Chain 4:
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 131888 seconds (Warm-up)
Chain 3: 508.108 seconds (Sampling)
Chain 3: 132396 seconds (Total)
Chain 3:
# manually stopped here after ~3 days
So I wanted to ask for advice here:
Thanks!
Can anybody help?
Many thanks in advance!
@scholz might be able to help you :)
Hey, I currently donāt have time to review your custom code but could you try if one of our custom families from bayesfam work for you to see if itās your general setup or the custom family code?