Obtaining predictions from multilevel model trace

Hi all, I’m trying to run a multilevel model to estimate an electricity consumption time series for a building, but I’m having issues in generating the predictions once I obtained the trace.

My goal is to obtain hourly electricity consumption values depending on the moment of the day (daypart) and on the outdoor temperature, so I formulated this partial pooling model where I have an intercept and a temperature slope which are pooled among the different parts of the day.

with pm.Model(coords=coords) as partial_pooling:

daypart_idx = pm.Data(“daypart_idx”, daypart, dims=“obs_id”)
temperature = pm.Data(“temperature”, outdoor_temp, dims=“obs_id”)

#Hyperpriors:
a = pm.Normal(“a”, mu=0.0, sigma=10.0)
sigma_a = pm.Exponential(“sigma_a”, 1.0)
bt = pm.Normal(“bt”, mu=0.0, sigma=1.0)
sigma_bt = pm.Exponential(“sigma_bt”, 0.5)

#Varying intercepts:
a_daypart = pm.Normal(“a_cluster”, mu=a, sigma=sigma_a, dims=“Daypart”)
#Varying slopes:
bt_daypart = pm.Normal(“bt_cluster”, mu=bt, sigma=sigma_bt, dims=“Daypart”)

#Electricity prediction
mu = a_daypart[daypart_idx] + bt_daypart[daypart_idx] * temperature
#Model error:
sigma = pm.Exponential(“sigma”, 1.0)
#Likelihood
y = pm.Normal(“y”, mu, sigma=sigma, observed=log_electricity, dims=“obs_id”)

with partial_pooling:
partial_pooling_trace = pm.sample(random_seed=RANDOM_SEED)
partial_pooling_idata = az.from_pymc3(partial_pooling_trace)

My problem is that now I can’t seem to find the way to get from the estimated parameter distributions to the predictions. For example my goal would be to estimate what would the electricity consumption be at a certain hour of the day and given a certain outdoor temperature value.

All the notebook examples I found define a model, then obtain the trace by fitting the model using sampling or some other methods. None of the notebooks explains how to use that trace to obtain predictions. My idea is that this should happen through some kind of sampling from the obtained trace, but I’m not sure if I’m missing something important here.

Thanks a lot in advance for the help!

1 Like

Yes, you do need to sample: from the posterior predictive using pm.sample_posterior_predictive(). Note that you also need to set the value you want to predict for. You can do this using pm.Data container so that during prediction you can change those values. See here for more information: Prior and Posterior Predictive Checks — PyMC3 3.10.0 documentation

2 Likes

Thank you very much for the answer Thomas. I was indeed able to calculate my predictions by using the pm.sample_posterior_predictive() function as described in the notebook you linked.

I’m now trying to calculate HDIs for my predictions using the az.hdi() function, but I’ve been encountering different issues. I would like to be able to estimate an HDI interval for my predictions, and my idea is that in order to obtain that I should calculate HDIs for each of the posterior distributions of the parameters in my model, and then combine them as specified by the model. But I’m not sure if this is the correct way, and I’m also not sure of how to obtain said HDIs, since when I try to run for example

az.hdi(posterior_predictive['a_daypart'])

I get the error:

UserWarning: More chains (4000) than draws (6). Passed array should have shape (chains, draws, *shape)

which I think is related to the fact that my posterior predictive is really just 1 chain with 4000 draws, while 6 is the amount of different ‘dayparts’, but arviz doesn’t get that.

Overall I’m just really confused about how the arviz data structure works and how to calculate HDIs in general and I’d be really grateful if someone could shed some light or point to any resources that can help me!

It looks like ArviZ is just unhappy with the data you put into it. You can reshape the array contained in posterior_predictive['a_daypart'] to fit the indicated shape it expects. I’m guessing that currently, posterior_predictive['a_daypart'].shape is something like (4000, 6), so I think you can reshape it to (1, 4000, 6) by just doing posterior_predictive['a_daypart'][None, :, :]. I hope this helps!

Somehow I missed the initial question. ArviZ expects the array data to either have labeled dimensions with two of them being chain and draw or to follow the convention chain, draw, *shape if dimensions are not labeled.

To get the result of sample_posterior_predictive to generate arrays following ArviZ convention, you can use the keep_size argument.

There are also some examples of using hdi with multilevel models, even calculating hdi on groupby objects in the “radon” notebook: A Primer on Bayesian Methods for Multilevel Modeling — PyMC3 3.10.0 documentation (hdi is used several times so I am not linking to any specific section, search hdi within the page to see all the occurences)