simpsons-paradox

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.13%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22arviz%3D%3D0.23.4%22%2C%0A%23%20%20%20%20%20%22graphviz%3D%3D0.21%22%2C%0A%23%20%20%20%20%20%22marimo%3E%3D0.20.4%22%2C%0A%23%20%20%20%20%20%22matplotlib%3D%3D3.10.8%22%2C%0A%23%20%20%20%20%20%22numpy%3D%3D2.4.2%22%2C%0A%23%20%20%20%20%20%22pandas%3D%3D3.0.1%22%2C%0A%23%20%20%20%20%20%22seaborn%3D%3D0.13.2%22%2C%0A%23%20%20%20%20%20%22xarray%3D%3D2026.2.0%22%2C%0A%23%20%5D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.1%22%0Aapp%20%3D%20marimo.App(auto_download%3D%5B%22html%22%5D)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell%0Aasync%20def%20_()%3A%0A%20%20%20%20%23%20From%20github.com%2Fmarimo-team%2Fgallery-examples%20notebooks%2Fanalysis%2Fsimpsons-paradox.py%20(MIT)%2C%0A%20%20%20%20%23%20changed%20to%20sample%20in%20the%20page%20with%20pymcwasm.%0A%20%20%20%20%23%20numba%20has%20no%20wasm%20build%2C%20so%20pytensor%20and%20pymc%20install%20without%20their%20declared%20deps.%0A%20%20%20%20import%20micropip%0A%0A%20%20%20%20SITE%20%3D%20%22https%3A%2F%2Fhabakan.github.io%2Fpymcwasm%2Flite%2Ffiles%22%0A%20%20%20%20await%20micropip.install(%5B%22setuptools%22%2C%20%22filelock%22%2C%20%22cachetools%22%2C%20%22typing-extensions%22%2C%20%22cloudpickle%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22cons%22%2C%20%22logical-unification%22%2C%20%22etuples%22%2C%20%22miniKanren%22%2C%20%22rich%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22threadpoolctl%22%2C%20%22scipy%22%5D)%0A%20%20%20%20%23%20marimo%20installs%20the%20newest%20arviz%20from%20the%20imports%3B%20this%20notebook%20is%20written%20for%200.23.%0A%20%20%20%20await%20micropip.install(%22arviz%3D%3D0.23.4%22%2C%20reinstall%3DTrue)%0A%20%20%20%20await%20micropip.install(%5B%22pytensor%3D%3D2.38.3%22%2C%20%22pymc%3D%3D5.28.1%22%5D%2C%20deps%3DFalse)%0A%20%20%20%20await%20micropip.install(SITE%20%2B%20%22%2Fpymcwasm-0.1.0-py3-none-any.whl%22%2C%20deps%3DFalse)%0A%0A%20%20%20%20import%20pytensor%0A%20%20%20%20pytensor.config.mode%20%3D%20%22FAST_COMPILE%22%0A%20%20%20%20pytensor.config.linker%20%3D%20%22py%22%0A%20%20%20%20import%20pymcwasm.linker%0A%20%20%20%20await%20pymcwasm.linker.load(SITE%20%2B%20%22%2Ftapewasm%2Fpkg%2Ftapewasm.js%22)%0A%20%20%20%20wasm_ready%20%3D%20True%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Aasync%20def%20_(mo)%3A%0A%20%20%20%20%23%20graphviz%20renders%20by%20running%20%60dot%60%2C%20which%20the%20page%20cannot%3B%20viz.js%20is%20the%20same%20Graphviz%20in%20wasm.%0A%20%20%20%20from%20pyodide.code%20import%20run_js%0A%0A%20%20%20%20_svg%20%3D%20await%20run_js(%22%22%22import(%22https%3A%2F%2Fcdn.jsdelivr.net%2Fnpm%2F%40viz-js%2Fviz%403.11.0%2F%2Besm%22)%0A%20%20%20%20%20%20.then(m%20%3D%3E%20m.instance()).then(viz%20%3D%3E%20src%20%3D%3E%20viz.renderString(src%2C%20%7Bformat%3A%20%22svg%22%7D))%22%22%22)%0A%0A%20%20%20%20def%20dot(g)%3A%0A%20%20%20%20%20%20%20%20return%20mo.Html(_svg(g.source))%0A%0A%20%20%20%20return%20(dot%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Simpson's%20paradox%0A%0A%20%20%20%20%3E%20This%20notebook%20was%20originally%20found%20on%20the%20PyMC%20docs%20%5Bhere%5D(https%3A%2F%2Fwww.pymc.io%2Fprojects%2Fexamples%2Fen%2Flatest%2Fcausal_inference%2FGLM-simpsons-paradox.html)%20written%20by%20Benjamin%20T.%20Vincent%2C%20%5BDOI%5D(https%3A%2F%2Fzenodo.org%2Frecords%2F18677308).%20It%20was%20adapted%20for%20marimo.%0A%20%20%20%20%3E%20This%20copy%20of%20%5Bmarimo's%20gallery%20version%5D(https%3A%2F%2Fgithub.com%2Fmarimo-team%2Fgallery-examples%2Fblob%2Fmain%2Fnotebooks%2Fanalysis%2Fsimpsons-paradox.py)%20samples%20in%20the%20page%2C%20with%20PyTensor's%20%60mode%3D%22WASM%22%60%20from%20%5Bpymcwasm%5D(https%3A%2F%2Fgithub.com%2Fhabakan%2Fpymcwasm).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%5BSimpson's%20Paradox%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSimpson%2527s_paradox)%20describes%20a%20situation%20where%20there%20might%20be%20a%20negative%20relationship%20between%20two%20variables%20within%20a%20group%2C%20but%20when%20data%20from%20multiple%20groups%20are%20combined%2C%20that%20relationship%20may%20disappear%20or%20even%20reverse%20sign.%20The%20gif%20below%20(from%20the%20Simpson's%20Paradox%20%5BWikipedia%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSimpson%2527s_paradox)%20page)%20demonstrates%20this%20very%20nicely.%0A%0A%20%20%20%20!%5B%5D(https%3A%2F%2Fupload.wikimedia.org%2Fwikipedia%2Fcommons%2Ff%2Ffb%2FSimpsons_paradox_-_animation.gif)%0A%0A%20%20%20%20Another%20way%20of%20describing%20this%20is%20that%20we%20wish%20to%20estimate%20the%20causal%20relationship%20%24x%20%5Crightarrow%20y%24.%20The%20seemingly%20obvious%20approach%20of%20modelling%20%60y%20~%201%20%2B%20x%60%20will%20lead%20us%20to%20conclude%20(in%20the%20situation%20above)%20that%20increasing%20%24x%24%20causes%20%24y%24%20to%20decrease%20(see%20Model%201%20below).%20However%2C%20the%20relationship%20between%20%24x%24%20and%20%24y%24%20is%20confounded%20by%20a%20group%20membership%20variable%20%24group%24.%20This%20group%20membership%20variable%20is%20not%20included%20in%20the%20model%2C%20and%20so%20the%20relationship%20between%20%24x%24%20and%20%24y%24%20is%20biased.%20If%20we%20now%20factor%20in%20the%20influence%20of%20%24group%24%2C%20in%20some%20situations%20(e.g.%20the%20image%20above)%20this%20can%20lead%20us%20to%20completely%20reverse%20the%20sign%20of%20our%20estimate%20of%20%24x%20%5Crightarrow%20y%24%2C%20now%20estimating%20that%20increasing%20%24x%24%20causes%20%24y%24%20to%20_increase_.%0A%0A%20%20%20%20In%20short%2C%20this%20'paradox'%20(or%20simply%20omitted%20variable%20bias)%20can%20be%20resolved%20by%20assuming%20a%20causal%20DAG%20which%20includes%20how%20the%20main%20predictor%20variable%20_and_%20group%20membership%20(the%20confounding%20variable)%20influence%20the%20outcome%20variable.%20We%20demonstrate%20an%20example%20where%20we%20_don't_%20incorporate%20group%20membership%20(so%20our%20causal%20DAG%20is%20wrong%2C%20or%20in%20other%20words%20our%20model%20is%20misspecified%3B%20Model%201).%20We%20then%20show%202%20ways%20to%20resolve%20this%20by%20including%20group%20membership%20as%20causal%20influence%20upon%20%24x%24%20and%20%24y%24.%20This%20is%20shown%20in%20an%20unpooled%20model%20(Model%202)%20and%20a%20hierarchical%20model%20(Model%203).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20arviz%20as%20az%0A%20%20%20%20import%20graphviz%20as%20gr%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20pymc%20as%20pm%0A%20%20%20%20import%20seaborn%20as%20sns%0A%20%20%20%20import%20xarray%20as%20xr%0A%0A%20%20%20%20return%20az%2C%20gr%2C%20np%2C%20pd%2C%20plt%2C%20pm%2C%20sns%2C%20xr%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20np%2C%20plt)%3A%0A%20%20%20%20az.style.use(%22arviz-darkgrid%22)%0A%20%20%20%20figsize%20%3D%20%5B12%2C%204%5D%0A%20%20%20%20plt.rcParams%5B%22figure.figsize%22%5D%20%3D%20figsize%0A%20%20%20%20rng%20%3D%20np.random.default_rng(1234)%0A%20%20%20%20return%20figsize%2C%20rng%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Generate%20data%0A%0A%20%20%20%20This%20data%20generation%20was%20influenced%20by%20this%20%5Bstackexchange%5D(https%3A%2F%2Fstats.stackexchange.com%2Fquestions%2F479201%2Funderstanding-simpsons-paradox-with-random-effects)%20question.%20It%20will%20comprise%20observations%20from%20%24G%3D5%24%20groups.%20The%20data%20is%20constructed%20such%20that%20there%20is%20a%20negative%20relationship%20between%20%24x%24%20and%20%24y%24%20within%20each%20group%2C%20but%20when%20all%20groups%20are%20combined%2C%20the%20relationship%20is%20positive.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20rng)%3A%0A%20%20%20%20def%20generate()%3A%0A%20%20%20%20%20%20%20%20group_list%20%3D%20%5B'one'%2C%20'two'%2C%20'three'%2C%20'four'%2C%20'five'%5D%0A%20%20%20%20%20%20%20%20trials_per_group%20%3D%2020%0A%20%20%20%20%20%20%20%20group_intercepts%20%3D%20rng.normal(0%2C%201%2C%20len(group_list))%0A%20%20%20%20%20%20%20%20group_slopes%20%3D%20np.ones(len(group_list))%20*%20-0.5%0A%20%20%20%20%20%20%20%20group_mx%20%3D%20group_intercepts%20*%202%0A%20%20%20%20%20%20%20%20group%20%3D%20np.repeat(group_list%2C%20trials_per_group)%0A%20%20%20%20%20%20%20%20subject%20%3D%20np.concatenate(%5Bnp.ones(trials_per_group)%20*%20i%20for%20i%20in%20np.arange(len(group_list))%5D).astype(int)%0A%20%20%20%20%20%20%20%20intercept%20%3D%20np.repeat(group_intercepts%2C%20trials_per_group)%0A%20%20%20%20%20%20%20%20slope%20%3D%20np.repeat(group_slopes%2C%20trials_per_group)%0A%20%20%20%20%20%20%20%20mx%20%3D%20np.repeat(group_mx%2C%20trials_per_group)%0A%20%20%20%20%20%20%20%20_x%20%3D%20rng.normal(mx%2C%201)%0A%20%20%20%20%20%20%20%20y%20%3D%20rng.normal(intercept%20%2B%20(_x%20-%20mx)%20*%20slope%2C%201)%0A%20%20%20%20%20%20%20%20data%20%3D%20pd.DataFrame(%7B'group'%3A%20group%2C%20'group_idx'%3A%20subject%2C%20'x'%3A%20_x%2C%20'y'%3A%20y%7D)%0A%20%20%20%20%20%20%20%20return%20(data%2C%20group_list)%0A%20%20%20%20data%2C%20group_list%20%3D%20generate()%0A%20%20%20%20return%20data%2C%20group_list%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20To%20follow%20along%2C%20it%20is%20useful%20to%20clearly%20understand%20the%20form%20of%20the%20data.%20This%20is%20%5Blong%20form%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FWide_and_narrow_data)%20data%20(also%20known%20as%20narrow%20data)%20in%20that%20each%20row%20represents%20one%20observation.%20We%20have%20a%20%60group%60%20column%20which%20has%20the%20group%20label%2C%20and%20an%20accompanying%20numerical%20%60group_idx%60%20column.%20This%20is%20very%20useful%20when%20it%20comes%20to%20modelling%20as%20we%20can%20use%20it%20as%20an%20index%20to%20look%20up%20group-level%20parameter%20estimates.%20Finally%2C%20we%20have%20our%20core%20observations%20of%20the%20predictor%20variable%20%60x%60%20and%20the%20outcome%20%60y%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data)%3A%0A%20%20%20%20data%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20And%20we%20can%20visualise%20this%20as%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20plt%2C%20sns)%3A%0A%20%20%20%20fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%206))%0A%20%20%20%20sns.scatterplot(data%3Ddata%2C%20x%3D'x'%2C%20y%3D'y'%2C%20hue%3D'group'%2C%20ax%3D_ax)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20rest%20of%20the%20notebook%20will%20cover%20different%20ways%20that%20we%20can%20analyse%20this%20data%20using%20linear%20models.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Model%201%3A%20Pooled%20regression%0A%0A%20%20%20%20First%20we%20examine%20the%20simplest%20model%20-%20plain%20linear%20regression%20which%20pools%20all%20the%20data%20and%20has%20no%20knowledge%20of%20the%20group%2Fmulti-level%20structure%20of%20the%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20From%20a%20causal%20perspective%2C%20this%20approach%20embodies%20the%20belief%20that%20%24x%24%20causes%20%24y%24%20and%20that%20this%20relationship%20is%20constant%20across%20all%20groups%2C%20or%20groups%20are%20simply%20not%20considered.%20This%20can%20be%20shown%20in%20the%20causal%20DAG%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dot%2C%20gr)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20g%20%3D%20gr.Digraph()%0A%20%20%20%20g.node(name%3D%22x%22%2C%20label%3D%22x%22)%0A%20%20%20%20g.node(name%3D%22y%22%2C%20label%3D%22y%22)%0A%20%20%20%20g.edge(tail_name%3D%22x%22%2C%20head_name%3D%22y%22)%0A%20%20%20%20dot(g)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20could%20describe%20this%20model%20mathematically%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cbeta_0%2C%20%5Cbeta_1%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%205)%20%5C%5C%0A%20%20%20%20%5Csigma%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%5Cmu_i%20%26%3D%20%5Cbeta_0%20%2B%20%5Cbeta_1%20x_i%20%5C%5C%0A%20%20%20%20y_i%20%26%5Csim%20%5Ctext%7BNormal%7D(%5Cmu_i%2C%20%5Csigma)%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20can%20also%20express%20Model%201%20in%20Wilkinson%20notation%20as%20%60y%20~%201%20%2B%20x%60%20which%20is%20equivalent%20to%20%60y%20~%20x%60%20as%20the%20intercept%20is%20included%20by%20default.%0A%0A%20%20%20%20*%20The%20%601%60%20term%20corresponds%20to%20the%20intercept%20term%20%24%5Cbeta_0%24.%0A%20%20%20%20*%20The%20%60x%60%20term%20corresponds%20to%20the%20slope%20term%20%24%5Cbeta_1%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20So%20now%20we%20can%20express%20this%20as%20a%20PyMC%20model.%20We%20can%20notice%20how%20closely%20the%20model%20syntax%20mirrors%20the%20mathematical%20notation%20above.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20pm)%3A%0A%20%20%20%20with%20pm.Model()%20as%20model1%3A%0A%20%20%20%20%20%20%20%20_%CE%B20%20%3D%20pm.Normal('%CE%B20'%2C%200%2C%20sigma%3D5)%0A%20%20%20%20%20%20%20%20_%CE%B21%20%3D%20pm.Normal('%CE%B21'%2C%200%2C%20sigma%3D5)%0A%20%20%20%20%20%20%20%20_sigma%20%3D%20pm.Gamma('sigma'%2C%202%2C%202)%0A%20%20%20%20%20%20%20%20_x%20%3D%20pm.Data('x'%2C%20data.x%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20_%CE%BC%20%3D%20pm.Deterministic('%CE%BC'%2C%20_%CE%B20%20%2B%20_%CE%B21%20*%20_x%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20pm.Normal('y'%2C%20mu%3D_%CE%BC%2C%20sigma%3D_sigma%2C%20observed%3Ddata.y%2C%20dims%3D'obs_id')%0A%20%20%20%20return%20(model1%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20And%20we%20can%20visualize%20the%20DAG%20which%20can%20be%20a%20useful%20way%20to%20check%20that%20our%20model%20is%20correctly%20specified.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dot%2C%20model1%2C%20pm)%3A%0A%20%20%20%20dot(pm.model_to_graphviz(model1))%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Conduct%20inference%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(model1%2C%20pm%2C%20rng)%3A%0A%20%20%20%20with%20model1%3A%0A%20%20%20%20%20%20%20%20idata1%20%3D%20pm.sample(random_seed%3Drng%2C%20cores%3D1%2C%20compile_kwargs%3D%7B%22mode%22%3A%20%22WASM%22%7D)%0A%20%20%20%20return%20(idata1%2C)%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20idata1)%3A%0A%20%20%20%20az.plot_trace(idata1%2C%20var_names%3D%5B%22~%CE%BC%22%5D)%3B%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Visualisation%0A%0A%20%20%20%20First%20we'll%20define%20a%20handy%20predict%20function%20which%20will%20do%20out%20of%20sample%20predictions%20for%20us.%20This%20will%20be%20handy%20when%20it%20comes%20to%20visualising%20the%20model%20fits.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20pm%2C%20rng)%3A%0A%20%20%20%20def%20predict(model%3A%20pm.Model%2C%20idata%3A%20az.InferenceData%2C%20predict_at%3A%20dict)%20-%3E%20az.InferenceData%3A%0A%20%20%20%20%20%20%20%20%22%22%22Do%20posterior%20predictive%20inference%20at%20a%20set%20of%20out%20of%20sample%20points%20specified%20by%20%60predict_at%60.%22%22%22%0A%20%20%20%20%20%20%20%20with%20model%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pm.set_data(predict_at)%0A%20%20%20%20%20%20%20%20%20%20%20%20idata.extend(pm.sample_posterior_predictive(idata%2C%20var_names%3D%5B%22y%22%2C%20%22%CE%BC%22%5D%2C%20random_seed%3Drng))%0A%20%20%20%20%20%20%20%20return%20idata%0A%0A%20%20%20%20return%20(predict%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20And%20now%20let's%20use%20that%20%60predict%60%20function%20to%20do%20out%20of%20sample%20predictions%20which%20we%20will%20use%20for%20visualisation.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20idata1%2C%20model1%2C%20np%2C%20predict)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-output%0A%20%20%20%20xi%20%3D%20np.linspace(data.x.min()%2C%20data.x.max()%2C%2020)%0A%20%20%20%20idata1_1%20%3D%20predict(model%3Dmodel1%2C%20idata%3Didata1%2C%20predict_at%3D%7B'x'%3A%20xi%7D)%0A%20%20%20%20return%20idata1_1%2C%20xi%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Finally%2C%20we%20can%20now%20visualise%20the%20model%20fit%20to%20data%2C%20and%20our%20posterior%20in%20parameter%20space.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20data%2C%20idata1_1%2C%20plt%2C%20xi%2C%20xr)%3A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%206))%0A%0A%20%20%20%20def%20plot_band(xi%2C%20var%3A%20xr.DataArray%2C%20ax%2C%20color%3A%20str)%3A%0A%20%20%20%20%20%20%20%20_ax.plot(xi%2C%20var.mean(%5B'chain'%2C%20'draw'%5D)%2C%20color%3Dcolor)%0A%20%20%20%20%20%20%20%20az.plot_hdi(xi%2C%20var%2C%20hdi_prob%3D0.6%2C%20color%3Dcolor%2C%20fill_kwargs%3D%7B'alpha'%3A%200.2%2C%20'linewidth'%3A%200%7D%2C%20ax%3D_ax)%0A%20%20%20%20%20%20%20%20az.plot_hdi(xi%2C%20var%2C%20hdi_prob%3D0.95%2C%20color%3Dcolor%2C%20fill_kwargs%3D%7B'alpha'%3A%200.2%2C%20'linewidth'%3A%200%7D%2C%20ax%3D_ax)%0A%0A%20%20%20%20def%20_plot(idata%3A%20az.InferenceData)%3A%0A%20%20%20%20%20%20%20%20fig%2C%20_ax%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12%2C%204))%0A%20%20%20%20%20%20%20%20_ax%5B0%5D.scatter(data.x%2C%20data.y%2C%20color%3D'k')%0A%20%20%20%20%20%20%20%20plot_band(xi%2C%20idata.posterior_predictive.%CE%BC%2C%20ax%3D_ax%5B0%5D%2C%20color%3D'k')%0A%20%20%20%20%20%20%20%20_ax%5B0%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Conditional%20mean')%0A%20%20%20%20%20%20%20%20_ax%5B1%5D.scatter(data.x%2C%20data.y%2C%20color%3D'k')%0A%20%20%20%20%20%20%20%20plot_band(xi%2C%20idata.posterior_predictive.y%2C%20ax%3D_ax%5B1%5D%2C%20color%3D'k')%0A%20%20%20%20%20%20%20%20_ax%5B1%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Posterior%20predictive%20distribution')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.scatter(az.extract(idata%2C%20var_names%3D%5B'%CE%B21'%5D)%2C%20az.extract(idata%2C%20var_names%3D%5B'%CE%B20'%5D)%2C%20color%3D'k'%2C%20alpha%3D0.01%2C%20rasterized%3DTrue)%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.set(xlabel%3D'slope'%2C%20ylabel%3D'intercept'%2C%20title%3D'Parameter%20space')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axhline(y%3D0%2C%20c%3D'k')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axvline(x%3D0%2C%20c%3D'k')%0A%0A%0A%20%20%20%20_plot(idata1_1)%0A%20%20%20%20_ax%0A%20%20%20%20return%20(plot_band%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20plot%20on%20the%20left%20shows%20the%20data%20and%20the%20posterior%20of%20the%20**conditional%20mean**.%20For%20a%20given%20%24x%24%2C%20we%20get%20a%20posterior%20distribution%20of%20the%20model%20(i.e.%20of%20%24%5Cmu%24).%0A%0A%20%20%20%20The%20plot%20in%20the%20middle%20shows%20the%20conditional%20**posterior%20predictive%20distribution**%2C%20which%20gives%20a%20statement%20about%20the%20data%20we%20expect%20to%20see.%20Intuitively%2C%20this%20can%20be%20understood%20as%20not%20only%20incorporating%20what%20we%20know%20of%20the%20model%20(left%20plot)%20but%20also%20what%20we%20know%20about%20the%20distribution%20of%20error.%0A%0A%20%20%20%20The%20plot%20on%20the%20right%20shows%20our%20posterior%20beliefs%20in%20**parameter%20space**.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20One%20of%20the%20clear%20things%20about%20this%20analysis%20is%20that%20we%20have%20credible%20evidence%20that%20%24x%24%20and%20%24y%24%20are%20_positively_%20correlated.%20We%20can%20see%20this%20from%20the%20posterior%20over%20the%20slope%20(see%20right%20hand%20panel%20in%20the%20figure%20above)%20which%20we%20isolate%20in%20the%20plot%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20idata1_1)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20_ax%20%3D%20az.plot_posterior(idata1_1.posterior%5B'%CE%B21'%5D%2C%20ref_val%3D0)%0A%20%20%20%20_ax.set(title%3D'Model%201%20strongly%20suggests%20a%20positive%20slope'%2C%20xlabel%3D'%24%5C%5Cbeta_1%24')%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Model%202%3A%20Unpooled%20regression%20with%20counfounder%20included%0A%0A%20%20%20%20We%20will%20use%20the%20same%20data%20in%20this%20analysis%2C%20but%20this%20time%20we%20will%20use%20our%20knowledge%20that%20data%20come%20from%20groups.%20From%20a%20causal%20perspective%20we%20are%20exploring%20the%20notion%20that%20both%20%24x%24%20and%20%24y%24%20are%20influenced%20by%20group%20membership.%20This%20can%20be%20shown%20in%20the%20causal%20directed%20acyclic%20graph%20(%5BDAG%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FDirected_acyclic_graph))%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dot%2C%20gr)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20g_1%20%3D%20gr.Digraph()%0A%20%20%20%20g_1.node(name%3D'x'%2C%20label%3D'x')%0A%20%20%20%20g_1.node(name%3D'g'%2C%20label%3D'group')%0A%20%20%20%20g_1.node(name%3D'y'%2C%20label%3D'y')%0A%20%20%20%20g_1.edge(tail_name%3D'x'%2C%20head_name%3D'y')%0A%20%20%20%20g_1.edge(tail_name%3D'g'%2C%20head_name%3D'x')%0A%20%20%20%20g_1.edge(tail_name%3D'g'%2C%20head_name%3D'y')%0A%20%20%20%20dot(g_1)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20So%20we%20can%20see%20that%20%24group%24%20is%20a%20%5Bconfounding%20variable%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FConfounding).%20So%20if%20we%20are%20trying%20to%20discover%20the%20causal%20relationship%20of%20%24x%24%20on%20%24y%24%2C%20we%20need%20to%20account%20for%20the%20confounding%20variable%20%24group%24.%20Model%201%20did%20not%20do%20this%20and%20so%20arrived%20at%20the%20wrong%20conclusion.%20But%20Model%202%20will%20account%20for%20this.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20More%20specifically%20we%20will%20essentially%20fit%20independent%20regressions%20to%20data%20within%20each%20group.%20This%20could%20also%20be%20described%20as%20an%20unpooled%20model.%20We%20could%20describe%20this%20model%20mathematically%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cvec%7B%5Cbeta_0%7D%2C%20%5Cvec%7B%5Cbeta_1%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%205)%20%5C%5C%0A%20%20%20%20%5Csigma%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%5Cmu_i%20%26%3D%20%5Cvec%7B%5Cbeta_0%7D%5Bg_i%5D%20%2B%20%5Cvec%7B%5Cbeta_1%7D%5Bg_i%5D%20x_i%20%5C%5C%0A%20%20%20%20y_i%20%26%5Csim%20%5Ctext%7BNormal%7D(%5Cmu_i%2C%20g_i)%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Where%20%24g_i%24%20is%20the%20group%20index%20for%20observation%20%24i%24.%20So%20the%20parameters%20%24%5Cvec%7B%5Cbeta_0%7D%24%20and%20%24%5Cvec%7B%5Cbeta_1%7D%24%20are%20now%20length%20%24G%24%20vectors%2C%20not%20scalars.%20And%20the%20%24%5Bg_i%5D%24%20acts%20as%20an%20index%20to%20look%20up%20the%20group%20for%20the%20%24i%5E%5Ctext%7Bth%7D%24%20observation.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3A%3A%3A%7Bnote%7D%0A%20%20%20%20We%20can%20also%20express%20this%20Model%202%20in%20Wilkinson%20notation%20as%20%60y%20~%200%20%2B%20g%20%2B%20x%3Ag%60.%0A%0A%20%20%20%20*%20The%20%60g%60%20term%20captures%20the%20group%20specific%20intercept%20%24%5Cbeta_0%5Bg_i%5D%24%20parameters.%0A%20%20%20%20*%20The%20%600%60%20means%20we%20do%20not%20have%20a%20global%20intercept%20term%2C%20leaving%20the%20group%20specific%20intercepts%20to%20be%20the%20only%20intercepts.%0A%20%20%20%20*%20The%20%60x%3Ag%60%20term%20captures%20group%20specific%20slope%20%24%5Cbeta_1%5Bg_i%5D%24%20parameters.%0A%20%20%20%20%3A%3A%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Let's%20express%20Model%202%20with%20PyMC%20code.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20group_list%2C%20pm)%3A%0A%20%20%20%20coords%20%3D%20%7B'group'%3A%20group_list%7D%0A%20%20%20%20with%20pm.Model(coords%3Dcoords)%20as%20model2%3A%0A%20%20%20%20%20%20%20%20_%CE%B20%20%3D%20pm.Normal('%CE%B20'%2C%200%2C%20sigma%3D5%2C%20dims%3D'group')%0A%20%20%20%20%20%20%20%20_%CE%B21%20%3D%20pm.Normal('%CE%B21'%2C%200%2C%20sigma%3D5%2C%20dims%3D'group')%0A%20%20%20%20%20%20%20%20_sigma%20%3D%20pm.Gamma('sigma'%2C%202%2C%202)%0A%20%20%20%20%20%20%20%20_x%20%3D%20pm.Data('x'%2C%20data.x%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20g_2%20%3D%20pm.Data('g'%2C%20data.group_idx%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20_%CE%BC%20%3D%20pm.Deterministic('%CE%BC'%2C%20_%CE%B20%5Bg_2%5D%20%2B%20_%CE%B21%5Bg_2%5D%20*%20_x%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20pm.Normal('y'%2C%20mu%3D_%CE%BC%2C%20sigma%3D_sigma%2C%20observed%3Ddata.y%2C%20dims%3D'obs_id')%0A%20%20%20%20return%20coords%2C%20model2%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20By%20plotting%20the%20DAG%20for%20this%20model%20it%20is%20clear%20to%20see%20that%20we%20now%20have%20individual%20intercept%20and%20slope%20parameters%20for%20each%20of%20the%20groups.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dot%2C%20model2%2C%20pm)%3A%0A%20%20%20%20dot(pm.model_to_graphviz(model2))%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Conduct%20inference%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20model2%2C%20pm%2C%20rng)%3A%0A%20%20%20%20with%20model2%3A%0A%20%20%20%20%20%20%20%20idata2%20%3D%20pm.sample(random_seed%3Drng%2C%20cores%3D1%2C%20compile_kwargs%3D%7B%22mode%22%3A%20%22WASM%22%7D)%0A%0A%20%20%20%20az.plot_trace(idata2%2C%20var_names%3D%5B%22~%CE%BC%22%5D)%3B%0A%20%20%20%20return%20(idata2%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Visualisation%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20np)%3A%0A%20%20%20%20_n_points%20%3D%2010%0A%20%20%20%20_n_groups%20%3D%20len(data.group.unique())%0A%20%20%20%20xi_1%20%3D%20np.concatenate(%5Bnp.linspace(group%5B1%5D.x.min()%2C%20group%5B1%5D.x.max()%2C%20_n_points)%20for%20group%20in%20data.groupby('group_idx')%5D)%0A%20%20%20%20g_3%20%3D%20np.concatenate(%5B%5Bi%5D%20*%20_n_points%20for%20i%20in%20range(_n_groups)%5D).astype(int)%0A%20%20%20%20predict_at%20%3D%20%7B'x'%3A%20xi_1%2C%20'g'%3A%20g_3%7D%0A%20%20%20%20return%20g_3%2C%20predict_at%2C%20xi_1%0A%0A%0A%40app.cell%0Adef%20_(idata2%2C%20model2%2C%20predict%2C%20predict_at)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-output%0A%20%20%20%20idata2_1%20%3D%20predict(model%3Dmodel2%2C%20idata%3Didata2%2C%20predict_at%3Dpredict_at)%0A%20%20%20%20return%20(idata2_1%2C)%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20data%2C%20g_3%2C%20group_list%2C%20idata2_1%2C%20plot_band%2C%20plt%2C%20xi_1)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20def%20_plot(idata)%3A%0A%20%20%20%20%20%20%20%20fig%2C%20_ax%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12%2C%204))%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(len(group_list))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B0%5D.scatter(data.x%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20data.y%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20plot_band(xi_1%5Bg_3%20%3D%3D%20i%5D%2C%20idata.posterior_predictive.%CE%BC.isel(obs_id%3Dg_3%20%3D%3D%20i)%2C%20ax%3D_ax%5B0%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B1%5D.scatter(data.x%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20data.y%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20plot_band(xi_1%5Bg_3%20%3D%3D%20i%5D%2C%20idata.posterior_predictive.y.isel(obs_id%3Dg_3%20%3D%3D%20i)%2C%20ax%3D_ax%5B1%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20_ax%5B0%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Conditional%20mean')%0A%20%20%20%20%20%20%20%20_ax%5B1%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Posterior%20predictive%20distribution')%0A%20%20%20%20%20%20%20%20for%20i%2C%20_%20in%20enumerate(group_list)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B2%5D.scatter(az.extract(idata%2C%20var_names%3D'%CE%B21')%5Bi%2C%20%3A%5D%2C%20az.extract(idata%2C%20var_names%3D'%CE%B20')%5Bi%2C%20%3A%5D%2C%20color%3Df'C%7Bi%7D'%2C%20alpha%3D0.01%2C%20rasterized%3DTrue%2C%20zorder%3D2)%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.set(xlabel%3D'slope'%2C%20ylabel%3D'intercept'%2C%20title%3D'Parameter%20space')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axhline(y%3D0%2C%20c%3D'k')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axvline(x%3D0%2C%20c%3D'k')%0A%20%20%20%20%20%20%20%20return%20_ax%0A%20%20%20%20_plot(idata2_1)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20In%20contrast%20to%20Model%201%2C%20when%20we%20consider%20groups%20we%20can%20see%20that%20now%20the%20evidence%20points%20toward%20_negative_%20relationships%20between%20%24x%24%20and%20%24y%24.%20We%20can%20see%20that%20from%20the%20negative%20slopes%20in%20the%20left%20and%20middle%20panels%20of%20the%20plot%20above%2C%20as%20well%20as%20from%20the%20majority%20of%20the%20posterior%20samples%20for%20the%20slope%20parameter%20being%20negative%20in%20the%20left%20panel%20above.%0A%0A%20%20%20%20The%20plot%20below%20takes%20a%20closer%20look%20at%20the%20group%20level%20slope%20parameters.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20figsize%2C%20idata2_1)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20_ax%20%3D%20az.plot_forest(idata2_1.posterior%5B'%CE%B21'%5D%2C%20combined%3DTrue%2C%20figsize%3Dfigsize)%0A%20%20%20%20_ax%5B0%5D.set(title%3D'Model%202%20suggests%20negative%20slopes%20for%20each%20group'%2C%20xlabel%3D'%24%5C%5Cbeta_1%24'%2C%20ylabel%3D'Group')%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Model%203%3A%20Partial%20pooling%20model%20with%20confounder%20included%0A%0A%20%20%20%20Model%203%20assumes%20the%20same%20causal%20DAG%20as%20model%202%20(see%20above).%20However%2C%20we%20can%20go%20further%20and%20incorporate%20more%20knowledge%20about%20the%20structure%20of%20our%20data.%20Rather%20than%20treating%20each%20group%20as%20entirely%20independent%2C%20we%20can%20use%20our%20knowledge%20that%20these%20groups%20are%20drawn%20from%20a%20population-level%20distribution.%20We%20could%20formalise%20this%20as%20saying%20that%20the%20group-level%20slopes%20and%20intercepts%20are%20modeled%20as%20deflections%20from%20a%20population-level%20slope%20and%20intercept%2C%20respectively.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20And%20we%20could%20describe%20this%20model%20mathematically%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%20%20%5Cbeta_0%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%205)%20%5C%5C%0A%20%20%20%20%20%20%5Cbeta_1%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%205)%20%5C%5C%0A%20%20%20%20%20%20p_%7B0%5Csigma%7D%2C%20p_%7B1%5Csigma%7D%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%20%20%5Cvec%7Bu_0%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%20p_%7B0%5Csigma%7D)%20%5C%5C%0A%20%20%20%20%20%20%5Cvec%7Bu_1%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%20p_%7B1%5Csigma%7D)%20%5C%5C%0A%20%20%20%20%20%20%5Csigma%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%20%20%5Cmu_i%20%26%3D%20%5Coverset%7B%5Ctext%7Bintercept%7D%7D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cbig(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Cunderset%7B%5Ctext%7Bpop%7D%7D%7B%5Cbeta_0%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5C%3B%20%7B%2B%7D%20%5C%3B%20%5Cunderset%7B%5Ctext%7Bgroup%7D%7D%7B%5Cvec%7Bu_0%7D%7B%5Bg_i%5D%7D%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cbig)%0A%20%20%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20%20%20%5C%3B%20%7B%2B%7D%20%5C%3B%20%5Coverset%7B%5Ctext%7Bslope%7D%7D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cbig(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Cunderset%7B%5Ctext%7Bpop%7D%7D%7B%5Cbeta_1%20%5Ccdot%20x_i%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5C%3B%20%7B%2B%7D%20%5C%3B%20%5Cunderset%7B%5Ctext%7Bgroup%7D%7D%7B%5Cvec%7Bu_1%7D%7B%5Bg_i%5D%7D%20%5Ccdot%20x_i%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cbig)%0A%20%20%20%20%20%20%20%20%20%20%7D%20%5C%5C%20%5C%5C%0A%20%20%20%20%20%20y_i%20%26%5Csim%20%5Ctext%7BNormal%7D(%5Cmu_i%2C%20%5Csigma)%0A%20%20%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%24%24%0A%0A%0A%20%20%20%20where%0A%20%20%20%20*%20%24%5Cbeta_0%24%20and%20%24%5Cbeta_1%24%20are%20the%20population%20level%20intercepts%20and%20slopes%2C%20respectively.%0A%20%20%20%20*%20%24%5Cvec%7Bu_0%7D%24%20and%20%24%5Cvec%7Bu_1%7D%24%20are%20group%20level%20deflections%20from%20the%20population%20level%20parameters.%0A%20%20%20%20*%20%24p_%7B0%5Csigma%7D%2C%20p_%7B1%5Csigma%7D%24%20are%20the%20standard%20deviations%20of%20the%20intercept%20and%20slope%20deflections%20and%20can%20be%20thought%20of%20as%20'shrinkage%20parameters'.%0A%20%20%20%20%20%20*%20In%20the%20limt%20of%20%24p_%7B0%5Csigma%7D%2C%20p_%7B1%5Csigma%7D%20%5Crightarrow%20%5Cinfty%24%20we%20recover%20the%20unpooled%20model.%0A%20%20%20%20%20%20*%20In%20the%20limit%20of%20%24p_%7B0%5Csigma%7D%2C%20p_%7B1%5Csigma%7D%20%5Crightarrow%200%24%20we%20recover%20the%20pooled%20model.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20can%20also%20express%20this%20Model%203%20in%20Wilkinson%20notation%20as%20%601%20%2B%20x%20%2B%20(1%20%2B%20x%20%7C%20g)%60.%0A%0A%20%20%20%20*%20The%20%601%60%20captures%20the%20global%20intercept%2C%20%24%5Cbeta_0%24.%0A%20%20%20%20*%20The%20%60x%60%20captures%20the%20global%20slope%2C%20%24%5Cbeta_1%24.%0A%20%20%20%20*%20The%20%60(1%20%2B%20x%20%7C%20g)%60%20term%20captures%20group%20specific%20terms%20for%20the%20intercept%20and%20slope.%0A%20%20%20%20%20%20*%20%601%20%7C%20g%60%20captures%20the%20group%20specific%20intercept%20deflections%20%24%5Cvec%7Bu_0%7D%24%20parameters.%0A%20%20%20%20%20%20*%20%60x%20%7C%20g%60%20captures%20the%20group%20specific%20slope%20deflections%20%24%5Cvec%7Bu_1%7D%5Bg_i%5D%24%20parameters.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(coords%2C%20data%2C%20pm)%3A%0A%20%20%20%20with%20pm.Model(coords%3Dcoords)%20as%20model3%3A%0A%20%20%20%20%20%20%20%20_%CE%B20%20%3D%20pm.Normal('%CE%B20'%2C%200%2C%205)%0A%20%20%20%20%20%20%20%20_%CE%B21%20%3D%20pm.Normal('%CE%B21'%2C%200%2C%205)%0A%20%20%20%20%20%20%20%20intercept_sigma%20%3D%20pm.Gamma('intercept_sigma'%2C%202%2C%202)%0A%20%20%20%20%20%20%20%20slope_sigma%20%3D%20pm.Gamma('slope_sigma'%2C%202%2C%202)%0A%20%20%20%20%20%20%20%20u0%20%3D%20pm.Normal('u0'%2C%200%2C%20intercept_sigma%2C%20dims%3D'group')%0A%20%20%20%20%20%20%20%20u1%20%3D%20pm.Normal('u1'%2C%200%2C%20slope_sigma%2C%20dims%3D'group')%0A%20%20%20%20%20%20%20%20_sigma%20%3D%20pm.Gamma('sigma'%2C%202%2C%202)%0A%20%20%20%20%20%20%20%20_x%20%3D%20pm.Data('x'%2C%20data.x%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20g_4%20%3D%20pm.Data('g'%2C%20data.group_idx%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20_%CE%BC%20%3D%20pm.Deterministic('%CE%BC'%2C%20_%CE%B20%20%2B%20u0%5Bg_4%5D%20%2B%20(_%CE%B21%20*%20_x%20%2B%20u1%5Bg_4%5D%20*%20_x)%2C%20dims%3D'obs_id')%0A%20%20%20%20%20%20%20%20pm.Normal('y'%2C%20mu%3D_%CE%BC%2C%20sigma%3D_sigma%2C%20observed%3Ddata.y%2C%20dims%3D'obs_id')%0A%20%20%20%20return%20(model3%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20DAG%20of%20this%20model%20highlights%20the%20scalar%20population%20level%20parameters%20%24%5Cbeta_0%24%20and%20%24%5Cbeta_1%24%20and%20the%20group%20level%20parameters%20%24%5Cvec%7Bu_0%7D%24%20and%20%24%5Cvec%7Bu_1%7D%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dot%2C%20model3%2C%20pm)%3A%0A%20%20%20%20dot(pm.model_to_graphviz(model3))%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20For%20the%20sake%20of%20completeness%2C%20there%20is%20another%20equivalent%20way%20to%20write%20this%20model.%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20p_%7B0%5Cmu%7D%2C%20p_%7B1%5Cmu%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(0%2C%205)%20%5C%5C%0A%20%20%20%20p_%7B0%5Csigma%7D%2C%20p_%7B1%5Csigma%7D%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%5Cvec%7B%5Cbeta_0%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(p_%7B0%5Cmu%7D%2C%20p_%7B0%5Csigma%7D)%20%5C%5C%0A%20%20%20%20%5Cvec%7B%5Cbeta_1%7D%20%26%5Csim%20%5Ctext%7BNormal%7D(p_%7B1%5Cmu%7D%2C%20p_%7B1%5Csigma%7D)%20%5C%5C%0A%20%20%20%20%5Csigma%20%26%5Csim%20%5Ctext%7BGamma%7D(2%2C%202)%20%5C%5C%0A%20%20%20%20%5Cmu_i%20%26%3D%20%5Cvec%7B%5Cbeta_0%7D%5Bg_i%5D%20%2B%20%20%5Cvec%7B%5Cbeta_1%7D%5Bg_i%5D%20%5Ccdot%20x_i%20%5C%5C%0A%20%20%20%20y_i%20%26%5Csim%20%5Ctext%7BNormal%7D(%5Cmu_i%2C%20%5Csigma)%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5Cvec%7B%5Cbeta_0%7D%24%20and%20%24%5Cvec%7B%5Cbeta_1%7D%24%20are%20the%20group-level%20parameters.%20These%20group%20level%20parameters%20can%20be%20thought%20of%20as%20being%20sampled%20from%20population%20level%20intercept%20distribution%20%24%5Ctext%7BNormal%7D(p_%7B0%5Cmu%7D%2C%20p_%7B0%5Csigma%7D)%24%20and%20population%20level%20slope%20distribution%20%24%5Ctext%7BNormal%7D(p_%7B1%5Cmu%7D%2C%20p_%7B1%5Csigma%7D)%24.%20So%20these%20distributions%20would%20represent%20what%20we%20might%20expect%20to%20observe%20for%20some%20as%20yet%20unobserved%20group.%0A%0A%20%20%20%20However%2C%20this%20formulation%20of%20the%20model%20does%20not%20so%20neatly%20map%20on%20to%20the%20Wilkinson%20notation.%20For%20this%20reason%2C%20we%20have%20chosen%20to%20present%20the%20model%20in%20the%20form%20given%20above.%20For%20an%20interesting%20discussion%20on%20this%20topic%2C%20see%20%5BDiscussion%20%23808%5D(https%3A%2F%2Fgithub.com%2Fbambinos%2Fbambi%2Fdiscussions%2F808)%20in%20the%20%5B%60bambi%60%5D(https%3A%2F%2Fgithub.com%2Fbambinos%2Fbambi)%20repository.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20hierarchical%20model%20we%20are%20considering%20contains%20a%20simplification%20in%20that%20the%20population%20level%20slope%20and%20intercept%20are%20assumed%20to%20be%20independent.%20It%20is%20possible%20to%20relax%20this%20assumption%20and%20model%20any%20correlation%20between%20these%20parameters%20by%20using%20a%20multivariate%20normal%20distribution.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20In%20one%20sense%20this%20move%20from%20Model%202%20to%20Model%203%20can%20be%20seen%20as%20adding%20parameters%2C%20and%20therefore%20increasing%20model%20complexity.%20However%2C%20in%20another%20sense%2C%20adding%20this%20knowledge%20about%20the%20nested%20structure%20of%20the%20data%20actually%20provides%20a%20constraint%20over%20parameter%20space.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Conduct%20inference%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(model3%2C%20pm%2C%20rng)%3A%0A%20%20%20%20with%20model3%3A%0A%20%20%20%20%20%20%20%20idata3%20%3D%20pm.sample(target_accept%3D0.95%2C%20random_seed%3Drng%2C%20cores%3D1%2C%20compile_kwargs%3D%7B%22mode%22%3A%20%22WASM%22%7D)%0A%20%20%20%20return%20(idata3%2C)%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20idata3)%3A%0A%20%20%20%20az.plot_trace(idata3%2C%20var_names%3D%5B%22~%CE%BC%22%5D)%3B%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Visualise%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20idata3%2C%20model3%2C%20np%2C%20predict)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-output%0A%20%20%20%20_n_points%20%3D%2010%0A%20%20%20%20_n_groups%20%3D%20len(data.group.unique())%0A%20%20%20%20xi_2%20%3D%20np.concatenate(%5Bnp.linspace(group%5B1%5D.x.min()%2C%20group%5B1%5D.x.max()%2C%20_n_points)%20for%20group%20in%20data.groupby('group_idx')%5D)%0A%20%20%20%20g_5%20%3D%20np.concatenate(%5B%5Bi%5D%20*%20_n_points%20for%20i%20in%20range(_n_groups)%5D).astype(int)%0A%20%20%20%20predict_at_1%20%3D%20%7B'x'%3A%20xi_2%2C%20'g'%3A%20g_5%7D%0A%20%20%20%20idata3_1%20%3D%20predict(model%3Dmodel3%2C%20idata%3Didata3%2C%20predict_at%3Dpredict_at_1)%0A%20%20%20%20return%20g_5%2C%20idata3_1%2C%20xi_2%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20data%2C%20g_5%2C%20group_list%2C%20idata3_1%2C%20plot_band%2C%20plt%2C%20xi_2)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20def%20_plot(idata)%3A%0A%20%20%20%20%20%20%20%20fig%2C%20_ax%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12%2C%204))%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(len(group_list))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B0%5D.scatter(data.x%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20data.y%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20plot_band(xi_2%5Bg_5%20%3D%3D%20i%5D%2C%20idata.posterior_predictive.%CE%BC.isel(obs_id%3Dg_5%20%3D%3D%20i)%2C%20ax%3D_ax%5B0%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B1%5D.scatter(data.x%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20data.y%5Bdata.group_idx%20%3D%3D%20i%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20plot_band(xi_2%5Bg_5%20%3D%3D%20i%5D%2C%20idata.posterior_predictive.y.isel(obs_id%3Dg_5%20%3D%3D%20i)%2C%20ax%3D_ax%5B1%5D%2C%20color%3Df'C%7Bi%7D')%0A%20%20%20%20%20%20%20%20_ax%5B0%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Conditional%20mean')%0A%20%20%20%20%20%20%20%20_ax%5B1%5D.set(xlabel%3D'x'%2C%20ylabel%3D'y'%2C%20title%3D'Posterior%20predictive%20distribution')%0A%20%20%20%20%20%20%20%20for%20i%2C%20_%20in%20enumerate(group_list)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%5B2%5D.scatter(az.extract(idata%2C%20var_names%3D'%CE%B21')%20%2B%20az.extract(idata%2C%20var_names%3D'u1')%5Bi%2C%20%3A%5D%2C%20az.extract(idata%2C%20var_names%3D'%CE%B20')%20%2B%20az.extract(idata%2C%20var_names%3D'u0')%5Bi%2C%20%3A%5D%2C%20color%3Df'C%7Bi%7D'%2C%20alpha%3D0.01%2C%20rasterized%3DTrue%2C%20zorder%3D2)%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.set(xlabel%3D'slope'%2C%20ylabel%3D'intercept'%2C%20title%3D'Parameter%20space')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axhline(y%3D0%2C%20c%3D'k')%0A%20%20%20%20%20%20%20%20_ax%5B2%5D.axvline(x%3D0%2C%20c%3D'k')%0A%20%20%20%20%20%20%20%20return%20_ax%0A%20%20%20%20_ax%20%3D%20_plot(idata3_1)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20panel%20on%20the%20right%20shows%20the%20group%20level%20posterior%20of%20the%20slope%20and%20intercept%20parameters%20as%20a%20contour%20plot.%20We%20can%20also%20just%20plot%20the%20marginal%20distribution%20below%20to%20see%20how%20much%20belief%20we%20have%20in%20the%20slope%20being%20less%20than%20zero.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(az%2C%20figsize%2C%20idata2_1)%3A%0A%20%20%20%20%23%20Cell%20tags%3A%20hide-input%0A%20%20%20%20_ax%20%3D%20az.plot_forest(idata2_1.posterior%5B'%CE%B21'%5D%2C%20combined%3DTrue%2C%20figsize%3Dfigsize)%5B0%5D%0A%20%20%20%20_ax.set(title%3D'Model%203%20suggests%20negative%20slopes%20for%20each%20group'%2C%20xlabel%3D'%24%5C%5Cbeta_1%24'%2C%20ylabel%3D'Group')%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Summary%0A%20%20%20%20Using%20Simpson's%20paradox%2C%20we've%20walked%20through%203%20different%20models.%20The%20first%20is%20a%20simple%20linear%20regression%20which%20treats%20all%20the%20data%20as%20coming%20from%20one%20group.%20This%20amounts%20to%20a%20causal%20DAG%20asserting%20that%20%24x%24%20causally%20influences%20%24y%24%20and%20%24%5Ctext%7Bgroup%7D%24%20was%20ignored%20(i.e.%20assumed%20to%20be%20causally%20unrelated%20to%20%24x%24%20or%20%24y%24).%20We%20saw%20that%20this%20lead%20us%20to%20believe%20the%20regression%20slope%20was%20positive.%0A%0A%20%20%20%20While%20that%20is%20not%20necessarily%20wrong%2C%20it%20is%20paradoxical%20when%20we%20see%20that%20the%20regression%20slopes%20for%20the%20data%20_within_%20a%20group%20is%20negative.%0A%0A%20%20%20%20This%20paradox%20is%20resolved%20by%20updating%20our%20causal%20DAG%20to%20include%20the%20group%20variable.%20This%20is%20what%20we%20did%20in%20the%20second%20and%20third%20models.%20Model%202%20was%20an%20unpooled%20model%20where%20we%20essentially%20fit%20separate%20regressions%20for%20each%20group.%0A%0A%20%20%20%20Model%203%20assumed%20the%20same%20causal%20DAG%2C%20but%20adds%20the%20knowledge%20that%20each%20of%20these%20groups%20are%20sampled%20from%20an%20overall%20population.%20This%20added%20the%20ability%20to%20make%20inferences%20not%20only%20about%20the%20regression%20parameters%20at%20the%20group%20level%2C%20but%20also%20at%20the%20population%20level.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Authors%0A%20%20%20%20*%20Authored%20by%20%5BBenjamin%20T.%20Vincent%5D(https%3A%2F%2Fgithub.com%2Fdrbenvincent)%20in%20July%202021%0A%20%20%20%20*%20Updated%20by%20%5BBenjamin%20T.%20Vincent%5D(https%3A%2F%2Fgithub.com%2Fdrbenvincent)%20in%20April%202022%0A%20%20%20%20*%20Updated%20by%20%5BBenjamin%20T.%20Vincent%5D(https%3A%2F%2Fgithub.com%2Fdrbenvincent)%20in%20February%202023%20to%20run%20on%20PyMC%20v5%0A%20%20%20%20*%20Updated%20to%20use%20%60az.extract%60%20by%20%5BBenjamin%20T.%20Vincent%5D(https%3A%2F%2Fgithub.com%2Fdrbenvincent)%20in%20February%202023%20(%5Bpymc-examples%23522%5D(https%3A%2F%2Fgithub.com%2Fpymc-devs%2Fpymc-examples%2Fpull%2F522))%0A%20%20%20%20*%20Updated%20by%20%5BBenjamin%20T.%20Vincent%5D(https%3A%2F%2Fgithub.com%2Fdrbenvincent)%20in%20September%202024%20(%5Bpymc-examples%23697%5D(https%3A%2F%2Fgithub.com%2Fpymc-devs%2Fpymc-examples%2Fpull%2F697)%20and%20%5Bpymc-examples%23709%5D(https%3A%2F%2Fgithub.com%2Fpymc-devs%2Fpymc-examples%2Fpull%2F709))%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A