Hierarchical Structural Equation Models in PyMC

[NathanielF]

Notebook proposal

Title: Hierarchical Structural Equation Models and Contemporary Bayesian Workflows

Why should this notebook be added to pymc-examples?

We have one example of a CFA and SEM and i think this is good as far as it goes, but it veers primarily pedagogical and slowly introduces the models with a slightly cumbersome implementation that focuses on (a) a Normal likelihood presupposing independence in the likelihood and relying on latent covariance structures to recover correlation structure and (b) doesn't make use of the recursive solving of mediation relationships among the latent variables. I would like to demonstrate how we move beyond the more pedagogical implementation, to a production ready specification that includes multivariate normal likelihoods, recursive solving of the direct and indirect effects and hierarchical structure to take advantage of the Bayesian setting. I'll also align the notation used in the model implementation. Since the other notebook is also quite long, i propose a supplementary one.

Suggested categories:

• Level: Advanced
• Diataxis type:

Related notebooks

https://www.pymc.io/projects/examples/en/latest/case_studies/CFA_SEM.html

with pm.Model(coords=coords) as sem_model_v3:
    
    # --- Factor loadings ---
    lambdas_1 = make_lambda('indicators_1', 'lambdas1', priors=[1, .5])
    lambdas_2 = make_lambda('indicators_2', 'lambdas2', priors=[1, .5])
    lambdas_3 = make_lambda('indicators_3', 'lambdas3', priors=[1, .5])
    lambdas_4 = make_lambda('indicators_4', 'lambdas4', priors=[1, .5])

    Lambda = pt.zeros((12, 4))
    Lambda = pt.set_subtensor(Lambda[0:3, 0], lambdas_1)
    Lambda = pt.set_subtensor(Lambda[3:6, 1], lambdas_2)
    Lambda = pt.set_subtensor(Lambda[6:9, 2], lambdas_3)
    Lambda = pt.set_subtensor(Lambda[9:12, 3], lambdas_4)
    Lambda = pm.Deterministic('Lambda', Lambda)

    latent_dim = len(coords['latent'])

    sd_dist = pm.Exponential.dist(1.0, shape=latent_dim)
    chol, _, _ = pm.LKJCholeskyCov("chol_cov", n=latent_dim, eta=2, sd_dist=sd_dist, compute_corr=True)
    gamma = pm.MvNormal("gamma", 0, chol=chol, dims=("obs", "latent"))

    B = make_B()
    I = pt.eye(latent_dim)
    eta = pm.Deterministic("eta", pt.slinalg.solve(I - B + 1e-8 * I, gamma.T).T)  

    mu = pt.dot(eta, Lambda.T)

    ## Error Terms
    Psi = make_Psi('indicators')
    _ = pm.MvNormal('likelihood', mu=mu, cov=Psi, observed=observed_data)

References

If applicable, references and material that could help in writing the notebook.

[NathanielF]

@fonnesbeck might be interested.