{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "a6613312-ce57-4bd1-b080-03036ee0ca62", "metadata": {}, "outputs": [], "source": [ "import pymc as pm\n", "import numpy as np\n", "import arviz as az\n", "\n", "%load_ext lab_black" ] }, { "cell_type": "markdown", "id": "44a739d2-4f9f-4554-a2aa-06d9ac902cd4", "metadata": {}, "source": [ "# 5. Meta-analysis via Hierarchical Models*\n", "\n", "Adapted from [Unit 7: rats_nocentering.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit7/rats_nocentering.odc).\n", "\n", "Data for x can be found [here](https://raw.githubusercontent.com/areding/6420-pymc/main/data/rats_nocenteringx.txt), and Y [here](https://raw.githubusercontent.com/areding/6420-pymc/main/data/rats_nocenteringY.txt)." ] }, { "cell_type": "code", "execution_count": 2, "id": "a25c7103-cf18-4532-bc1d-5070ce2489f8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((5,), (30, 5))" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = np.loadtxt(\"../data/rats_nocenteringy.txt\")\n", "x = np.loadtxt(\"../data/rats_nocenteringx.txt\")\n", "x.shape, y.shape" ] }, { "cell_type": "markdown", "id": "cf2b3eab-42bb-4873-a01c-b9323a6372aa", "metadata": {}, "source": [ "This example is taken from Gelfand et al (1990), and concerns 30 young rats whose weights were measured weekly for five weeks." ] }, { "cell_type": "code", "execution_count": 3, "id": "ff75becd-7aa5-41b8-84e4-55e5d050b9b6", "metadata": { "tags": [ "hide-output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [mu_alpha, tau_alpha, mu_beta, tau_beta, tau, alpha, beta]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [16000/16000 00:27<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 27 seconds.\n" ] } ], "source": [ "with pm.Model() as m:\n", " mu_alpha = pm.Normal(\"mu_alpha\", 0, tau=1e-6)\n", " tau_alpha = pm.Gamma(\"tau_alpha\", 0.001, 0.001)\n", " mu_beta = pm.Normal(\"mu_beta\", 0, tau=1e-6)\n", " tau_beta = pm.Gamma(\"tau_beta\", 0.001, 0.001)\n", "\n", " tau = pm.Gamma(\"tau\", 0.001, 0.001)\n", " sigma = pm.Deterministic(\"sigma\", 1 / tau**0.5)\n", "\n", " alpha = pm.Normal(\"alpha\", mu_alpha, tau=tau_alpha, shape=(30, 1))\n", " beta = pm.Normal(\"beta\", mu_beta, tau=tau_beta, shape=(30, 1))\n", "\n", " mu = alpha + beta * x\n", " pm.Normal(\"y\", mu, tau=tau, observed=y)\n", "\n", " trace = pm.sample(3000)" ] }, { "cell_type": "code", "execution_count": 4, "id": "4994c305-2970-4e5d-be8b-5fd39597bb74", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%mcse_meanmcse_sdess_bulkess_tailr_hat
mu_alpha106.5622.294102.169110.8870.0200.01413847.09366.01.0
mu_beta6.1860.1065.9826.3820.0010.00113038.09634.01.0
tau_alpha0.0100.0040.0040.0170.0000.0009855.08919.01.0
tau_beta4.3131.4761.7556.9120.0140.01012379.09436.01.0
tau0.0270.0040.0200.0350.0000.0007883.09493.01.0
sigma6.1420.4675.2817.0290.0050.0047883.09493.01.0
\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \\\n", "mu_alpha 106.562 2.294 102.169 110.887 0.020 0.014 13847.0 \n", "mu_beta 6.186 0.106 5.982 6.382 0.001 0.001 13038.0 \n", "tau_alpha 0.010 0.004 0.004 0.017 0.000 0.000 9855.0 \n", "tau_beta 4.313 1.476 1.755 6.912 0.014 0.010 12379.0 \n", "tau 0.027 0.004 0.020 0.035 0.000 0.000 7883.0 \n", "sigma 6.142 0.467 5.281 7.029 0.005 0.004 7883.0 \n", "\n", " ess_tail r_hat \n", "mu_alpha 9366.0 1.0 \n", "mu_beta 9634.0 1.0 \n", "tau_alpha 8919.0 1.0 \n", "tau_beta 9436.0 1.0 \n", "tau 9493.0 1.0 \n", "sigma 9493.0 1.0 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(trace, var_names=[\"~alpha\", \"~beta\"])" ] }, { "cell_type": "code", "execution_count": 5, "id": "d66f026f-6824-413d-b7ab-74eee42e2af3", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "cluster30 x 1\n", "\n", "30 x 1\n", "\n", "\n", "cluster30 x 5\n", "\n", "30 x 5\n", "\n", "\n", "\n", "sigma\n", "\n", "sigma\n", "~\n", "Deterministic\n", "\n", "\n", "\n", "mu_alpha\n", "\n", "mu_alpha\n", "~\n", "Normal\n", "\n", "\n", "\n", "alpha\n", "\n", "alpha\n", "~\n", "Normal\n", "\n", "\n", "\n", "mu_alpha->alpha\n", "\n", "\n", "\n", "\n", "\n", "tau_alpha\n", "\n", "tau_alpha\n", "~\n", "Gamma\n", "\n", "\n", "\n", "tau_alpha->alpha\n", "\n", "\n", "\n", "\n", "\n", "mu_beta\n", "\n", "mu_beta\n", "~\n", "Normal\n", "\n", "\n", "\n", "beta\n", "\n", "beta\n", "~\n", "Normal\n", "\n", "\n", "\n", "mu_beta->beta\n", "\n", "\n", "\n", "\n", "\n", "tau_beta\n", "\n", "tau_beta\n", "~\n", "Gamma\n", "\n", "\n", "\n", "tau_beta->beta\n", "\n", "\n", "\n", "\n", "\n", "tau\n", "\n", "tau\n", "~\n", "Gamma\n", "\n", "\n", "\n", "tau->sigma\n", "\n", "\n", "\n", "\n", "\n", "y\n", "\n", "y\n", "~\n", "Normal\n", "\n", "\n", "\n", "tau->y\n", "\n", "\n", "\n", "\n", "\n", "alpha->y\n", "\n", "\n", "\n", "\n", "\n", "beta->y\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pm.model_to_graphviz(model=m)" ] }, { "cell_type": "code", "execution_count": 6, "id": "e36efc90-d365-41cd-b7f9-f9a36c179132", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Wed Mar 22 2023\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.0\n", "IPython version : 8.9.0\n", "\n", "pytensor: 2.10.1\n", "\n", "arviz: 0.14.0\n", "numpy: 1.24.2\n", "pymc : 5.1.2\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv -p pytensor" ] }, { "cell_type": "code", "execution_count": null, "id": "0ad643e0-2e96-4894-93c5-6d8a16b75989", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.0" } }, "nbformat": 4, "nbformat_minor": 5 }