leroux_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running on PyMC v4.1.3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/conor/miniconda3/envs/pymc_non_dev/lib/python3.10/site-packages/pymc/sampling_jax.py:36: UserWarning: This module is experimental.\n",
      "  warnings.warn(\"This module is experimental.\")\n"
     ]
    }
   ],
   "source": [
    "import aesara \n",
    "import aesara.tensor as at\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np \n",
    "import pymc as pm \n",
    "import xarray as xr \n",
    "import pandas as pd \n",
    "import pymc.sampling_jax\n",
    "import numpyro \n",
    "\n",
    "print(f\"Running on PyMC v{pm.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data: loading in the Scottish lip cancer data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scot_cancer = pd.read_csv(pm.get_data(\"scotland_lips_cancer.csv\"))\n",
    "\n",
    "# observed cancer counts\n",
    "y = df_scot_cancer[\"CANCER\"].values\n",
    "\n",
    "# number of observations\n",
    "N = len(y)\n",
    "\n",
    "# expected cancer counts E for each county: this is calculated using age-standardized rates of the local population\n",
    "E = df_scot_cancer[\"CEXP\"].values\n",
    "logE = np.log(E)\n",
    "\n",
    "# proportion of the population engaged in agriculture, forestry, or fishing\n",
    "x = df_scot_cancer[\"AFF\"].values / 10.0\n",
    "\n",
    "# spatial adjacency information: column `ADJ` contains list entries which are preprocessed to obtain adj as a list of lists\n",
    "adj = (\n",
    "    df_scot_cancer[\"ADJ\"].apply(lambda x: [int(val) for val in x.strip(\"][\").split(\",\")]).to_list()\n",
    ")\n",
    "\n",
    "# change to Python indexing (i.e. -1)\n",
    "for i in range(len(adj)):\n",
    "    for j in range(len(adj[i])):\n",
    "        adj[i][j] = adj[i][j] - 1\n",
    "\n",
    "# storing the adjacency matrix as a two-dimensional np.array\n",
    "adj_matrix = np.zeros((N, N), dtype=\"int32\")\n",
    "\n",
    "for area in range(N):\n",
    "    adj_matrix[area, adj[area]] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getting the necessary information to use for a spatial model specification \n",
    "N_edges = (adj_matrix == 1).sum()\n",
    "node1 = np.where(adj_matrix == 1)[0] \n",
    "node2 = np.where(adj_matrix == 1)[1] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = {\"num_areas\": N, \n",
    "          \"num_neighbours\": N_edges}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Leroux model\n",
    "\n",
    "We now use a Leroux model such that the random effect structure for each area $i$ is $\\psi_i=\\phi_i$, where the vector $\\boldsymbol{\\phi}$ is distributed according to \n",
    "\\begin{align}\n",
    "\\boldsymbol{\\phi}\\sim\\text{Normal}(0, [\\tau^2\\rho(\\mathbf{D}-\\mathbf{W}) + (1-\\rho)\\mathbf{I})\\big ]^{-1}). \n",
    "\\end{align} \n",
    "When $\\rho=1$, the Leroux prior is equal to an ICAR prior, and when $\\rho=0$, it is equal to an independent random effect. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating the pairwise specification for the Leroux model\n",
    "def pairwise_diff_leroux(rho, phi, node1, node2):\n",
    "    return -0.5 * rho * ((phi[node1]-phi[node2]) ** 2).sum()\n",
    "# creating the square sum potential specification \n",
    "def square_sum(rho, phi):\n",
    "    return -0.5 * (1-rho) * (phi ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model(coords={\"num_areas\": np.arange(N)}) as leroux_model:\n",
    "    # precision priors, transform to standard deviation \n",
    "    tau_phi = pm.Gamma(\"tau_phi\", alpha=1, beta=1)\n",
    "    # transform to standard deviation \n",
    "    sigma_phi = pm.Deterministic(\"sigma_phi\", 1/at.sqrt(tau_phi))\n",
    "    # spatial smoothing effect prior \n",
    "    rho = pm.Beta(\"rho\",  alpha=1, beta=1)\n",
    "\n",
    "    # spatial random effect\n",
    "    phi = pm.Flat(\"phi\", dims=\"num_areas\")\n",
    "    pm.Potential(\"spatial_diff\", pairwise_diff_leroux(rho, phi, node1, node2))\n",
    "    pm.Potential(\"square_sum\", square_sum(rho, phi))\n",
    "\n",
    "    # constraint on the spatial random effect \n",
    "    zero_constraint = pm.Normal.dist(mu=0.0, sigma=np.sqrt(0.001))\n",
    "    pm.Potential(\"zero_sum_phi\", pm.logp(zero_constraint, pm.math.sum(phi))) \n",
    "\n",
    "    # regression coefficient priors \n",
    "    beta0 = pm.Normal(\"beta0\", mu=0, sigma=5)\n",
    "    beta1 = pm.Normal(\"beta1\", mu=0, sigma=5)\n",
    "    \n",
    "    # linear predictor \n",
    "    eta = pm.Deterministic(\"eta\", logE + beta0 + beta1*x + sigma_phi*phi, dims=\"num_areas\") \n",
    "\n",
    "    # likelihood\n",
    "    obs = pm.Poisson(\"obs\", at.exp(eta), observed=y, dims=\"num_areas\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "INFO:pymc:Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "INFO:pymc:Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "INFO:pymc:Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [tau_phi, rho, phi, beta0, beta1]\n",
      "INFO:pymc:NUTS: [tau_phi, rho, phi, beta0, beta1]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 04:05<00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/conor/miniconda3/envs/pymc_non_dev/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/conor/miniconda3/envs/pymc_non_dev/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/conor/miniconda3/envs/pymc_non_dev/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/conor/miniconda3/envs/pymc_non_dev/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 260 seconds.\n",
      "INFO:pymc:Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 260 seconds.\n"
     ]
    }
   ],
   "source": [
    "with leroux_model:\n",
    "    idata = pm.sample()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "f977084116a5d06351f0969c43dbc528856b1012f23a32e53dbf92034586142c"
  },
  "kernelspec": {
   "display_name": "Python 3.10.5 ('pymc_non_dev')",
   "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.10.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}