# PyMC3 posterior prediction example; how to incorporate a matrix of betas?
**URL:** https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251
**Category:** Questions
**Created:** [April 15, 2021, 7:59pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251 "2021-04-15T19:59:29Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 15, 2021, 7:59pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/1 "2021-04-15T19:59:29Z")
While working through some examples, I am stuck on how to incorporate a matrix of covariates in this example: [Prior and Posterior Predictive Checks — PyMC3 3.10.0 documentation](https://docs.pymc.io/notebooks/posterior_predictive.html)
`mu_pp = (ppc["a"] + ppc["b"] * predictor_scaled[:, None]).T`
Let’s say that I also generated `ppc['covariates']` which is a shape `[samples, # of covariates (56 in my case)]` and I want to incorporate that into the equation above. This runs into broadcasting issues which I am not familiar with solving:
`operands could not be broadcast together with shapes (2000,) (2000,56)`
Does every object then need to be broadcasted against this shape, `[2000, 56]`?
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 15, 2021, 9:46pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/2 "2021-04-15T21:46:56Z")
Hi,
I would recommend using `InferenceData` so that xarray takes care of all broadcasting automatically. There are some examples of this in the [rugby analytics example](https://nbviewer.jupyter.org/github/pymc-devs/pymc-examples/blob/main/examples/case_studies/rugby_analytics.ipynb) and in the [multilevel model example](https://docs.pymc.io/notebooks/multilevel_modeling.html). This second example has examples of xarray usage but with the posterior variables, not with posterior predictive. Once you have generated the inferencedata, there are no practical differences in operating with posterior, prior or posterior predictive.
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 15, 2021, 11:19pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/3 "2021-04-15T23:19:52Z")
Thanks so much, I will dig into these and get more accustom to the pm.Data containers!
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 16, 2021, 1:14am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/4 "2021-04-16T01:14:22Z")
I am actually still wracking my brain over this; if I made a MWE would you consider showing how to plot the predicted relationship between the predictor and the outcome when using `idata` and many control parameters in addition to the treatment?
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 16, 2021, 2:57pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/5 "2021-04-16T14:57:28Z")
Attached is a synthetic example showing where I am having trouble! (Table for the moment that a linear model is not the best fit for Likert-scale data).
```auto
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import arviz as az
import theano
# generate Likert y; 1-7
obs_y = np.random.randint(7, size=(100, 1)) + 1
# generate treatment and control group
treatment = np.random.randint(3, size=(100,))
# convert to one-hot;
treatment = pd.get_dummies(data=treatment, drop_first=True)
# generate many categorical controls
covariates = np.random.randint(2, size=(100, 56))
# coords for pymc3
coords = {
'obs_id': np.arange(len(obs_y))
}
# specify model
with pm.Model(coords=coords) as main_lm_model:
'''
linear model replication
'''
# priors
alpha = pm.Normal('alpha', mu=0, sigma=1)
treatment_betas = pm.Normal("treatment_betas", mu=0, sigma=1, shape=treatment.shape[1])
covariates_betas = pm.Normal("covariates_betas", mu=0, sigma=1, shape=covariates.shape[1])
# model error
sigma = pm.HalfNormal("sigma", sigma=1)
# matrix-dot products
m1 = pm.math.matrix_dot(treatment, treatment_betas)
# i separate these here so its easier to just look at these values only
# through the var_names argument for different checks so I dont have to look
# at 56 covariates and draw their chains, etc
m2 = pm.math.matrix_dot(covariates, covariates_betas)
# expected value of y
mu = alpha + m1 + m2
# Likelihood: Normal
y = pm.Normal("y", mu=mu, sigma=sigma, observed=obs_y, dims='obs_id')
# set step and sampler
step = pm.NUTS([alpha, treatment_betas, covariates_betas, sigma], target_accept=0.9)
# prior analysis
with main_lm_model:
prior_pc = pm.sample_prior_predictive(500)
# Inference button (TM)!
with main_lm_model:
lm_trace = pm.sample(draws=1000,
step=step,
init='jitter+adapt_diag',
cores=4,
tune=500, # burn in
return_inferencedata=False) # seemingly have to do this to get idata?
# ppc
with main_lm_model:
ppc = pm.fast_sample_posterior_predictive(trace=lm_trace,
random_seed=1,
# have to specify the params here unlike prior
# or will only return 'y'
var_names=['alpha', 'treatment_betas', 'covariates_betas', 'y'])
# build idata - this will fail if return_inferencedata=True
lm_idata = az.from_pymc3(trace=lm_trace,
model=main_lm_model,
prior=prior_pc,
posterior_predictive=ppc)
# stack posterior samples
stacked = lm_idata.posterior.stack(draws=("chain", "draw"))
# try to calculate mu across the posterior like:
# https://docs.pymc.io/notebooks/posterior_predictive.html
mu_pp = ( stacked['alpha'] + stacked['covariates_betas'] + (stacked['treatment_betas'] * treatment) )
# yields dataframe coercion errors
```
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 16, 2021, 3:59pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/6 "2021-04-16T15:59:24Z")
Quick answer without playing with the code yet. treatment is still a dataframe/series, not an xarray object, so this should not work. Two main approches come to mind:
Try doing `treatment.to_xarray()`, pandas and xarray allow to automatically convert data from one to the other.
```auto
# specify model
with pm.Model(coords=coords) as main_lm_model:
t = pm.Data("treatment", treatment, dims=("obs_id", "treatment_covariate")
...
# then get treatment as
lm_idata.constant_data["treatment"]
```
Note that you are currently using `*` when I think you want a matrix product, not a pointwise product. You should be able to use `np.linalg.matmul` directly, as it acts on the last dimensions only and batches over any initial dimension as necessary (chain and draw in this case).
> [@localh](#):
>
> `# build idata - this will fail if return_inferencedata=True`
Yes, if you get an inferencedata from sample you need to use `extend` like it’s done in the rugby example.
side comment, you already have `alpha`, `treatment_betas`… in the posterior, “sampling” them in the posterior predictive is actually not doing anything.
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 16, 2021, 4:50pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/7 "2021-04-16T16:50:14Z")
Thanks for your time and help! It is much appreciated in helping me get the right structure so that it may be used for all other problems I have.
I am running into issues using pm.Data as it is not liking my dimensions, but I am not sure why:
```auto
# generate Likert y; 1-7
obs_y = np.random.randint(7, size=(100, 1)) + 1
# generate treatment and control group
treatment = np.random.randint(3, size=(100,))
# convert to one-hot;
treatment = pd.get_dummies(data=treatment, drop_first=True)
# generate many categorical controls
covariates = np.random.randint(2, size=(100, 56))
# treatment to xarray
treatment = treatment.to_xarray().to_array() # tried both; to_xarray()
treatment.shape # (2, 100)
# coords for pymc3
coords = {
'obs_id': np.arange(len(obs_y)),
'treatment_betas': treatment # (2, 100)
}
# specify model
with pm.Model(coords=coords) as main_lm_model:
t = pm.Data('treatment_betas', treatment, dims=('treatment_betas'))
# Length of `dims` must match the dimensions of the dataset. (actual 1 != expected 2)
```
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 16, 2021, 6:47pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/8 "2021-04-16T18:47:23Z")
It’s important to distinguish between model parameters, constant data and dimensions/coordinates. This post should help [PyMC3 with labeled coords and dims | Oriol unraveled](https://oriolabril.github.io/oriol_unraveled/python/arviz/pymc3/xarray/2020/09/22/pymc3-arviz.html). I’ll try to take the example amd write an executable version with it at some point of the weekend, this is the best I can offer for now
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 16, 2021, 6:53pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/9 "2021-04-16T18:53:49Z")
Thanks so much; I have been trying a bunch of different solutions but will try this as well!
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 16, 2021, 8:37pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/10 "2021-04-16T20:37:31Z")
Here is my (_working_) stab at working with the examples (again tabling that a LM is not “right” here).
I would be super interested in your feedback on this general framework setup for acquiring the different objects!
```auto
# generate Likert y; 1-7
obs_y = np.random.randint(7, size=(100,)) + 1
obs_y = pd.DataFrame(obs_y)
# generate treatment and control group
treatment = np.random.randint(3, size=(100,))
# convert to one-hot;
treatment = pd.get_dummies(data=treatment, drop_first=True)
treatment.columns = ['treat_{}'.format(x) for x in range(2)]
# generate many categorical controls
covariates = np.random.randint(2, size=(100, 56))
covariates = pd.DataFrame(covariates) # to df
covariates.columns = ['var_{}'.format(x) for x in range(56)]
# coords for pymc3
coords = {
'obs_id': np.arange(len(obs_y)),
'treatment': treatment.columns,
'covariates': covariates.columns
}
# specify model
with pm.Model(coords=coords) as main_lm_model:
# data containers
covariates_data = pm.Data('covariates_data', covariates, dims=('obs_id', 'covariates'))
treatment_data = pm.Data('treatment_data', treatment, dims=('obs_id', 'treatment'))
# priors
alpha = pm.Normal('alpha', mu=0, sd=1)
covariates_betas = pm.Normal("covariates_betas", mu=0, sd=1, dims='covariates')
treatment_betas = pm.Normal("treatment_betas", mu=0, sd=1, dims='treatment')
# model error
sigma = pm.HalfNormal("sigma", sigma=1)
# matrix-dot products
m1 = pm.math.matrix_dot(treatment, treatment_betas)
m2 = pm.math.matrix_dot(covariates, covariates_betas)
# expected value of y
mu = alpha + m1 + m2
# Likelihood: Normal
y = pm.Normal("y",
mu=mu,
sigma=sigma,
observed=obs_y.values.reshape(-1,),
dims='obs_id')
# set sampler
step = pm.NUTS([alpha, treatment_betas, covariates_betas, sigma], target_accept=0.9)
# Inference button (TM)!
with main_lm_model:
lm_trace = pm.sample(draws=1000,
step=step,
init='jitter+adapt_diag',
cores=4,
tune=500, # burn in
return_inferencedata=True) # seemingly have to do this to get idata?
# prior analysis
with main_lm_model:
prior_pc = pm.sample_prior_predictive()
# posterior analysis
with main_lm_model:
# ppc
ppc = pm.fast_sample_posterior_predictive(trace=lm_trace,
random_seed=1,
# actually have to specify these
# or you cant get them for calcs
var_names=['y', 'alpha', 'treatment_betas', 'covariates_betas', 'sigma']
)
# turn to idata
ppc = az.from_pymc3(posterior_predictive=ppc)
# seems like it already stacks chains + draws?
ppc.posterior_predictive['covariates_betas'].shape # (1, 4000, 56)
# plot ppc
az.plot_ppc(data=ppc, num_pp_samples=100);
# try to calculate mu across the posterior like:
# https://docs.pymc.io/notebooks/posterior_predictive.html
treatment_mm = np.matmul(
# remove chain dim
np.squeeze(ppc.posterior_predictive['treatment_betas'].values),
lm_trace.constant_data['treatment_data'].T.values
)
covariate_mm = np.matmul(
# remove chain dim
np.squeeze(ppc.posterior_predictive['covariates_betas'].values),
lm_trace.constant_data['covariates_data'].T.values
)
# broadcast alpha
mu_pp = np.squeeze(ppc.posterior_predictive['alpha'].values)[:, None] + treatment_mm + covariate_mm
mu_pp.shape # 4000 samples, 100 obs
mu_pp_mean = mu_pp.mean(0) # get mean down dim 0
_, ax = plt.subplots()
ax.plot(treatment.iloc[:, 0], obs_y.values.reshape(-1), "o", ms=4, alpha=0.4, label="Data")
ax.plot(treatment.iloc[:, 0], mu_pp_mean, label="Mean outcome", alpha=0.6)
az.plot_hpd(
treatment.iloc[:, 0],
mu_pp,
ax=ax,
fill_kwargs={"alpha": 0.8, "label": "Mean outcome 94% HPD"},
)
az.plot_hpd(
treatment.iloc[:, 0],
ppc.posterior_predictive['y'],
ax=ax,
fill_kwargs={"alpha": 0.8, "color": "#a1dab4", "label": "Outcome 94% HPD"},
)
```
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 17, 2021, 1:32am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/11 "2021-04-17T01:32:31Z")
> [@localh](#):
>
> ```auto
> with main_lm_model:
> lm_trace = pm.sample(draws=1000,
> step=step,
> init='jitter+adapt_diag',
> cores=4,
> tune=500, # burn in
> return_inferencedata=True) # seemingly have to do this to get idata?
>
> ```
If you plan to sample both posterior and posterior predictive virtually at the same time, you can use `False` and call `az.from_pymc3` with all the arguments.
> [@localh](#):
>
> ```auto
> with main_lm_model:
> # ppc
> ppc = pm.fast_sample_posterior_predictive(trace=lm_trace,
> random_seed=1,
> # actually have to specify these
> # or you cant get them for calcs
> var_names=['y', 'alpha', 'treatment_betas', 'covariates_betas', 'sigma']
> )
> # turn to idata
> ppc = az.from_pymc3(posterior_predictive=ppc)
>
> ```
`alpha`, `treatment_betas`… are variables of the posterior. You already have their posterior samples in `lm_trace.posterior["alpha"]`. `pm.sample_posterior_predictive` won’t sample them again, it will use the values provided in `lm_trace`. You have to specify them in order to get them to appear here because the default is to sample only from the variables in the posterior predictive. The posterior and posterior predictive distributions are different, the first is a distribution over the parameter space, the second over the observations space.
As a side note, you can get the values with the right shape with
```auto
ppc = pm.sample_post...(..., keep_size=True)
az.from_dict(posterior_predictive=ppc, dims={"y": ["obs_id"]}, coords=coords)
```
> [@localh](#):
>
> ```auto
> treatment_mm = np.matmul(
> # remove chain dim
> np.squeeze(ppc.posterior_predictive['treatment_betas'].values),
> lm_trace.constant_data['treatment_data'].T.values
> )
>
> ```
From the numpy docs, it looks like neither `np.dot` nor `matmul` are the right choice for this situation, so I’d go with einsum. Moreover, I have seen that linalg operations on dataarrays work but loose are converted to arrays in the process, so to maintain that information we need to use `xr.apply_ufunc`.
```python
import xarray as xr
post = lm_trace.posterior
const = lm_trace.constant_data
post["treatment_mm"] = xr.apply_ufunc(
np.einsum, "...i,...i", post["treatment_betas"], const["treatment_data"], input_core_dims=[[], ["treatment"], ["treatment"]])
```
Using `input_core_dims` will ensure that `treatment` is the last dimension of the array, which is the one we are dotting over with einsum thanks to the ellipsis. The same thing would have to be done for the covariates. This is a bit verbose, but the good thing is that now you can do:
```auto
mu = post["alpha"] + post["treatment_mm"] + post["covariate_mm"]
```
And everything will be broadcast together automatically. Another very cool thing, even it’s probably not useful in most cases, is that all the xarray computation is now independent of the order and the number of dimensions, everything works automagically. You could add an extra `replication` dimension if you perform the same treatment to the same subjects multiple times if that makes any sense, and the code would still work.
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 17, 2021, 3:24pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/12 "2021-04-17T15:24:50Z")
This is incredible help; I just went through it all and played with all on my end. It is making so much more sense! I am excited to port it over to a “real” linear model setup with a continuous DV and then extend from there.
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 18, 2021, 9:49pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/13 "2021-04-18T21:49:26Z")
Last few questions on this topic if I may (here with actual regression data)!
I am curious about a few things:
1. Is stacking our InferenceData items recommended and if so, is this a recommended way? `.posterior.stack(draws=("chain", "draw"))` e.g., `chain, draw` or `draw, chain`?
2. I am curious why my expected value of y is not plotting a straight line akin to official example. The results I have also seem very counterintuitive to a VERY certain HDI result (see image below). ([http://docs.pymc.io/notebooks/posterior\_predictive.html](http://docs.pymc.io/notebooks/posterior_predictive.html)).
Thanks again for all your time. I have spent most the weekend playing with pymc3; its super addicting and very useful.
MWE:
```auto
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import arviz as az
import theano
import xarray as xr
from sklearn.datasets import make_regression
X, target, coef = make_regression(n_samples=100,
bias=1.0,
n_informative=3,
n_features=3,
noise=2.5,
random_state=1,
coef=True)
target = target
treatment = X[:, 0]
covariates = X[:, 1:]
# coords
coords = {
# dim 1: length of data set
'obs_id': np.arange(len(target)),
# dim 2: treatment var
'treatment': ['treatment'],
# dim 3: covariates
'covariates': ['covariate_1', 'covariate_2']
}
# specify model
with pm.Model(coords=coords) as sk_lm:
# data
treatment_data = pm.Data('treatment_data', treatment, dims=('obs_id'))
covariates_data = pm.Data('covariates_data', covariates, dims=('obs_id', 'covariates'))
# priors
# intentionally off the mark
alpha = pm.Normal('alpha', mu=0, sigma=1)
treatment_beta = pm.Normal("treatment_beta", mu=30, sigma=2)
covariates_betas = pm.Normal('covariates_betas', mu=[66, 33], sigma=2, dims=('covariates'))
# model error
sigma = pm.HalfNormal("sigma", sigma=1)
# matrix-dot products
m1 = pm.math.matrix_dot(covariates, covariates_betas)
# expected value of y
mu = alpha + (treatment_beta * treatment) + m1
# Likelihood: Normal
y = pm.Normal("y",
mu=mu,
sigma=sigma,
observed=target,
dims='obs_id')
# set sampler
step = pm.NUTS([alpha, treatment_beta, covariates_betas, sigma], target_accept=0.9)
# prior analysis
with sk_lm:
prior_pc = pm.sample_prior_predictive()
# Inference button (TM)!
with sk_lm:
lm_trace = pm.sample(draws=1500,
step=step,
init='jitter+adapt_diag',
cores=4,
tune=500, # burn in
return_inferencedata=False)
# posterior analysis
with sk_lm:
ppc = pm.fast_sample_posterior_predictive(trace=lm_trace,
random_seed=1,
)
# generate InferenceData
with sk_lm:
lm_idata = az.from_pymc3(
trace=lm_trace,
prior=prior_pc,
posterior_predictive=ppc,
coords=coords
)
# Posterior Analysis
# extract posterior data for all parameters; stack chains + draws
post = lm_idata.posterior.stack(draws=("chain", "draw"))
# extract the data used to build the model
const = lm_idata.constant_data
# check data
post['treatment_beta'].shape # (6000, )
# matrix multiply covariate matrix by posterior covariates betas matrix
post["covariate_mm"] = xr.apply_ufunc(np.einsum, "...i,...i",
post["covariates_betas"],
const["covariates_data"],
input_core_dims=[[],
["covariates"],
["covariates"]])
```
## Area of Concern
```auto
# roughly following http://docs.pymc.io/notebooks/posterior_predictive.html
# generate expected value of y; alpha + beta1 + covariates
mu = post["alpha"] + (post['treatment_beta'] * const['treatment_data']) + post["covariate_mm"]
# mean outcome plot
_, ax = plt.subplots(figsize=(8,8))
# why is this not a straight line? also seems counterintuitive to the HDI results
ax.plot(np.sort(treatment), np.sort(mu.mean('draws')), label="Mean outcome", alpha=0.6)
# plot data
ax.plot(treatment, target, "o", ms=4, alpha=0.4, label="Data")
# plot 94% HPD for E(y)
# arviz.plot_hdi has superceded arvis.plot_hpd?
az.plot_hdi(x=lm_idata.constant_data.treatment_data, # observed treatment
y=mu, # expected value of y
ax=ax,
smooth=False,
hdi_prob=0.94,
fill_kwargs={"alpha": 0.8, "label": "Mean outcome 94% HPD"})
# plot 94% HPD for posterior possible unobserved y values
az.plot_hdi(x=lm_idata.constant_data.treatment_data, # observed treatment data
y=stacked_posterior_predictive['y'].T, # set of possible unobserved y
ax=ax,
smooth=False,
hdi_prob=0.94,
fill_kwargs={"alpha": 0.8, "color": "#a1dab4", "label": "Outcome 94% HPD"})
ax.set_xlabel("Predictor")
ax.set_ylabel("Outcome")
ax.set_title("Posterior predictive checks")
ax.legend(ncol=2, fontsize=10);
```
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 18, 2021, 10:37pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/14 "2021-04-18T22:37:24Z")
> [@localh](#):
>
> Is stacking our InferenceData items recommended and if so, is this a recommended way? `.posterior.stack(draws=("chain", "draw"))` e.g., `chain, draw` or `draw, chain` ?
It should not really matter. Some plots and other functions will need to have the original chain, draw shape, but you can always unstack. I personally prefer not stacking unless necessary, but as we say in catalan “cada mestre té el seu llibre” (literal translation “every teacher has their own book”)
> [@localh](#):
>
> I am curious why my expected value of y is not plotting a straight line akin to official example. The results I have also seem very counterintuitive to a VERY certain HDI result (see image below).
There are several issues in the sample code.
> [@localh](#):
>
> ```auto
> # why is this not a straight line? also seems counterintuitive to the HDI results
> ax.plot(np.sort(treatment), np.sort(mu.mean('draws')), label="Mean outcome", alpha=0.6)
>
> ```
This is sorting `treatment` and `mu` independently! You therefore have no guarantees that the 1-1 relation between treatment and their corresponding mu is maintained. You can use `np.argsort` and handle broadcasting, `:`… but again xarray can help you with that. Take a look for example at code cell 21 of [https://docs.pymc.io/notebooks/multilevel\_modeling.html](https://docs.pymc.io/notebooks/multilevel_modeling.html).
> [@localh](#):
>
> ```auto
> # arviz.plot_hdi has superceded arvis.plot_hpd?
> az.plot_hdi(x=lm_idata.constant_data.treatment_data, # observed treatment
> y=mu, # expected value of y
>
> ```
This is not sorted at all, unlike the plot above, so you should **not** expect those to be comparable at all, ever. You are also plotting the HDI for `mu`, not for `y`, so you are plotting the hdi of the _mean_ of your predictions, not of the predictions themselves. These are very different things. See for example code cell 16 in [https://docs.pymc.io/notebooks/multilevel\_modeling.html](https://docs.pymc.io/notebooks/multilevel_modeling.html).
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 19, 2021, 3:03am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/15 "2021-04-19T03:03:59Z")
Thanks so much! I have been digging into the linked notebook and stumbled across another problem. It requires the treatment variable to be inside the `InferenceData` object.
`level_labels = pooled_idata.posterior.Level[pooled_idata.constant_data.floor_idx]`
I really would like to keep my treatment beta separate from the covariates and find it convenient to do so; especially in higher dimensional orders. The problem is, I cannot get my treatment coords inside the `InferenceData` object when it is a single beta parameter (coords appears to be only for multi-dimensional items).
```auto
# coords
coords = {
# dim 1: length of data set
'obs_id': np.arange(len(target)),
# dim 2: treatment var
'treatment': ['treatment'],
# dim 3: covariates
'covariates': ['covariate_1', 'covariate_2']
}
# specify model
with pm.Model(coords=coords) as sk_lm:
# data
treatment_data = pm.Data('treatment_data', treatment, dims=('obs_id'))
covariates_data = pm.Data('covariates_data', covariates, dims=('obs_id', 'covariates'))
# priors
alpha = pm.Normal('alpha', mu=0, sigma=1)
treatment_beta = pm.Normal("treatment_beta", mu=30, sigma=2, dims=('treatment')) # fails
covariates_betas = pm.Normal('covariates_betas', mu=[66, 33], sigma=2, dims=('covariates'))
```
If I include a dims argument (which would solve the issue) and link it to the coords, then:
```auto
# prior analysis
with sk_lm:
prior_pc = pm.sample_prior_predictive()
```
Throws a value error, which after searching, seems to be associated with the shape of `mu` and it does not like specifying a shape value size of 1.
`ValueError: 0-dimensional argument does not have enough dimensions for all core dimensions ('i_1_0',)`
How might I account for a single parameter such that it can be linked to the coords and therefore available under `InferenceData.posterior`?
I have added a category manually like so in the meantime:
`lm_idata.posterior.coords['treatment'] = ['treatment']`
# Edit
When I get home I hope to try using an idata kwarg, drawing on: [Redirecting to Oriol Unraveled](https://oriolabril.github.io/oriol_unraveled/python/arviz/pymc3/xarray/2020/09/22/pymc3-arviz.html)
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 20, 2021, 1:36am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/16 "2021-04-20T01:36:17Z")
> [@localh](#):
>
> The problem is, I cannot get my treatment coords inside the `InferenceData` object when it is a single beta parameter (coords appears to be only for multi-dimensional items).
It looks like you are using an array with one dimension of length `obs_id` so you effectively can’t use `obs_id, treatment` as dims, your array only has one dimension. You should be using an array with two dimensions `obs_id, treatment` and have the treatment dimension have length 1. An easy way to do that from a 1d treatment array is `treatment[:, None]`
As for the error, I don’t really understand it, `treatment_beta` looks fine to me and I don’t know what `i_1_0` is.
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 20, 2021, 3:11am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/17 "2021-04-20T03:11:55Z")
Thanks again for your patience and help!
I think I have come pretty far and I would be curious to see your thoughts; I think the last issue remaining is that I cannot quite figure out your last recommendation about the dimensions. Could you offer a quick code example? I tried playing with substituting `treatment[:, None]` in a few areas to no avail.
## Setup
```auto
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import arviz as az
import theano
import xarray as xr
from sklearn.datasets import make_regression
X, target, coef = make_regression(n_samples=100,
bias=1.0,
n_informative=3,
n_features=3,
noise=2.5,
random_state=1,
coef=True)
target = target
treatment = X[:, 0]
covariates = X[:, 1:]
# coords
coords = {
# dim 1: treatment var
'treatment': ['treatment_1'],
# dim 2: covariates
'covariates': ['covariate_1', 'covariate_2'],
# dim 3: len of df
'obs_id': np.arange(len(target))
}
# specify model
with pm.Model(coords=coords) as sk_lm:
# data
treatment_data = pm.Data('treatment_data', treatment, dims=('obs_id'))
covariates_data = pm.Data('covariates_data', covariates, dims=('obs_id', 'covariates'))
# priors
alpha = pm.Normal('alpha', mu=0, sigma=1)
treatment_beta = pm.Normal("treatment_beta", mu=44, sigma=2) # really does not want dims=('treatment')
covariates_betas = pm.Normal('covariates_betas', mu=[87, 58], sigma=2, dims=('covariates'))
# model error
sigma = pm.HalfNormal("sigma", sigma=1)
# matrix-dot products
m1 = pm.math.matrix_dot(covariates, covariates_betas)
# expected value of y
mu = alpha + (treatment_beta * treatment) + m1
# Likelihood: Normal
y = pm.Normal("y",
mu=mu,
sigma=sigma,
observed=target,
dims='obs_id')
# set sampler
step = pm.NUTS([alpha, treatment_beta, covariates_betas, sigma], target_accept=0.9)
# Inference button (TM)!
lm_trace = pm.sample(draws=1000,
step=step,
init='jitter+adapt_diag',
cores=4,
tune=500, # burn in
return_inferencedata=False)
# prior analysis
prior_pc = pm.sample_prior_predictive()
# posterior predictive
ppc = pm.fast_sample_posterior_predictive(trace=lm_trace,
random_seed=1,
)
# inference data
lm_idata = az.from_pymc3(
trace=lm_trace,
prior=prior_pc,
posterior_predictive=ppc,
)
```
## Posterior Analysis
```auto
# Posterior Analysis
post = lm_idata.posterior
# extract the data used to build the model
const = lm_idata.constant_data
# generate expected value of y given treatment; alpha + beta1
mu = (post["alpha"] + (post['treatment_beta'] * const['treatment_data'])) # (4, 1000, 100)
# argsort init
argsort1 = np.argsort(treatment)
# manually generate HDI information
hdi_data_yhat = az.hdi(lm_idata.posterior_predictive['y'])
hdi_data_eY = az.hdi(mu)
# Mean Outcome Plot
# start plot
_, ax = plt.subplots(figsize=(8,8))
# plot data
ax.plot(treatment, target, "o", ms=4, alpha=0.4, label="Data")
# mean outcome; # E(y) for one unit change in alpha + b1*x1
ax.plot(treatment, mu.mean(dim=("chain", "draw")), label="Mean outcome", alpha=0.6)
# plot 94% HPD for predicted y-hats; given trained model
# arviz.plot_hdi has superceded arvis.plot_hpd?
az.plot_hdi(x=treatment, # observed treatment
hdi_data=hdi_data_eY, # expected 94% HPD for y-hat
ax=ax,
hdi_prob=0.94,
fill_kwargs={"alpha": 0.8, "label": "Mean outcome 94% HPD"})
# plot 94% HPD for posterior possible unobserved y values
az.plot_hdi(x=treatment, # observed treatment data
hdi_data=hdi_data_yhat, # set of possible unobserved y
ax=ax,
hdi_prob=0.94,
fill_kwargs={"alpha": 1.0, "color": "#a1dab4", "label": "Outcome 94% HPD"})
ax.set_xlabel("Predictor")
ax.set_ylabel("Outcome")
ax.set_title("Posterior predictive checks")
ax.legend(ncol=2, fontsize=10);
```
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 20, 2021, 10:11pm UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/18 "2021-04-20T22:11:19Z")
Should be something like
```auto
treatment_data = pm.Data('treatment_data', treatment[:, None], dims=('obs_id', "treatment"))
```
---
### Author:  [@localh](https://discourse.pymc.io/u/localh)
#### Post date: [April 21, 2021, 2:46am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/19 "2021-04-21T02:46:18Z")
And that, of course worked! Thank you so much for all the help – I feel like I owe you money; its been great to learn!
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [April 21, 2021, 10:31am UTC](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251/20 "2021-04-21T10:31:40Z")
You can donate some money via NumFOCUS if you found PyMC3, ArviZ and my help useful, it will help us continue improving the libraries even more
- Donate to ArviZ: [Donate for Arviz](https://numfocus.org/donate-to-arviz)
- Donate to PyMC3: [Donate to PyMC3](https://numfocus.org/donate-to-pymc3)
[Next page](https://discourse.pymc.io/t/pymc3-posterior-prediction-example-how-to-incorporate-a-matrix-of-betas/7251.md?page=2)