Different hyperparameters for random effects mean

[ericb-ambevtech]

Hello there,

In light of the discussion #808 , I would like to know if it's possible to further choose the hyperparameters for each individual mean in a mixed effects model. For example, in the Pigs dataset, if I wanted to influence the random intercept of the first two pigs, but in different directions:

1|Pig ~ Normal(mu: [-1, +1, 0, ..., 0], sigma: HalfNormal(sigma: 133.0346)).

I'm aware that this may violate the theory behind random effects modeling in the sense that the estimates are draws of a distribution. But still, is there a way to "tweak" the model specification to set this influence while having the pooling behavior from the mixed effects model?

Thanks!

[zwelitunyiswa]

I think you would need to assign fixed effects for each pig and then assign hyper-parameters. So think like weight ~ 1 + pig then pig ~ Normal(mu: [-1, 1, 0, ..., 0]), , sigma: HalfNormal(sigma: 133.0346))

[tomicapretto]

@ericb-ambevtech, you can do the following:

import bambi as bmb
import numpy as np
import statsmodels.api as sm

data = sm.datasets.get_rdataset("dietox", "geepack").data

n_pigs = data["Pig"].nunique() # 72
mu_vector = np.array([-1, 1] + [0] * (n_pigs - 2)) # Subtract 2 since you set it for two pigs

priors = {
    "1|Pig": bmb.Prior("Normal", mu=mu_vector, sigma=bmb.Prior("HalfNormal", sigma=1))
}
model = bmb.Model("Weight ~ Time + (Time|Pig)", data, priors=priors)
model

       Formula: Weight ~ Time + (Time|Pig)
        Family: gaussian
          Link: mu = identity
  Observations: 861
        Priors: 
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: 60.7258, sigma: 133.0346)
            Time ~ Normal(mu: 0.0, sigma: 18.1283)
        
        Group-level effects
            1|Pig ~ Normal(mu: [-1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
              0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
              0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
              0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.], sigma:
                  HalfNormal(sigma: 1.0))
            Time|Pig ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 18.1283))
        
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 24.9644)

However, make sure you do know what are the first two groups (pigs in this case).

You can explore it like this

model.components["mu"].design.group

GroupEffectsMatrix with shape (861, 144)
Terms:  
  1|Pig  
    kind: intercept
    groups: ['4601', '4602', '4603', '4605', '4641', '4643', '4645', '4756', '4757', '4759', '4760',
      '4813', '4814', '4815', '4817', '4854', '4856', '4857', '4858', '5389', '5392', '5497', '5500',
      '5501', '5502', '5524', '5527', '5528', '5578', '5581', '5582', '5850', '5851', '5852', '5862',
      '5865', '5866', '6055', '6056', '6057', '6058', '6207', '6208', '6211', '6284', '6287', '6288',
      '6430', '6432', '6433', '6909', '6910', '6912', '8049', '8050', '8051', '8053', '8139', '8141',
      '8142', '8144', '8191', '8192', '8193', '8195', '8269', '8270', '8271', '8273', '8437', '8439',
      '8442']
    columns: 0:72
  Time|Pig  
    kind: numeric
    groups: ['4601', '4602', '4603', '4605', '4641', '4643', '4645', '4756', '4757', '4759', '4760',
      '4813', '4814', '4815', '4817', '4854', '4856', '4857', '4858', '5389', '5392', '5497', '5500',
      '5501', '5502', '5524', '5527', '5528', '5578', '5581', '5582', '5850', '5851', '5852', '5862',
      '5865', '5866', '6055', '6056', '6057', '6058', '6207', '6208', '6211', '6284', '6287', '6288',
      '6430', '6432', '6433', '6909', '6910', '6912', '8049', '8050', '8051', '8053', '8139', '8141',
      '8142', '8144', '8191', '8192', '8193', '8195', '8269', '8270', '8271', '8273', '8437', '8439',
      '8442']
    columns: 72:144

To access the actual design matrix do 'np.array(this_obj)'

And you can see they're order from smallest to largest... but anyway, always double-check it