Equation for a random intercept negative binomial model with interaction

Hello everyone, I would like to express a negative binomial model model with an interaction as an equation.
Here is the code used with brms:

mod3 <- brm(formula = y~ x1* x2+ (1 | id), data = data,
            family = negbinomial(), prior = prior(normal(0,5), class = b))

And here’s my attempt to express it as an equation:

y_{i} \sim NegativeBinomial(μ_{i},ϕ),\\ μ_{i} = α_{0}+α_{i} + \beta_{1}X1+\beta_{2}X2+\beta_{3}X1X2\\ α_{0} \sim Normal(0,5)\\ \beta_{1} \sim Normal(0,5)\\ \beta_{2} \sim Normal(0,5)\\ \beta_{3} \sim Normal(0,5)\\ α_{i}∼Normal(0,σ_{i})\\ σ_{i}∼Student \textrm{-}t\textrm{+}(3,0,2.5)\\ ϕ∼Gamma(0.01,0.01)

Does this sound right to you?

Thank you in advance!


Looks mostly right, but I think there will only be one \sigma that is the SD of the \alpha's, rather than each \alpha_i having it’s own \sigma_i. So you could write that for each i

\alpha_i \sim Normal(0, \sigma_\alpha)

You can double check whether the model estimated a posterior distribution of one sigma or many posteriors of many sigmas. But your model specification to me looks like it should just estimate the one, which is fine and what you want for partial pooling of information across the \alpha parameters. I also don’t remember what the default priors on \sigma and \phi are, but what you put could be right.

Actually, I think \alpha_0 may also have a different prior that can be specified using class = "Intercept". So I think currently it would be using the default prior on the intercept rather than the same one specified for the \beta's.

Lately, I’ve been trying to note the non-centered parameterization explicity:

\begin{aligned} y_{i} & \sim NegativeBinomial(μ_{i},ϕ),\\ μ_{i} & = α_{0} + α_{ID,i} + \beta_{1}X1_i+\beta_{2}X2_i+\beta_{3}X1_iX2_i \\ \alpha_0 & \sim \operatorname{Normal}(location,scale) \\ \alpha_{ID,i} & = \sigma_{ID} * \zeta_{ID,i} \\ \sigma_{ID }& \sim \operatorname{Student } \textrm{-}t\textrm{+}(3,0,2.5)\\ \zeta_{ID} & \sim \operatorname{Normal}(0,1) \\ \phi & \sim \operatorname{Gamma}(0.001, 0.001) \\ \beta_1 & \sim \operatorname{Normal}(0, 5) \\ \beta_2 & \sim \operatorname{Normal}(0, 5) \\ \beta_3 & \sim \operatorname{Normal}(0, 5) \\ \end{aligned}

1 Like

[edit: fat-fingered before this was fully edited.]

@andrewgelman has professed liking to write the models this way, too, though I haven’t seen him do it in practice (not that anyone can keep up with his paper output).

Personally, I prefer to just see the model written out naturally (and to put markers in superscripts to save subscripts for indexing and to set them in roman shape because they’re constants). I feel it’s a shortcoming of Stan that we can’t fit the centered parameterizations and the reparameterization always feels like a technical detail when discussing the model.

For example, here I find it easier to follow written this way:

\zeta^\textrm{id} \sim \textrm{normal}(0, \sigma^\textrm{id}).

Now if the goal is to define the space the sampler operates over, then you have to not only do this kind of manual transform, you also have to be careful to transform to an unconstrained parameterization and include the change-of-variables adjustments. Even there, though, I’d be inclined to write the model out in natural form first and then explicitly transform to something easier to sample, for example saying we sample over \zeta^\textrm{id} / \sigma^\textrm{id} rather than over \zeta^\textrm{id} and sample over \log(\sigma^\textrm{id}) rather than over \sigma^\textrm{id}.