Hi Martin, first of all thanks a lot for your reply. I tried to implement your solution but without much success. I will try to describe my model a bit better so that you can have a clearer idea of what I’m trying to achieve and maybe provide some guidance.

The model I’m trying to build is one that should be able to predict the daily energy consumption timeseries of an individual building, using as predictors the two variables I mentioned in the main post, and pooling on a variable that is representative of the shape of the electricity load of the building on a certain day, and that I call here ‘cluster’ since it is detected with a clustering algorithm.

I am quite unexperienced with multilevel modeling, so I started building my models by following the radon notebook, gradually increasing the complexity. If we wanted to create a parallel with the radon, we could say that:

- the 85 radon readings would be replaced by 365 electricity consumption readings for the year I’m trying to predict,
- the floor_idx predictor variable would be in my case the two temperatures (heating and cooling),
- the different Counties would be the ‘clusters’ detected, that represent the load shape of each of the 365 days I’m trying to predict (in this case there are 7 different clusters that repeat throughout the year)

The last model I managed to compute and make work, is one where both the intercept and the temperature coefficients are varying according to the cluster variable:

```
with pm.Model(coords=coords) as varying_intercept_and_temp:
cluster_idx = pm.Data("cluster_idx", cluster, dims="obs_id")
heating_temp = pm.Data("heating_temp", temp_dep_h, dims="obs_id")
cooling_temp = pm.Data("cooling_temp", temp_dep_c, dims="obs_id")
# Hyperpriors:
a = pm.Normal("a", mu=0.0, sigma=10.0)
sigma_a = pm.Exponential("sigma_a", 1.0)
bh = pm.Normal("bh", mu=0.0, sigma=1.0)
sigma_bh = pm.Exponential("sigma_bh", 0.5)
bc = pm.Normal("bc", mu=0.0, sigma=1.0)
sigma_bc = pm.Exponential("sigma_bc", 0.5)
# Varying intercepts:
a_cluster = pm.Normal("a_cluster", mu=a, sigma=sigma_a, dims="Cluster")
# Varying slopes:
bh_cluster = pm.Normal("bh_cluster", mu=bh, sigma=sigma_bh, dims="Cluster")
bc_cluster = pm.Normal("bc_cluster", mu=bc, sigma=sigma_bc, dims="Cluster")
# Electricity prediction
mu = a_cluster[cluster_idx] + bh_cluster[cluster_idx] * heating_temp + bc_cluster[cluster_idx] * cooling_temp
# Model error:
sigma = pm.Exponential("sigma", 1.0)
y = pm.Normal("y", mu, sigma=sigma, observed=log_electricity, dims="obs_id")
```

Now what I would like to know is if there is a way (and if would make sense) to create a new model where the intercept and the 2 slopes covary, in a similar way to the covariation model described in the radon notebook, but now with one intercept and 2 slopes, instead of just one intercept and one predictor variable.

Hope I managed to describe my problem clearly enough (not really easy to sum up in a forum post, as it is the subject of a scientific article, but I tried my best).

Thanks again for your help!