15_hierarchical_modelling.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Hierarchical modelling\n",
        "\n",
        "## Why are linear models not enough?\n",
        "\n",
        "The linear models which we considered earlier have a set of assumptions behind them which are not always realistic. Such assumptions may enable analytical derivations in some cases, but these simplifications might be too limiting when we need to model real-life data. \n",
        "\n",
        "Epidemiological data is often complex, and we need to have appropriate tools to model complexity.\n",
        "\n",
        "`````{admonition} Group Task\n",
        ":class: tip\n",
        "Discuss with your neighbours how exactly the assumption of indendepdent errors $\\epsilon_i$ in linear regression $y_i = X_i \\beta + \\epsilon_i, \\quad \\epsilon_i \\sim \\mathcal{N}(0, \\sigma^2)$ helps in the classical setting. What does it enable us to compute analytically?\n",
        "`````\n",
        "\n",
        "One of the simplifying assumptions is <font color='orange'>homoskedasticity</font>, i.e. the assumption of <font color='orange'>equal or similar variances in different groups.</font>\n",
        "\n",
        "Sometimes it is indeed appropriate, but very often it is not.\n",
        "\n",
        "Here is a wider list of assumptions made by linear models:\n",
        "- homoskedasticity: equal or similar variances in different groups,\n",
        "- no error in predictors $X$,\n",
        "- no missing data,\n",
        "- normally distributed errors,\n",
        "- observations $y_i$ are independent.\n",
        "\n",
        "What we would often like to do, however, is \n",
        "- to model variance,\n",
        "- to capture errors in variables,\n",
        "- to allow for missing data,\n",
        "- use generalised linear models (GLMs),\n",
        "- use spatial and/or temporal error structure imposing correlation structure.\n",
        "\n",
        "## Hierarchies in data and parameters\n",
        "\n",
        "Hierarchical structures are commonly found in both natural data and statistical models. These hierarchies can represent various levels of organization or grouping within the data, and incorporating them into Bayesian inference can provide more accurate and insightful results. Such approach to modelling allows to account for <font color='orange'>different sources of variation</font> in the data.\n",
        "\n",
        "The hierarchical structure often arises in problems such as multiparameter models, where parameters can be regarded or connected in some way, or models of complex phenomena with different levels of hierarchy in the data. \n",
        "\n",
        "```{margin}\n",
        "Hierarchical models are also often called <font color='orange'>multilevel models</font>.\n",
        "```\n",
        "Generally, a hierarchical structure can be written if the joint distribution of the parameters can be decomposed to a <font color='orange'>series of conditional distributions</font>. The term <font color='orange'>hierarchical models</font> refers to a general set of modelling principles rather than to a specific family of models.\n",
        "\n",
        "Priors of model parameters are the first level of hierarchy. Priors of model parameters can depend on other parameters, and new priors (<font color='orange'>hyperpriors</font>) are defined over these prior parameters. In this case, the hyperpriors would be the second level of hierarchy. Following this scheme, the structure can be extended to more levels of hierarchy. In principle, there is not a limited level of hierarchy.\n",
        "\n",
        "Let us illustrate in an equation the posterior distribution of a model parameter with conditional structure and two levels of hierarchy. The prior distribution of model parameter $\\theta$, $p(\\theta|a)$, depends on the parameter $a$; the hyperprior $p(a|b)$ depends on a fixed value $b$ (hyper-parameter). \n",
        "\n",
        "$$\n",
        "p(\\theta|y) \\propto p(y|\\theta)p(\\theta|a)p(a|b).\n",
        "$$\n",
        "\n",
        "The sampling order in this example follows the computational graph: \n",
        "\n",
        "$$\n",
        "\\begin{array}{c}\n",
        "b\\\\ \\downarrow\\\\ a\\\\ \\downarrow\\\\ \\theta\\\\ \\downarrow\\\\ y\n",
        "\\end{array}\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here is a more specific example of hierarchical structure:\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\text{likelihood: } \\quad &p(y_i | \\mu_i, \\sigma) &&= \\mathcal{N}(y_i | \\mu_i, \\sigma), \\\\\n",
        "\\text{prior: } \\quad &p(\\mu_i) &&= \\mathcal{N}(\\mu_i | 0, \\sigma_\\mu)\\\\\n",
        "\\text{prior: } \\quad &p(\\sigma| a) &&= \\mathcal{HalfCauchy}(\\sigma|a)\\\\\n",
        "\\text{hyperprior: } &p(\\sigma_\\mu) &&= \\mathcal{HalfCauchy}(\\sigma_\\mu|b)\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Due to their flexibility, hierarchical models allow us to <font color='orange'>model variability in the parameters</font> of the model, <font color='orange'>partition variability</font> more explicitly into multiple terms, <font color='orange'>\"borrow strength\"</font> across groups.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Levels of pooling\n",
        "\n",
        "Hierarchical models exist on the continuum of two extreme cases: complete pooling and no pooling. In-between the two extremes is partial pooling.\n",
        "\n",
        "Let's explore each of these approaches and provide Numpyro code examples."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Complete Pooling\n",
        "\n",
        "In the \"complete pooling\" approach, all data points are treated as if they belong to a single group or population, and the model estimates a single set of parameters for the entire dataset. \n",
        "\n",
        "A pooled model implies that the data are sampled from the same model. This ignores all variation among the units being sampled. All observations share common parameter $\\theta$:\n",
        "\n",
        "```{margin}\n",
        "Image credit:  Chris Fonnesbeck\n",
        "```\n",
        "![](assets/pooled_model.png)\n",
        "\n",
        "Making predictions from a complete pooling model is straightforward."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/opt/anaconda3/envs/aims/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
            "  from .autonotebook import tqdm as notebook_tqdm\n"
          ]
        }
      ],
      "source": [
        "import numpyro\n",
        "import numpyro.distributions as dist\n",
        "from numpyro.infer import MCMC, NUTS\n",
        "\n",
        "import jax\n",
        "from jax import random\n",
        "import jax.numpy as jnp\n",
        "\n",
        "rng_key = random.PRNGKey(678)\n",
        "\n",
        "import pandas as pd\n",
        "\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [],
      "source": [
        "# artificial data for us to work with\n",
        "data = jnp.array([10, 12, 9, 11, 8]) # remember to turn data into a jnp array\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
            "        mu      9.92      0.77      9.94      8.64     11.13    460.42      1.00\n",
            "     sigma      1.65      0.53      1.55      0.93      2.37    388.83      1.00\n",
            "\n",
            "Number of divergences: 0\n"
          ]
        }
      ],
      "source": [
        "# Model\n",
        "def complete_pooling_model(data):\n",
        "    mu = numpyro.sample(\"mu\", dist.Normal(0, 10))\n",
        "    sigma = numpyro.sample(\"sigma\", dist.Exponential(1))\n",
        "    obs = numpyro.sample(\"obs\", dist.Normal(mu, sigma), obs=data)\n",
        "\n",
        "# Inference\n",
        "nuts_kernel = NUTS(complete_pooling_model)\n",
        "mcmc = MCMC(nuts_kernel, num_samples=1000, num_warmup=500, progress_bar=False)\n",
        "mcmc.run(rng_key, data)\n",
        "\n",
        "# Note how many mu-s and sigma-s are estimated\n",
        "mcmc.print_summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This approach assumes that there is no variation between data points, which can be overly restrictive when there is actual heterogeneity in the data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### No Pooling\n",
        "\n",
        "In the \"no pooling\" approach, each data point is treated independently without any grouping or hierarchical structure. This approach assumes that there is no shared information between data points, which can be overly simplistic when there is underlying structure or dependencies in the data.\n",
        "\n",
        "```{margin}\n",
        "Image credit:  Chris Fonnesbeck\n",
        "```\n",
        "![](assets/unpooled_model.png)\n",
        "\n",
        "Making predictions for new data from a no pooling model is impossible since it assumes no relashionship between the previously observed and new data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
            "      mu_0      9.98      1.06     10.00      8.17     11.35    338.50      1.00\n",
            "      mu_1     11.71      1.77     11.97      8.72     14.16    304.16      1.00\n",
            "      mu_2      8.85      1.35      8.92      6.53     10.57    419.85      1.00\n",
            "      mu_3     10.81      1.37     10.92      8.91     12.81    212.94      1.00\n",
            "      mu_4      7.87      1.23      7.96      6.34      9.82    115.95      1.00\n",
            "   sigma_0      0.94      0.82      0.67      0.08      2.02    259.38      1.00\n",
            "   sigma_1      1.15      1.11      0.81      0.06      2.53    117.35      1.00\n",
            "   sigma_2      1.10      0.99      0.85      0.07      2.29    242.58      1.01\n",
            "   sigma_3      1.06      0.89      0.84      0.08      2.20    102.34      1.00\n",
            "   sigma_4      1.01      0.91      0.80      0.05      2.12    147.39      1.00\n",
            "\n",
            "Number of divergences: 120\n"
          ]
        }
      ],
      "source": [
        "\n",
        "# Model\n",
        "def no_pooling_model(data):\n",
        "    for i, obs in enumerate(data):\n",
        "        mu_i = numpyro.sample(f\"mu_{i}\", dist.Normal(0, 10))\n",
        "        sigma_i = numpyro.sample(f\"sigma_{i}\", dist.Exponential(1))\n",
        "        numpyro.sample(f\"obs_{i}\", dist.Normal(mu_i, sigma_i), obs=data[i])\n",
        "\n",
        "# Inference\n",
        "nuts_kernel = NUTS(no_pooling_model)\n",
        "mcmc = MCMC(nuts_kernel, num_samples=1000, num_warmup=500, progress_bar=False)\n",
        "mcmc.run(rng_key, data)\n",
        "\n",
        "# Note how many mu-s and sigma-s are estimated\n",
        "mcmc.print_summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Partial Pooling\n",
        "\n",
        "In the \"partial pooling\" approach, the data is grouped into distinct categories or levels, and each group has its own set of parameters. However, these parameters are constrained by a shared distribution, allowing for both individual variation within groups and shared information across groups.\n",
        "\n",
        "```{margin}\n",
        "Image credit:  Chris Fonnesbeck\n",
        "```\n",
        "![](assets/partial_pooled_model.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Heart attack meta-analysis\n",
        "\n",
        "Consider the data from a meta-analysis of heart attack data (it is from Draper, \"Combining Information: Statistical Issues and Opportunities for Research\", 1992).\n",
        "\n",
        "Each data point represents the outcome of a study on post-heart attack survivorship. Each study involved administering aspirin to some victims immediately after the heart attack, while others did not receive aspirin. The `y` values denote the differences in mean survivorship observed in each study. Moreover, each study provided a standard deviation, calculated based on the relative sizes of the two groups within the study."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Data\n",
        "y = jnp.array([2.77, 2.50, 1.84, 2.56, 2.31, -1.15])\n",
        "sd = jnp.array([1.65, 1.31, 2.34, 1.67, 1.98, 0.90])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can build a model\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "y_i &\\sim \\mathcal{N}(\\theta_i, s_i),\\\\\n",
        "\\theta_i &\\sim \\mathcal{N}(\\mu, \\tau),\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "where $y_i$ is the mean of each study, and $s_i$ is the standard deviation for each study. Parameters $\\mu$ and $\\tau$ themselves can have priors\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\mu &\\sim \\mathcal{N}(\\bar{y}, 10s_y),\\\\\n",
        "\\tau &\\sim \\mathcal{HalfCauchy}(5s_y).\n",
        "\\end{align*}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "                mean       std    median      5.5%     94.5%     n_eff     r_hat\n",
            "        mu      1.54      1.14      1.48     -0.24      3.27    555.96      1.00\n",
            "       tau      1.86      1.18      1.64      0.23      3.16    168.80      1.00\n",
            "  theta[0]      2.17      1.29      2.10      0.19      4.27    642.40      1.00\n",
            "  theta[1]      2.08      1.13      1.94      0.28      3.77    567.66      1.00\n",
            "  theta[2]      1.66      1.55      1.52     -0.91      4.07    610.41      1.00\n",
            "  theta[3]      1.98      1.32      1.92     -0.18      3.96    562.49      1.00\n",
            "  theta[4]      1.83      1.45      1.77     -0.65      4.06    885.60      1.00\n",
            "  theta[5]     -0.37      0.96     -0.38     -1.98      1.02    193.60      1.00\n",
            "\n"
          ]
        }
      ],
      "source": [
        "def model(y, sd):\n",
        "    # hyperpriors\n",
        "    mu = numpyro.sample('mu', dist.Normal(jnp.mean(y), 10*jnp.std(y)))\n",
        "    tau = numpyro.sample('tau', dist.HalfCauchy(scale=5*jnp.std(y)))\n",
        "    \n",
        "    # prior\n",
        "    theta = numpyro.sample('theta', dist.Normal(mu, tau), sample_shape=(len(y),))\n",
        "\n",
        "    # Likelihood\n",
        "    numpyro.sample('obs', dist.Normal(theta, sd), obs=y)\n",
        "\n",
        "# Inference\n",
        "nuts_kernel = NUTS(model)\n",
        "mcmc = MCMC(nuts_kernel, num_samples=1000, num_warmup=500, progress_bar=False)\n",
        "mcmc.run(jax.random.PRNGKey(0), y=y, sd=sd)\n",
        "\n",
        "# Get posterior samples\n",
        "posterior_samples = mcmc.get_samples()\n",
        "\n",
        "# Print summary statistics\n",
        "numpyro.diagnostics.print_summary(posterior_samples, prob=0.89, group_by_chain=False)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "On the level of data, hierarchies usually are expressed in terms of groups:\n",
        "\n",
        "![](assets/pop_group.png)\n",
        "\n",
        "\n",
        "For example, $y_{2j}$ here could represent repeated measurements within the same group."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Sleep study\n",
        "\n",
        "```{margin}\n",
        "Belenky, Gregory, Nancy J. Wesensten, David R. Thorne, Maria L. Thomas, Helen C. Sing, Daniel P. Redmond, Michael B. Russo, and Thomas J. Balkin. 2003. “Patterns of Performance Degradation and Restoration During Sleep Restriction and Subsequent Recovery: A Sleep Dose-Response Study.” Journal of Sleep Research 12: 1–12.\n",
        "```\n",
        "The aim of the sleep study was to assess the imapct of sleep deprivation on reaction time. The `sleepstudy` dataset consists of records of 18 subjects. All subjects got a regular night’s sleep on “day 0” of the study, and were then restricted to 3 hours of sleep per night for the next 9 days. Each day, researchers recorded the subjects’ reaction times (in ms) on a series of tests."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Reaction</th>\n",
              "      <th>Days</th>\n",
              "      <th>Subject</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>249.5600</td>\n",
              "      <td>0</td>\n",
              "      <td>308</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>258.7047</td>\n",
              "      <td>1</td>\n",
              "      <td>308</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>250.8006</td>\n",
              "      <td>2</td>\n",
              "      <td>308</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>321.4398</td>\n",
              "      <td>3</td>\n",
              "      <td>308</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>356.8519</td>\n",
              "      <td>4</td>\n",
              "      <td>308</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Reaction  Days  Subject\n",
              "0  249.5600     0      308\n",
              "1  258.7047     1      308\n",
              "2  250.8006     2      308\n",
              "3  321.4398     3      308\n",
              "4  356.8519     4      308"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# To get the file using wget\n",
        "# !wget -O sleepstudy.csv https://raw.githubusercontent.com/elizavetasemenova/prob-epi/main/data/sleepstudy.csv\n",
        "# df = pd.read_csv('sleepstudy.csv')\n",
        "# Drop the colum\n",
        "# df = df.drop('Unnamed: 0', axis=1)\n",
        "# df.head()\n",
        "\n",
        "\n",
        "df = pd.read_csv('data/sleepstudy.csv')\n",
        "# Drop the column\n",
        "df = df.drop('Unnamed: 0', axis=1)\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`````{admonition} Task 22\n",
        ":class: tip\n",
        "\n",
        "- Draw a diagram showing the hierarchy in the data. What are the “groups” here?\n",
        "\n",
        "Complete pooling: \n",
        "- Suppose that we took a complete pooling approach to modelling `Reaction` time ($y$) by `Days` of sleep deprivation ($x$). Draw a diagram of complete pooling.\n",
        "- Ignore the subjects and fit the complete pooling model. Visualise the result.\n",
        "\n",
        "No pooling: \n",
        "- Construct and discuss separate scatterplots of Reaction by Days for each Subject.\n",
        "- Fit the no pooling model to the data.\n",
        "\n",
        "\n",
        "Partial pooling: \n",
        "- Plot only data for subjects `308` and `335`.\n",
        "- For these two subjects, provide a loose sketch of three separate trend lines corresponding to a completely pooled, not pooled, and partially pooled model \n",
        "- Fit a partially poolled model and visualise the results.\n",
        "\n",
        "`````"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`````{admonition} Task 23\n",
        ":class: tip\n",
        "\n",
        "```{margin}\n",
        "Cushny, Arthur R., and A. Roy Peebles. \"The action of optical isomers: II. Hyoscines.\" The Journal of Physiology (1905)\n",
        "```\n",
        "The dataset `sleep.csv` contains results of a small clinical trial.\n",
        "\n",
        "In this study, ten individuals were administered a sleep aid drug, followed by another sleep aid drug. The study tracked the additional hours of sleep each participant experienced under each drug regimen, relative to their baseline measurement.\n",
        "\n",
        "The dataset is relatively small, with a low sample size and only one predictor variable (the treatment). Luckily, the design of the experiment has a hierarchical structure.\n",
        "\n",
        "The data for this task can be loaded as follows:\n",
        "\n",
        "```\n",
        "df = pd.read_csv('data/sleep.csv')\n",
        "df.head()\n",
        "```\n",
        "\n",
        "- Plot the data by showing treatment ID on the x-axis, and the number of extra hours of sleep on the y-axis. Add a separate line for each individual.\n",
        "- Before fitting the model, what do you think the result will be: did drug 1 have an effect? did drug 2 have an effect?\n",
        "- Construct a hierarchical model to analyse this data (you might want to consult [this](https://rpubs.com/peopletrees/stan-sleep) analysis for an example).\n",
        "- What conclusions can you make from the analysis about drug 1 and drug 2?\n",
        "\n",
        "`````"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Random effects\n",
        "\n",
        "A common special case of hierarchical models are random effects.\n",
        "\n",
        "Consider the following model with a hierarchical mean and common variance modelling outcome $y_k$ at locations $k$:\n",
        "\n",
        "```{margin}\n",
        "We use $\\mathcal{IG}$ here as a short version of $\\mathcal{InverseGamma}$.\n",
        "```\n",
        "$$\n",
        "\\begin{align*}\n",
        "y_k \\sim \\mathcal{N}(\\mu_k, \\sigma^2)\\\\\n",
        "\\mu_k \\sim \\mathcal{N}(\\mu, \\tau^2)\\\\\n",
        "\\sigma^2 \\sim \\mathcal{IG}(s_1, s_2)\\\\\n",
        "\\mu \\sim \\mathcal{N}(\\mu_0, v_0)\\\\\n",
        "\\tau^2 \\sim \\mathcal{IG}(t_1, t_2)\\\\\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Often it is more convenient to write models in a slightly different form:\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "y_k \\sim \\mathcal{N}(\\mu_g + \\alpha_k, \\sigma^2)\\\\\n",
        "\\alpha_k \\sim \\mathcal{N}(0, \\tau^2)\\\\\n",
        "\\mu \\sim \\mathcal{N}(\\mu_0, v_0)\\\\\n",
        "\\tau^2 \\sim \\mathcal{IG}(t_1, t_2)\\\\\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Here $\\mu_g$ is the <font color='orange'>global mean</font> across all locations, and $\\alpha_k$ is a <font color='orange'>random effect</font> which shows how each location $k$ differs from the global mean.\n",
        "\n",
        "Due to this way of construction, random effects always have mean 0. \n",
        "\n",
        "Random effect variance attributes a portion of uncertainty to a specific source. \n",
        "\n",
        "The term <font color='orange'>random effects</font> is used to contrast the difference with <font color='orange'>fixed effects</font> in the linear predictor:\n",
        "\n",
        "$$\n",
        "\\eta  = \\beta X + \\epsilon\n",
        "$$\n",
        "\n",
        "Random effects are parts of the predictor in a model which would be different, if the experiment was replicated. And fixed effects are the parts of the model which would remain unchanged.\n",
        "\n",
        "Random effects often include some correlation structures, such as temporal random effect and spatial random effects. They attempt to account for the unexplained variance associated with a certain group due to everything that was not measured. \n",
        "\n",
        "Models which include both random and fixed effects are called <font color='orange'>mixed effects models</font>, and GLMs corresponding to them become GLMMs (<font color='orange'>generalised linear mixed models</font>)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Outro: comparison of traditional ML and Hierarchical Bayesian modelling\n",
        "\n",
        "|        | Traditional ML  | Hierarchical Bayesian modelling  |\n",
        "|---     |---|---|\n",
        "| Scale   | large  | small  |\n",
        "| Knowledge & structure | discard  | leverage  |\n",
        "| Data types | homogeneous   | heterogeneous   |\n",
        "| Approach | engineering-heavy, modelling light   | modelling-heavy   |"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.9.18"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}