How can I get the parameter values which define a given confidence interval?

Modeling
1

Hello! I need help on how to get the results I’m interested in from a brms model I fit.
Once again, my background in statistics is lackluster, so pardon me if I make some mistake while writing this post.

The data I’m working with is structured as follows:

> head(dat_teste)
  sample phase TS  gamma  tau
1 A02010     2 65 69.940 3033
2 A02010     2 65 38.470 2861
3 A02010     2 65 20.480 2726
4 A02010     2 65 10.870 2611

Where phase and TS are numerical discrete grouping factors, gamma and tau are my (x,y) data.

Here’s the model I fit using brms.

modelo_HB <- brms::bf(tau ~ a + b * (gamma^c),
                      a ~ 1 + (1|TS) + (1|TS:phase),
                      b ~ 1 + (1|TS) + (1|TS:phase),
                      c ~ 1 + (1|TS) + (1|TS:phase),
                      sigma ~ gamma + (1|TS),
                      nl = TRUE)

prior <- c(prior(uniform(10, 8000), class = 'b', nlpar = 'a', lb = 10, ub = 8000),
           prior(uniform(1, 1000), class = 'b', nlpar = 'b', lb = 1, ub = 1000),
           prior(uniform(10^-3, 1), class = 'b', nlpar = 'c', lb = 10^-3, ub = 1),
           prior(student_t(3, 0, 2.5), class = 'Intercept', dpar = 'sigma'))

As you can see, I’m trying to fit the following nonlinear equation to my data:

\tau = a + b\gamma^c

With a, b and c as the model parameters. The fit went decently well.

What I’m specially interested in is defining the credible intervals for the fit, for a given value of TS, and unknown value of phase, and defining the boundaries of the credible interval region in the (\gamma, \tau) plane using curves of the \tau = a + b\gamma^c format.

Simply put, for a given value of TS (a.k.a. for a given group) and a non specified value of phase, I want to have two (a,b,c) parameter vectors (one for the bottom limit, one for the upper limit) which define the boundaries of the credible interval region using the equation \tau = a + b\gamma^c.

How could I do this? Any help is greatly appreciated.

2

Howdy! Before addressing your question, I want to point out that uniform priors with hard upper and lower boundaries like those you specified for the a, b, and c parameters are generally not recommended. You might need some sort of bound, for example, if the parameter can only be positive, but in general it is best to set priors that don’t have these hard boundaries. The prior is a piece of information in your model - here, you are saying that parameter a has equal probability of being any value between 10 and 8,000, but it has zero probability of being 9.99999999 or 8,000.00000001. That kind of information seems highly unrealistic for most scenarios. Since you have decided on these boundaries, that tells me that you do know something about the parameters, unless the boundaries are completely arbitrary (in which case you likely would not have chosen them). Since you know something about the parameters, then it would be best to reflect this in some sort of distribution that places more probability density around the region that is more probable with tails that stretch into areas of low probability. For example, maybe you know that a should be somewhere around 4,000 but you think it could reasonably be between 10 and 8000. Then you might put a normal(4000, 2000) prior on it (and if it can’t be less than zero, then you add the lower bound of zero). One really revealing way to view the implications of your priors is by doing prior predictive checks. This can be done in brms by running your model with the sample_prior = "only" option, which will ignore the likelihood. You can then use the various pp_check() commands and options to view the implications of your prior on the scale of the outcome variable. I think that you will see that the information from your uniform priors imply some rather strange things about the model.

To answer your question - you can do this graphically with the conditional_effects() command in brms (see the documentation for extensive details and options, particularly as regards your group-level effects). You can also do it yourself manually by creating a data frame of the values of \gamma and TS that you want to make predictions on and either using the fitted() or the predict() functions in brms, with the appropriate re_formula = options, to obtain the credible intervals for \tau in the (\gamma, \tau) plane using the a, b, and c parameters from your model. See the documentation for fitted() and predict() to understand the difference between them, so that you can decide which it is that you are looking for. Based on your question it sounds like you are wanting fitted().

I hope that helps!

1 Like