# How to write a multilevel GLM poisson model in mathematical notation?

I have constructed a Poisson GLM (modeling it Bayesian in brms) as a part of a social network analysis, predicting counts of outgroup/ingroup connections over two different cohorts/years of a study programme. We are investigating whether the COVID-lockdown has had any effect on how many outgroup connections a student has in general, and we have some social network collected data for this on two different cohorts (years) of the same study, one affected by the lockdown, and one unaffected.

The model contains an interaction and varying effects, and I have a hard time grasping how to write up the model and the priors in mathematical notation.

My likelihood is, as mentioned, Poisson.

Our model is like this - in R-syntax:

log ( λ_i) = Count_of_Connections ~ 0 + Year:Group + (0 + Group | ID)

And the variables are:

• Count_of_Connections discrete count of friendships
• Year: Categorical variable of 2 levels, as we have two years 2018 (unaffected by lockdown) and 2019 (affected by lockdown)
• Group: Categorical variable of 2 levels, as we have Ingroup and Outgroup, which defines whether the Connection counted is one made to a person WITHIN the assigned studygroup the student has been tied to throughout the semester, or if its a Connection made to a person OUTSIDE the assigned studygroup (more likely when there is not a lockdown, because you meet people - outside the ones you normally see - randomly at physical lectures, parties, etc.)
• ID: A personal ID for every person

Question 1: How do I write up an interaction and varying effects (i.e. this model) in mathematical notation? Google does help a little, but it seems there are many ways to do this sort of thing and not an agreed upon one way.

If interested - I have defined the following priors (vague) currently to let the data speak (here in R-syntax):

priors_for_model_main ← c(

prior(normal(1.7, 0.5), class = b, coef = year2018),

prior(normal(1.7, 0.5), class = b, coef = year2019),

prior(normal(0.5, 0.2), class = b, coef = year2018:groupoutgroup),

prior(normal(0.5, 0.2), class = b, coef = year2019:groupoutgroup),

prior(lkj(5), class = cor),

prior(normal(0, 0.1), class = sd))

Question 2 for the interested: how would one write up these priors?

There is no problems in getting the model to run - this post is only for me to learn more about how to mathematically notate this. Appreciating all help given and apologies for the R-syntax - I really want to learn how to do this, as my course in Bayesian Statistics does not have the mathematical notation as a core focus but rather the application of the methodology.

• Why do you have the interaction Year:Group but neither of the main effects in the fixed part of your model ? – Robert Long May 22 at 18:59
• We are essentially conducting a paired t-test between the 2018-Out, 2019-Out, 2018-In and 2019-In in like a 2 by 2 design - but we have also got the model including the main effect of Year, and the model with only Year:Group outperformed the other in terms of LOO-IC. Is there some modelling/conceptual issue with only including the interaction? – Prams May 22 at 19:40
• You've cross-posted here and on the Stan forum: discourse.mc-stan.org/t/… . I don't know if cross-posting is officially discouraged or not, but please at least provide cross-links between cross-posts if you're going to do it, so that people don't waste effort answering a question that has already been resolved in another place. – Ben Bolker May 22 at 19:57
• "Is there some modelling/conceptual issue with only including the interaction?" Yes, see here for example: stats.stackexchange.com/questions/11009/… – Robert Long May 22 at 20:05
• Thanks for the link and info! I have already made the model log ( λ_i) = Count ~ 0 + Year + Year:Group + (0 + Group | ID) and its cousin log (λ_i) = Count ~ 0 + Year + Group + Year:Group + (0 + Group | ID), so the former with 1 main, the latter with both. Both perform worse than the interaction-only model in terms of LOO. Other than that, they show the same effects as this interaction-only model does. So while I agree it does make sense to include at least Year as main, the effects are unchanged, just shown in a different way - and the model is worse - so interaction-only model seemed better. – Prams May 22 at 20:57

I was going to recommend the equatiomatic package to you, but it looks like this model is a little bit beyond its capabilities at present. I would write this model as

$$\begin{split} y_i & \sim \textrm{Poisson}(\eta_i) \\ \eta_i & = \beta_{\textrm{year}(i):\textrm{group}(i)} + b_{\textrm{indiv}(i), \textrm{group}(i)} \\ \mathbf b_{\textrm{indiv}(i)} & \sim \textrm{MVN}(\mathbf 0, \mathbf\Sigma) \end{split}$$

where $$i$$ is an observation index; $$\textrm{year}(i)$$, $$\textrm{group}(i)$$, $$\textrm{indiv(i)}$$ refer to the identifying values associated with observation $$i$$; $$\mathbf\Sigma$$ is a 2 $$\times$$ 2 general positive definite matrix.

You could also write this out in a more abstract way ($$\eta_i = \mathbf X \boldsymbol \beta + \mathbf Z \mathbf b$$) (compactness at the cost of transparency), or in terms of indicator variables ($$\boldsymbol \eta = \beta_1 \cdot \textrm{I}_{\textrm{year}=2018, \textrm{group}=\textrm{in}} + \ldots$$)

As far as the distinction between all the different ways of writing out the fixed-effects model (~ Year * Group vs. 0 + Year:Group vs Year:Group vs ...); as long as the formula includes the interaction (i.e. has a : or a * in it somewhere ...), these models will all give identical likelihoods; the model matrices ($$\mathbf X$$) are all convertible via linear transformations. In principle, then, you could also specify priors that would make all of these models identical, but you'd have to specify joint priors on the $$\beta$$ parameters (e.g. a multivariate normal prior on the $$\boldsymbol \beta$$ vector). The default treatment contrasts in R (i.e. setting up the model as ~Year*Group generally make it easier to test comparisons, but make it harder to specify separate priors in a sensible way.

Because the likelihoods are all equivalent, the fact that models specified in different ways are giving different predictive accuracies must (logically) be caused by your priors.

• I'm a little puzzled by your coef specifications, which look like they go with a model specification of 0+Year+Year:Group rather than 0+Year*Group (compare colnames(model.matrix(~ 0+Year:Group, dd)) with colnames(model.matrix(~ 0+Year+Year:Group, dd))
• in terms of 95% prior intervals, the priors you've specified suggest a range of about 2-15 counts for ingroup connections for both years and a range of 10% - 250% increase in counts between ingroup and outgroup in both years. This is neutral in terms of differences across years, if that's what you intended ... but it does strongly suggest that there will be more outgroup than ingroup contacts (prior prob of outgroup<ingroup = pnorm(-0.5/0.2) = 0.006).