Bayesian Workflow with SEMs

[NathanielF]

Bayesian Workflow with SEMs

Related to proposal here
#806

• Notebook follows style guide https://docs.pymc.io/en/latest/contributing/jupyter_style.html
• PR description contains a link to the relevant issue:
  • a tracker one for existing notebooks (tracker issues have the "tracker id" label)
  • or a proposal one for new notebooks
• Check the notebook is not excluded from any pre-commit check: https://github.com/pymc-devs/pymc-examples/blob/main/.pre-commit-config.yaml

Helpful links

• https://github.com/pymc-devs/pymc-examples/blob/main/CONTRIBUTING.md

---

📚 Documentation preview 📚: https://pymc-examples--807.org.readthedocs.build/en/807/

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[NathanielF]

Apparent indexing issue between pymc versions 5.17 --> 5.30

Indexing trick works for pymc 5.17, but breaks on 5.30

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[15], [line 61](vscode-notebook-cell:?execution_count=15&line=61)
     [57](vscode-notebook-cell:?execution_count=15&line=57) priors = {"lambdas": [1, 0.5], "eta": 2, "B": [0, 0.5], "tau": [0, 1]}
     [59](vscode-notebook-cell:?execution_count=15&line=59) priors_wide = {"lambdas": [1, 5], "eta": 2, "B": [0, 5], "tau": [0, 10]}
---> [61](vscode-notebook-cell:?execution_count=15&line=61) sem_model_hierarchical_tight = make_hierarchical(priors, grp_idx)
     [62](vscode-notebook-cell:?execution_count=15&line=62) sem_model_hierarchical_wide = make_hierarchical(priors_wide, grp_idx)
     [64](vscode-notebook-cell:?execution_count=15&line=64) pm.model_to_graphviz(sem_model_hierarchical_tight)

Cell In[15], [line 52](vscode-notebook-cell:?execution_count=15&line=52)
     [50](vscode-notebook-cell:?execution_count=15&line=50)         Sigma_y.append(Sigma_y_g)
     [51](vscode-notebook-cell:?execution_count=15&line=51)     Sigma_y = pt.stack(Sigma_y)
---> [52](vscode-notebook-cell:?execution_count=15&line=52)     _ = pm.MvNormal("likelihood", mu=0, cov=Sigma_y[grp_idx], dims=('obs', 'indicators'))
     [54](vscode-notebook-cell:?execution_count=15&line=54) return sem_model_hierarchical

File ~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:310, in Distribution.__new__(cls, name, rng, dims, initval, observed, total_size, transform, *args, **kwargs)
    [307](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:307)     elif observed is not None:
    [308](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:308)         kwargs["shape"] = tuple(observed.shape)
--> [310](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:310) rv_out = cls.dist(*args, **kwargs)
    [312](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:312) rv_out = model.register_rv(
    [313](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:313)     rv_out,
    [314](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:314)     name,
   (...)
    [319](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:319)     initval=initval,
    [320](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:320) )
    [322](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/distribution.py:322) # add in pretty-printing support

File ~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:262, in MvNormal.dist(cls, mu, cov, tau, chol, lower, **kwargs)
    [259](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:259) @classmethod
    [260](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:260) def dist(cls, mu, cov=None, tau=None, chol=None, lower=True, **kwargs):
    [261](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:261)     mu = pt.as_tensor_variable(mu)
--> [262](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:262)     cov = quaddist_matrix(cov, chol, tau, lower)
    [263](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:263)     # PyTensor is stricter about the shape of mu, than PyMC used to be
    [264](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:264)     mu = pt.broadcast_arrays(mu, cov[..., -1])[0]

File ~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:123, in quaddist_matrix(cov, chol, tau, lower, *args, **kwargs)
    [121](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:121)     cov = pt.as_tensor_variable(cov)
    [122](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:122)     if cov.ndim != 2:
--> [123](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:123)         raise ValueError("cov must be two dimensional.")
    [124](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:124) elif tau is not None:
    [125](https://file+.vscode-resource.vscode-cdn.net/Users/nathanielforde/Documents/Github/pymc-examples/examples/case_studies/~/mambaforge/envs/pymc_examples_new/lib/python3.9/site-packages/pymc/distributions/multivariate.py:125)     tau = pt.as_tensor_variable(tau)

ValueError: cov must be two dimensional.

This works in pymc 5.17

alpha = np.random.normal(0, 1, size=(2, 4))
M = np.random.normal(0, 1, size=(2, 12, 4))
inv_I_minus_B = np.random.normal(0,1, size=(2, 4, 4))
Lambda = np.random.normal(0,1, size=(12, 4))

Sigma_y = np.random.normal(0, 1, size=(2, 12, 12))
print("Sigma_y shape", Sigma_y.shape)

print(np.matmul(Lambda, inv_I_minus_B).shape)

mu_y = np.matmul(alpha[:, None, :], M.transpose(0, 2, 1))[:, 0, :]
print("Mu_y shape", mu_y.shape)

[fonnesbeck]

I would ditch the Rhat plots -- all the action is in a tiny region just above 1.0, so most of the plot is irrelevant.

[NathanielF]

@ricardoV94 , not looking for a review, just wondering about the shape handling in the MvNormal after pymc 5.17. If you see above i have a hierarchical SEM model which uses an indexing trick to pass group specific covariance structures to the likelihood. While this works in 5.17 see above it breaks in 5.30... is that a bug, or intended behaviour. Do you know how i could replicate the results with 5.30+?

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

ricardoV94 commented on 2025-09-30T10:28:00Z
----------------------------------------------------------------

Line #9.    corr_values = [

Nit this is a terrible way to have the values in the notebook for a reader. Just tell black to ignore and let multiple values per line

[ricardoV94]

@NathanielF I didn't see any block failing in the notebook, can you give me a small snippet of code that is failing for you? You showed numpy code above

[NathanielF]

The notebook code works but its running on pymc 5.17, you can see in the watermark... if i change or update the version. I tried to 5.30 it breaks on cell which creates the hierarchical modelling and gives the trace back you see above. I can run the notebook tonight on 5.30 and push it to show you

But it should break on the cell that defines the hierarchical model. Just under the section heading "hierarchical model on structural components"...

[ricardoV94]

I see, let me take a quick look

[NathanielF]

Thank you!

[ricardoV94]

What do you mean by pymc 5.30, last release is 5.25

[ricardoV94]

Unrelated but please don't do this pytensor.config.cxx = "/usr/bin/clang++". It's specific to your setup, probably macos, but it means that someone with a different config is going to loose all their C caching.

You can use a .pytensorrc in your home directory instead: https://pytensor.readthedocs.io/en/latest/library/config.html

[NathanielF]

What do you mean by pymc 5.30, last release is 5.25

Oh, shoot, sorry i think this is my fault. I was using 5.3.0. I thought because the pymc-examples pixi installer had the > 5.16... it was 5.30.

Sorry, my bad

[ricardoV94]

No worries, problem goes away then?

Suggestion use the newer syntax for set_subtensor:

Lambda = pt.set_subtensor(Lambda[0:3, 0], lambdas_1)
# Equivalent
Lambda = Lambda[0:3, 0].set(lambdas_1)

[NathanielF]

Spent a good bit of time tightening this write up over the weekend. It should be in a good place for review now @fonnesbeck

Thanks!

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:56Z
----------------------------------------------------------------

Suggestion: Maybe explain to the uninitiated why the first loading is fixed to 1

NathanielF commented on 2025-10-21T09:21:41Z
----------------------------------------------------------------

Added some detail in a comment on the code. There is also a brief discussion in the text below.

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:57Z
----------------------------------------------------------------

Line #74.        r = pt.set_subtensor(r[3, 5], beta_params[0])

Hard coding indicies always makes me nervous. Maybe create global variables UFI=3 and FOR=5 and use them?

NathanielF commented on 2025-10-21T09:22:31Z
----------------------------------------------------------------

changed this in both the Lambda and the Psi. Used a lookup based off the coordinates. Should be safer and hopefully more meaningful.

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:57Z
----------------------------------------------------------------

Isn't the pseudo-observation matrix just conditional means?

NathanielF commented on 2025-10-21T09:23:03Z
----------------------------------------------------------------

Yes basically for each observation. Added a sentence stating that.

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:58Z
----------------------------------------------------------------

ussually -> usually

NathanielF commented on 2025-10-21T09:23:40Z
----------------------------------------------------------------

fixed.

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:59Z
----------------------------------------------------------------

proprely -> properly

NathanielF commented on 2025-10-21T09:23:53Z
----------------------------------------------------------------

fixed

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:15:59Z
----------------------------------------------------------------

observe is -> observe its

inn -> in

NathanielF commented on 2025-10-21T09:24:11Z
----------------------------------------------------------------

fixed. Thanks

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:16:00Z
----------------------------------------------------------------

we can happy -> we can be happy (?)

NathanielF commented on 2025-10-21T09:25:34Z
----------------------------------------------------------------

adjusted this sentence. "Happy" didn't seem right in retrospect

[review-notebook-app]

View / edit / reply to this conversation on ReviewNB

fonnesbeck commented on 2025-10-20T22:16:01Z
----------------------------------------------------------------

two-step of -> two step process of (?)

"deriving an interpretable representation of the latent attitude in their expressions" is awkward. Suggestion: "deriving interpretable representations of latent attitudes from expressed opinions"

NathanielF commented on 2025-10-21T09:26:14Z
----------------------------------------------------------------

fixed. Thanks

[NathanielF]

Added some detail in a comment on the code. There is also a brief discussion in the text below.

---

View entire conversation on ReviewNB

[NathanielF]

changed this in both the Lambda and the Psi. Used a lookup based off the coordinates. Should be safer and hopefully more meaningful.

---

View entire conversation on ReviewNB

[NathanielF]

Yes basically for each observation. Added a sentence stating that.

---

View entire conversation on ReviewNB

[NathanielF]

fixed.

---

View entire conversation on ReviewNB

[NathanielF]

fixed

---

View entire conversation on ReviewNB

[NathanielF]

fixed. Thanks

---

View entire conversation on ReviewNB

[NathanielF]

adjusted this sentence. "Happy" didn't seem right in retrospect

---

View entire conversation on ReviewNB

[NathanielF]

fixed. Thanks

---

View entire conversation on ReviewNB

[NathanielF]

Thanks @fonnesbeck ! Addressed those comments.