econmt-statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn about how to do statistics with code. We already saw some statistics in the chapter on probability and random processes: here we'll focus on computing basic statistics and using statistical tests. We'll make use of the excellent [*pingouin*](https://pingouin-stats.org/index.html) statistics package and its documentation for many of the examples and methods in this chapter {cite:ps}`vallat2018pingouin`. This chapter also draws on Open Intro Statistics {cite:ps}`diez2012openintro`.\n",
    "\n",
    "### Notation and basic definitions\n",
    "\n",
    "Greek letters, like $\\beta$, are the truth and represent parameters. Modified Greek letters are an estimate of the truth, for example $\\hat{\\beta}$. Sometimes Greek letters will stand in for vectors of parameters. Most of the time, upper case Latin characters such as $X$ will represent random variables (which could have more than one dimension). Lower case letters from the Latin alphabet denote realised data, for instance $x$ (which again could be multi-dimensional).  Modified Latin alphabet letters denote computations performed on data, for instance $\\bar{x} = \\frac{1}{n} \\displaystyle\\sum_{i} x_i$ where $n$ is number of samples. Parameters are given following a vertical bar, for example if $f(x|\\mu, \\sigma)$ is a probability density function, the vertical line indicates that its parameters are $\\mu$ and $\\sigma$. The set of distributions with densities $f_\\theta(x)$, $\\theta \\in \\Theta$ is called a parametric family, eg there is a family of different distributions that are parametrised by $\\theta$.\n",
    "\n",
    "A **statistic** $T(x)$ is a function of the data $x=(x_1, \\dots, x_n)$. \n",
    "\n",
    "An **estimator** of a parameter $\\theta$ is a function $T=T(x)$ which is used to estimate $\\theta$ based on observations of data. $T$ is an unbiased estimator if $\\mathbb{E}(T) = \\theta$.\n",
    "\n",
    "If $X$ has PDF $f(x|\\theta)$ then, given the observed value $x$ of $X$, the **likelihood** of $\\theta$ is defined by $\\text{lik}(\\theta) = f(x | \\theta)$. For independent and identically distributed observed values, then $\\text{lik}(\\theta) = f(x_1, \\dots, x_n| \\theta) = \\Pi_{i=1}^n f(x_i | \\theta)$. The $\\hat{\\theta}$ such that this function attains its maximum value is the **maximum likelihood estimator (MLE)** of $\\theta$.\n",
    "\n",
    "Given an MLE $\\hat{\\theta}$ of $\\theta$, $\\hat{\\theta}$ is said to be **consistent** if $\\mathbb{P}(\\hat{\\theta} - \\theta > \\epsilon) \\rightarrow 0$ as $n\\rightarrow \\infty$.\n",
    "\n",
    "An estimator *W* is **efficient** relative to another estimator $V$ if $\\text{Var}(W) < \\text{Var}(V)$.\n",
    "\n",
    "Let $\\alpha$ be the 'significance level' of a test statistic $T$.\n",
    "\n",
    "Let $\\gamma(X)$ and $\\delta(X)$ be two statistics satisfying $\\gamma(X) < \\delta(X)$ for all $X$. If on observing $X = x$, the inference can be made that $\\gamma(x) \\leq \\theta \\leq \\delta(x)$. Then $[\\delta(x), \\gamma(x)]$ is an **interval estimate** and $[\\delta(X), \\gamma(X)]$ is an **interval estimator**. The random interval (random because the *endpoints* are random variables) $[\\delta(X), \\gamma(X)]$ is called a $100\\cdot\\alpha \\%$ **confidence interval** for $\\theta$. Of course, there is a true $\\theta$, so either it is in this interval or it is not. But if the confidence interval was constructed many times over using samples, $\\theta$ would be contained within it $100\\cdot\\alpha \\%$ of the times.\n",
    "\n",
    "A **hypothesis test** is a conjecture about the distribution of one or more random variables, and a test of a hypothesis is a procedure for deciding whether or not to reject that conjecture. The **null hypothesis**, $H_0$, is only ever conservatively rejected and represents the default position. The **alternative hypothesis**, $H_1$, is the conclusion contrary to this.\n",
    "\n",
    "A type I error occurs when $H_0$ is rejected when it is true, ie when a *true* null hypothesis is rejected. Mistakenly failing to reject a false null hypothesis is called a type II error.\n",
    "\n",
    "\n",
    "In the most simple situations, the upper bound on the probability of a type I error is called the size or **significance level** of the *test*. The **p-value** of a random variable $X$ is the smallest value of the significance level (denoted $\\alpha$) for which $H_0$ would be rejected on the basis of seeing $x$. The p-value is sometimes called the significance level of $X$. The probability that a test will reject the null is called the power of the test. The probability of a type II error is equal to 1 minus the power of the test.\n",
    "\n",
    "\n",
    "Recall that there are two types of statistics out there: parametrised, eg by $\\theta$, and non-parametrised. The latter are often distribution free (ie don't involve a PDF) or don't require parameters to be specified.\n",
    "\n",
    "### Imports\n",
    "\n",
    "First we need to import the packages we'll be using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pingouin as pg\n",
    "import statsmodels.formula.api as smf\n",
    "from scipy import stats\n",
    "\n",
    "# Set seed for random numbers\n",
    "seed_for_prng = 78557\n",
    "prng = np.random.default_rng(\n",
    "    seed_for_prng\n",
    ")  # prng=probabilistic random number generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic statistics\n",
    "\n",
    "Let's start with computing the simplest statistics you can think of using some synthetic data. Many of the functions have lots of extra options that we won't explore here (like weights or normalisation); remember that you can see these using the `help()` method. \n",
    "\n",
    "We'll generate a vector with 100 entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(range(100))\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "from myst_nb import glue\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "dict_fns = {\n",
    "    \"mean\": np.mean(data),\n",
    "    \"std\": np.std(data),\n",
    "    \"mode\": stats.mode([0, 1, 2, 3, 3, 3, 5])[0],\n",
    "    \"median\": np.median(data),\n",
    "}\n",
    "\n",
    "for name, eval_fn in dict_fns.items():\n",
    "    glue(name, f\"{eval_fn:.1f}\")\n",
    "\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)\n",
    "# Plot settings\n",
    "plt.style.use(\"plot_style.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, let's see how some basic statistics are computed. The mean is `np.mean(data)=` {glue:}`mean`, the standard deviation is `np.std(data)=` {glue:}`std`, and the median is given by `np.median(data)= `{glue:}`median`. The mode is given by `stats.mode([0, 1, 2, 3, 3, 3, 5])[0]=` {glue:}`mode` (access the counts using `stats.mode(...)[1]`).\n",
    "\n",
    "Less famous quantiles than the median are given by, for example for $q=0.25$,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.quantile(data, 0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with **pandas**, **numpy** and **scipy** work on scalars, vectors, matrices, and tensors: you just need to specify the axis that you'd like to apply a function to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.fromfunction(lambda i, j: i + j, (3, 6), dtype=int)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(data, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that, for discrete data points, the $k$th (unnormalised) moment is\n",
    "\n",
    "$$\n",
    "\\hat{m}_k = \\frac{1}{n}\\displaystyle\\sum_{i=1}^{n} \\left(x_i - \\bar{x}\\right)^k\n",
    "$$\n",
    "\n",
    "To compute this use scipy's `stats.moment(a, moment=1)`. For instance for the kurtosis ($k=4$), it's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.moment(data, moment=4, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Covariances are found using `np.cov`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cov(np.array([[0, 1, 2], [2, 1, 0]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, as expected, the $C_{01}$ term is -1 as the vectors are anti-correlated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parametric tests\n",
    "\n",
    "Reminder: parametric tests assume that data are effectively drawn a probability distribution that can be described with fixed parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-sample t-test\n",
    "\n",
    "The one-sample t-test tells us whether a given parameter for the mean, i.e. a suspected $\\mu$, is likely to be consistent with the sample mean. The null hypothesis is that $\\mu = \\bar{x}$. Let's see an example using the default `alternative='two-sided'` option. Imagine we have data on the number of hours people spend working each day and we want to test the (alternative) hypothesis that $\\bar{x}$ is not $\\mu=$8 hours:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [8.5, 5.4, 6.8, 9.6, 4.2, 7.2, 8.8, 8.1]\n",
    "\n",
    "pg.ttest(x, 8).round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The returned object is a **pandas** dataframe.) We only have 8 data points, and so that is a great big confidence interval! It's worth remembering what a t-statistic and t-test really are. In this case, the statistic that is constructed to test whether the sample mean is different from a known parameter $\\mu$ is\n",
    "\n",
    "$$\n",
    "T = \\frac{\\sqrt{n}(\\bar{x}-\\mu)}{\\hat{\\sigma}} \\thicksim t_{n-1}\n",
    "$$\n",
    "\n",
    "where $t_{n-1}$ is the student's t-distribution and $n-1$ is the number of degrees of freedom. The $100\\cdot(1-\\alpha)\\%$ test interval in this case is given by\n",
    "\n",
    "$$\n",
    "1 - \\alpha = \\mathbb{P}\\left(-t_{n-1, \\alpha/2} \\leq \\frac{\\sqrt{n}(\\bar{x} - \\mu)}{\\hat{\\sigma}} \\leq t_{n-1,\\alpha/2}\\right)\n",
    "$$\n",
    "\n",
    "where we define $t_{n-1, \\alpha/2}$ such that $\\mathbb{P}(T > t_{n-1, \\alpha/2}) = \\alpha/2$. For $\\alpha=0.05$, implying confidence intervals of 95%, this looks like:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "\n",
    "\n",
    "def plot_t_stat(x, mu):\n",
    "    T = np.linspace(-7, 7, 500)\n",
    "    pdf_vals = st.t.pdf(T, len(x) - 1)\n",
    "\n",
    "    sigma_hat = np.sqrt(np.sum((x - np.mean(x)) ** 2) / (len(x) - 1))\n",
    "    actual_T_stat = (np.sqrt(len(x)) * (np.mean(x) - mu)) / sigma_hat\n",
    "\n",
    "    alpha = 0.05\n",
    "    T_alpha_over_2 = st.t.ppf(1.0 - alpha / 2, len(x) - 1)\n",
    "\n",
    "    interval_T = T[((T > -T_alpha_over_2) & (T < T_alpha_over_2))]\n",
    "    interval_y = pdf_vals[((T > -T_alpha_over_2) & (T < T_alpha_over_2))]\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(T, pdf_vals, label=f\"Student t: dof={len(x)-1}\", zorder=2)\n",
    "    ax.fill_between(\n",
    "        interval_T, 0, interval_y, alpha=0.2, label=\"95% interval\", zorder=1\n",
    "    )\n",
    "    ax.plot(\n",
    "        actual_T_stat,\n",
    "        st.t.pdf(actual_T_stat, len(x) - 1),\n",
    "        \"bo\",\n",
    "        ms=15,\n",
    "        label=r\"$\\sqrt{n} (\\bar{x} - \\mu) / \\hat{\\sigma}$\",\n",
    "        color=\"orchid\",\n",
    "        zorder=4,\n",
    "    )\n",
    "    ax.vlines(\n",
    "        actual_T_stat, 0, st.t.pdf(actual_T_stat, len(x) - 1), color=\"orchid\", zorder=3\n",
    "    )\n",
    "    ax.set_xlabel(\"Value of statistic T\")\n",
    "    ax.set_ylabel(\"PDF\")\n",
    "    ax.set_xlim(-7, 7)\n",
    "    ax.set_ylim(0.0, 0.4)\n",
    "    ax.legend(frameon=False)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "mu = 8\n",
    "plot_t_stat(x, mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we would reject the alternative hypothesis. You can see why from the plot; the test statistic we have constructed lies within the interval where we cannot reject the null hypothesis. $\\bar{x}-\\mu$ is close enough to zero to give us cause for concern. (You can also see from the plot why this is a two-tailed test: we don't care if $\\bar{x}$ is greater or less than $\\mu$, just that it's different--and so the test statistic could appear in either tail of the distribution for us to accept $H_1$.)\n",
    "\n",
    "We accept the null here, but about if there were many more data points? Let's try adding some generated data (pretend it is from making extra observations).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 'Observe' extra data\n",
    "extra_data = prng.uniform(5.5, 8.5, size=(30))\n",
    "# Add it in to existing vector\n",
    "x_prime = np.concatenate((np.array(x), extra_data), axis=None)\n",
    "# Run t-test\n",
    "pg.ttest(x_prime, 8).round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, what happened? Our extra observations have seen the confidence interval shrink considerably, and the p-value is effectively 0. There's a large negative t-statistic too. Unsurprisingly, as we chose a uniform distribution that only just included 8 but was centered on $(8-4.5)/2$ *and* we had more points, the test now rejects the null hypothesis that $\\mu=8$ . Because the alternative hypothesis is just $\\mu\\neq8$, and these tests are conservative, we haven't got an estimate of what the mean actually is; we just know that our test rejects that it's $8$.\n",
    "\n",
    "We can see this in a new version of the chart that uses the extra data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_t_stat(x_prime, mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our test statistic is safely outside the interval.\n",
    "\n",
    "#### Connection to linear regression\n",
    "\n",
    "Note that testing if $\\mu\\neq0$ is equivalent to having the alternative hypothesis that a single, non-zero scalar value is a good expected value for $x$, i.e. that $\\mathbb{E}(x) \\neq 0$. Which may sound familiar if you've run **linear regression** and, indeed, this t-test has an equivalent linear model! It's just regressing $X$ on a constant--a single, non-zero scalar value. In general, t-tests appear in linear regression to test whether any coefficient $\\beta \\neq 0$. \n",
    "\n",
    "We can see this connection by running a hypothesis test of whether the sample mean is not zero. Note the confidence interval, t-statistic, and p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg.ttest(x, 0).round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, as an alternative, regressing x on a constant, again noting the interval, t-stat, and p-value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(x, columns=[\"x\"])\n",
    "\n",
    "res = smf.ols(formula=\"x ~ 1\", data=df).fit()\n",
    "# Show only the info relevant to the intercept (there are no other coefficients)\n",
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many tests have an equivalent  linear model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other information provided by **Pingouin** tests\n",
    "\n",
    "We've covered the degrees of freedom, the T statistic, the p-value, and the confidence interval. So what's all that other gunk in our t-test? Cohen's d is a measure of whether the difference being measured in our test is large or not (this is important; you can have statistically significant differences that are so small as to be inconsequential). Cohen suggested that $d = 0.2$ be considered a 'small' effect size, 0.5 represents a 'medium' effect size and 0.8 a 'large' effect size. BF10 represents the Bayes factor, the ratio (given the data) of the likelihood of the alternative hypothesis relative to the null hypothesis. Values greater than unity therefore favour the alternative hypothesis. Finally, power is the achieved power of the test, which is $1 - \\mathbb{P}(\\text{type II error})$. A common default to have in mind is a power greater than 0.8."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two-sample t-test\n",
    "\n",
    "The two-sample t-test is used to determine if two population means are equal (with the null being that they *are* equal). Let's look at an example with synthetic data of equal length, which means we can use the *paired* version of this. We'll imagine we are looking at an intervention with a pre- and post- dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre = [5.5, 2.4, 6.8, 9.6, 4.2, 5.9]\n",
    "post = [6.4, 3.4, 6.4, 11.0, 4.8, 6.2]\n",
    "pg.ttest(pre, post, paired=True, alternative=\"two-sided\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we cannot reject the null hypothesis that the means are the same pre- and post-intervention."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pearson correlation\n",
    "\n",
    "The Pearson correlation coefficient measures the linear relationship between two datasets. Strictly speaking, it requires that each dataset be normally distributed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, cov = [4, 6], [(1, 0.5), (0.5, 1)]\n",
    "x, y = prng.multivariate_normal(mean, cov, 30).T\n",
    "# Compute Pearson correlation\n",
    "pg.corr(x, y).round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Welch's t-test\n",
    "\n",
    "In the case where you have two samples with unequal variances (or, really, unequal sample sizes too), Welch's t-test is appropriate. With `correction='true'`, it assumes that variances are not equal.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = prng.normal(loc=7, size=20)\n",
    "y = prng.normal(loc=6.5, size=15)\n",
    "pg.ttest(x, y, correction=\"true\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-way ANOVA\n",
    "\n",
    "Analysis of variance (ANOVA) is a technique for testing hypotheses about means, for example testing the equality of the means of $k>2$ groups. The model would be\n",
    "\n",
    "$$\n",
    "X_{ij} = \\mu_i + \\epsilon_{ij} \\quad j=1, \\dots, n_i \\quad i=1, \\dots, k.\n",
    "$$\n",
    "\n",
    "so that the $i$th group has $n_i$ observations. The null hypothesis of one-way ANOVA is that $H_0: \\mu_1 = \\mu_2 = \\dots = \\mu_k$, with the alternative hypothesis that this is *not* true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pg.read_dataset(\"mixed_anova\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the ANOVA\n",
    "pg.anova(data=df, dv=\"Scores\", between=\"Group\", detailed=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple pairwise t-tests\n",
    "\n",
    "There's a problem with running multiple t-tests: if you run enough of them, something is bound to come up as significant! As such, some *post-hoc* adjustments exist that correct for the fact that multiple tests are occurring simultaneously. In the example below, multiple pairwise comparisons are made between the scores by time group. There is a corrected p-value, `p-corr`, computed using the Benjamini/Hochberg FDR correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg.pairwise_ttests(\n",
    "    data=df,\n",
    "    dv=\"Scores\",\n",
    "    within=\"Time\",\n",
    "    subject=\"Subject\",\n",
    "    parametric=True,\n",
    "    padjust=\"fdr_bh\",\n",
    "    effsize=\"hedges\",\n",
    ").round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-way ANCOVA\n",
    "\n",
    "Analysis of covariance (ANCOVA) is a general linear model which blends ANOVA and regression. ANCOVA evaluates whether the means of a dependent variable (dv) are equal across levels of a categorical independent variable (between) often called a treatment, while statistically controlling for the effects of other continuous variables that are not of primary interest, known as covariates or nuisance variables (covar)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pg.read_dataset(\"ancova\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg.ancova(data=df, dv=\"Scores\", covar=\"Income\", between=\"Method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power calculations\n",
    "\n",
    "Often, it's quite useful to know what sample size is needed to avoid certain types of testing errors. **Pingouin** offers ways to compute effect sizes and test powers to help with these questions.\n",
    "\n",
    "As an example, let's assume we have a new drug (`x`) and an old drug (`y`) that are both intended to reduce blood pressure. The standard deviation of the reduction in blood pressure of those receiving the old drug is 12 units. The null hypothesis is that the new drug is no more effective than the new drug. But it will only be worth switching production to the new drug if it reduces blood pressure by more than 3 units versus the old drug. In this case, the effect size of interest is 3 units.\n",
    "\n",
    "Let's assume for a moment that the true difference is 3 units and we want to perform a test with $\\alpha=0.05$. The problem is that, for small differences in the effect, the distribution of effects under the null and the distribution of effects under the alternative have a great deal of overlap. So the chances of making a Type II error - accepting the null hypothesis when it is actually false - is quite high. Let's say we'd ideally have at most a 20% chance of making a Type II error: what sample size do we need?\n",
    "\n",
    "We can compute this, but we need an extra piece of information first: a normalised version of the effect size called Cohen's $d$. We need to transform the difference in means to compute this. For independent samples, $d$ is:\n",
    "\n",
    "$$ d = \\frac{\\overline{X} - \\overline{Y}}{\\sqrt{\\frac{(n_{1} - 1)\\sigma_{1}^{2} + (n_{2} - 1)\\sigma_{2}^{2}}{n_1 + n_2 - 2}}}.$$\n",
    "\n",
    "(If you have real data samples, you can compute this using `pg.compute_effsize()`.)\n",
    "\n",
    "For this case, $d$ is $-3/12 = -1/4$ if we assume the standard deviations are the same across the old (`y`) and new (`x`) drugs. So we will plug that $d$ in and look at a range of possible sample sizes along with a standard value for $alpha$ of 0.05. In the below `alternative=less` tests the alternative that `x` has a smaller mean than `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cohen_d = -0.25  # Fixed effect size\n",
    "sample_size_array = np.arange(1, 500, 50)  # Incrementing sample size\n",
    "# Compute the achieved power\n",
    "pwr = pg.power_ttest(\n",
    "    d=cohen_d,\n",
    "    n=sample_size_array,\n",
    "    alpha=0.05,\n",
    "    contrast=\"two-samples\",\n",
    "    alternative=\"less\",\n",
    ")\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(sample_size_array, pwr, \"ko-.\")\n",
    "ax.axhline(0.8, color=\"r\", ls=\":\")\n",
    "ax.set_xlabel(\"Sample size\")\n",
    "ax.set_ylabel(\"Power (1 - type II error)\")\n",
    "ax.set_title(\"Achieved power of a T-test\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this, we can see we need a sample size of at least 200 in order to have a power of 0.8. \n",
    "\n",
    "The `pg.power_ttest()` function takes any three of the four of `d`, `n`, `power`, and `alpha` (ie leave one of these out), and then returns what the missing parameter should be. We passed in `d`, `n`, and `alpha`, and so the `power` was returned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-parametric tests\n",
    "\n",
    "Reminder: non-parametrics tests do not make any assumptions about the distribution from which data are drawn or that it can be described by fixed parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wilcoxon Signed-rank Test\n",
    "\n",
    "This tests the null hypothesis that two related paired samples come from the same distribution. It is the non-parametric equivalent of the t-test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [20, 22, 19, 20, 22, 18, 24, 20, 19, 24, 26, 13]\n",
    "y = [38, 37, 33, 29, 14, 12, 20, 22, 17, 25, 26, 16]\n",
    "pg.wilcoxon(x, y, alternative=\"two-sided\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mann-Whitney U Test (aka Wilcoxon rank-sum test)\n",
    "\n",
    "The Mann–Whitney U test is a non-parametric test of the null hypothesis that it is equally likely that a randomly selected value from one sample will be less than or greater than a randomly selected value from a second sample. It is the non-parametric version of the two-sample T-test.\n",
    "\n",
    "Like many non-parametric **pingouin** tests, it can take values of tail that are 'two-sided', 'one-sided', 'greater', or 'less'. Below, we ask if a randomly selected value from `x` is greater than one from `y`, with the null that it is not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = prng.uniform(low=0, high=1, size=20)\n",
    "y = prng.uniform(low=0.2, high=1.2, size=20)\n",
    "pg.mwu(x, y, alternative=\"greater\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spearman Correlation\n",
    "\n",
    "The Spearman correlation coefficient is the Pearson correlation coefficient between the rank variables, and does not assume normality of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, cov = [4, 6], [(1, 0.5), (0.5, 1)]\n",
    "x, y = prng.multivariate_normal(mean, cov, 30).T\n",
    "pg.corr(x, y, method=\"spearman\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kruskal-Wallace\n",
    "\n",
    "The Kruskal-Wallis H-test tests the null hypothesis that the population median of all of the groups are equal. It is a non-parametric version of ANOVA. The test works on 2 or more independent samples, which may have different sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pg.read_dataset(\"anova\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg.kruskal(data=df, dv=\"Pain threshold\", between=\"Hair color\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Chi-Squared Test\n",
    "\n",
    "The chi-squared test is used to determine whether there is a significant difference between the expected frequencies and the observed frequencies in one or more categories. This test can be used to evaluate the quality of a categorical variable in a classification problem or to check the similarity between two categorical variables.\n",
    "\n",
    "There are two conditions for a chi-squared test:\n",
    "\n",
    "- Independence: Each case that contributes a count to the table must be independent of all the other cases in the table.\n",
    "\n",
    "- Sample size or distribution: Each particular case (ie cell count) must have at least 5 expected cases.\n",
    "\n",
    "Let's see an example from the **pingouin** docs: whether gender is a good predictor of heart disease. First, let's load the data and look at the gender split in total:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi_data = pg.read_dataset(\"chi2_independence\")\n",
    "chi_data[\"sex\"].value_counts(ascending=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If gender is *not* a predictor, we would expect a roughly similar split between those who have heart disease and those who do not. Let's look at the observed versus the expected split once we categorise by gender and 'target' (heart disease or not)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected, observed, stats = pg.chi2_independence(chi_data, x=\"sex\", y=\"target\")\n",
    "observed - expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have fewer in the 0, 0 and 1, 1 buckets than expected but more in the 0, 1 and 1, 0 buckets. Let's now see how the test interprets this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From these, it is clear we can reject the null and therefore it seems like gender is a good predictor of heart disease."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shapiro-Wilk Test for Normality\n",
    "\n",
    "Note that the null here is that the distribution *is* normal, so normality is only rejected when the p-value is sufficiently small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = prng.normal(size=20)\n",
    "pg.normality(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test can also be run on multiple variables in a dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pg.read_dataset(\"ancova\")\n",
    "pg.normality(df[[\"Scores\", \"Income\", \"BMI\"]], method=\"normaltest\").round(3)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3.10.12 ('codeforecon')",
   "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.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "c4570b151692b3082981c89d172815ada9960dee4eb0bedb37dc10c95601d3bd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}