Multilevel model with covariance of 3 variables

Questions
1

Hi all, I’m trying to build a multilevel model that should estimate the electricity consumption of a building based on:

  • an intercept that depends on the day of the week,
  • a temperature slope coefficient for days when cooling is needed,
  • a temperature slope coefficient for days when heating is needed.

I first estimated a model with varying intercepts and slopes, and it worked fine. I would now like to explore the possibility of adding the covariance between these 3 variables as a predictor. Unfortunately I only managed to find examples and notebooks to analyse the covariance between 2 variables in multilevel models.

Is it possible to estimate the covariance between 3 variables in pymc3? If not, would it be possible (and would it make sense) to add instead the individual covariance between the intercept and each of the two slopes?

Thanks in advance for your help

1 Like
2

Hi there,

I’m not sure I completely understand what you mean, but if I understand it right you would like to use multiple correlated random effects. The way I would probably approach that is by putting an LKJ prior on the covariance of the random effects. This example has only two correlated coefficients, but if instead of:

    chol, corr, stds = pm.LKJCholeskyCov("chol", n=2, eta=2.0, sd_dist=sd_dist, compute_corr=True)

you write

    chol, corr, stds = pm.LKJCholeskyCov("chol", n=3, eta=2.0, sd_dist=sd_dist, compute_corr=True)

you will get a LKJ prior on a covariance matrix for three random variables rather than 2. You can then use that Cholesky factor to put a prior on your coefficients, e.g.:

vec = pm.MvNormal("coefs", mu=np.zeros(3), chol=chol, shape=(n_buildings, 3))

this should give you n_buildings vectors, each of length 3, and estimating the covariance between these. Of course the mean doesn’t have to be zero, you could do something else.

Let me know if that makes sense or this is not what you were looking for!

1 Like
3

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!

4

Hi again,

sorry for the delayed response! Here’s what I would try, if I understand you correctly:

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)
   bh = pm.Normal("bh", mu=0.0, sigma=1.0)
   bc = pm.Normal("bc", mu=0.0, sigma=1.0)

   sd_dist = pm.Exponential.dist(0.5)

   chol, corr, stds = pm.LKJCholeskyCov("chol", n=3, eta=2.0, sd_dist=sd_dist, compute_corr=True)

   # Correlated varying intercept and slopes within clusters
   # Note that I don't really know how to use the dims argument, but you could replace
   # the shape bit with something using that I'm sure.
   coefs = pm.MvNormal("coefficients", mu=tt.stack([a, bh, bc]), chol=chol, shape=(n_clusters, 3))

   # You can now pick out the a, bh and bc if you like and they should be correlated within clusters
   a_cluster = coefs[:, 0]
   bh_cluster = coefs[:, 1]
   bc_cluster = coefs[:, 2]

   # 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")

Since I don’t have your data, I haven’t tested this, but I hope this makes sense and I’m not misinterpreting what you’re looking for.

1 Like
5

Thanks a lot for your help, Martin. The model you posted computes and provides pretty satisfactory results too.

I have two more questions, related to this model, that maybe you might help me with.

First, I was trying to write this model in mathematical form. Following this section of the radon notebook:

image
image2282×482 47.6 KB

I came up with the following:

image
image1338×832 34.8 KB

Do you think this mathematical representation is correct for the model you sent? What would the matrix P look like?

Also I’m trying to a obtain graphviz representation of this model, but for some reason when I try to use graphviz with covarying models, through the command

pm.model_to_graphviz(covarying_intercept_and_temp)

I get the following error:

AttributeError: 'float' object has no attribute 'name'

It’s a pretty vague error and I can’t seem to debug it. Graphviz works fine with every other model, but as soon as I introduce covarying terms this happens. Might be a long shot but maybe someone here encountered a similar issue at some point.

Thanks a lot for your help again!

1 Like
6

Hi Benedetto,

Your mathematical form looks good to me!

The matrix P is a correlation matrix. I think pymc3 should return you samples from chol_corr, which are the entries of this matrix. These entries will tell you how correlated the coefficients are. For example, the entry chol_corr[0, 1] will give you samples from the correlation between \alpha and \beta_h within clusters. If it’s high, then a cluster with a high \alpha will tend to also have a high \beta_h, for example. Similarly, chol_corr[0, 2] will give you samples from the correlation between \alpha and \beta_c. Let me know if that makes sense.

I’m afraid I haven’t used the graphviz functionality in pymc3, so I don’t know what’s causing that error – hopefully someone else knows!

Best,
Martin

1 Like
7

Yes. The same issue was reported in #4391 and fixed by @Spaak . PyMC3 3.11.0 release will include this fix :slight_smile:

2 Likes
8

Many thanks to both! Great to find support here.

9

Hello grillone and Martin,
I’m a beginner with pymc and am looking for some input.
My goal is to create probabilistic load forecasting for a low-voltage network.
I start with simple load profiles in csv.
I would then like to add further factors, for example PV grid feed-in.
Would it be possible to take a look at your entire code?
can? This would help me understand the structure.

Best regards
Tobi

10

Hi Tobi, I was researching on this as part of my PhD, we managed to publish a scientific paper out of it but I didn’t continue on this line of research. This said, I have a couple of links that you might find useful:

  1. The scientific paper: Energies | Free Full-Text | Baseline Energy Use Modeling and Characterization in Tertiary Buildings Using an Interpretable Bayesian Linear Regression Methodology

  2. The Github repo with the code I used for the paper: GitHub - bgrillone/Bayesian-Measurement-and-Verification

The Github repo is not well documented and was meant for my personal use only, but maybe you can have a look and figure out if it can help you with your project as well.

Good luck!