09 - Quantifying Scatter of Continuous Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantifying Scatter of Continuous Data\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In the realm of biostatistics, understanding the concept of scatter is crucial. Scatter refers to the **spread** of data points around the **central value**. It provides insights into the **variability** or diversity in a data set.\n",
    "\n",
    "In this chapter, we will explore various methods to quantify scatter, such as the **range**, **interquartile** range, **variance**, and **standard deviation**. These measures help us understand the **dispersion** in our data and are fundamental to many statistical tests and models.\n",
    "\n",
    "Here's a sneak peek into what we will cover with Python:\n",
    "\n",
    "1. *Range*: the simplest measure of scatter, defined as the difference between the maximum and minimum values in a dataset.\n",
    "2. *Interquartile Range (IQR)*: this measure gives us the range within which the central 50% of our data lies.\n",
    "3. *Variance* and *Standard Deviation*: these are more sophisticated measures of scatter that take into account how far each data point is from the mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "### Variance\n",
    "\n",
    "The **variance** is a measure of how much the data varies around its mean. There are two different definitions of variance.\n",
    "\n",
    "The **population variance** assumes that that we are working with the entire population and is defined as the average squared difference from the mean:\n",
    "\n",
    "$$\\operatorname{Var}_p(x) = \\sigma_p^2(x) = \\frac{1}{n}\\sum_{i = 1}^{n}(x_i - \\overline{x})^2$$\n",
    "\n",
    "More generally, the variance of a random variable $X$ is the expected value of the squared deviation from the mean of $X$, $\\mu =\\operatorname{E}[X]$:\n",
    "\n",
    "$$\\operatorname{Var}(X) = \\operatorname{E} \\left[(X - \\mu)^2\\right]$$\n",
    "\n",
    "The **sample variance** assumes that we are working with a sample and attempts to estimate the variance of a larger population by applying [*Bessel's correction*](https://en.wikipedia.org/wiki/Bessel%27s_correction) to account for potential sampling error. The sample variance is: \n",
    "\n",
    "$$\\operatorname{Var}_s(x)= \\sigma_s^2(x)  = \\frac{1}{n-1}\\sum_{i = 1}^{n}(x_i - \\overline{x})^2$$\n",
    "\n",
    "One can understand Bessel's correction as the degrees of freedom in the residuals vector (residuals, not errors, because the population mean is unknown): $(x_1 - \\bar{x}, \\dots, x_n - \\bar{x})$, where $\\bar{x}$ is the sample mean. While there are $n$ independent observations in the sample, there are only $n- 1$ independent residuals, as they sum to 0.\n",
    "\n",
    "We can see that $\\operatorname{Var}_p(x) = \\frac{n - 1}{n}\\operatorname{Var}_s(x)$, so as the data set gets larger, the sample variance and the population variance become less and less distinguishable, which intuitively makes sense.\n",
    "\n",
    "Note that variance can be decomposed to show that it is also the difference of the mean of squares to the mean squared:\n",
    "\n",
    "$$\\operatorname{Var}(x) = (\\frac{1}{n}\\sum_{i = 1}^{n}x_i^2) - \\overline{x}^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard Deviation\n",
    "\n",
    "Variance does not have intuitive scale relative to the data being studied, because we have used a 'squared distance metric', therefore we can square-root it to get a measure of 'deviance' on the same scale as the data. We call this the *standard deviation* $\\sigma(x)$, where $\\operatorname{Var}(x) = \\sigma^2(x)$\n",
    "\n",
    "The standard deviation accounts for the variation among the values, with an estimation of the spread of the distribution:\n",
    "\n",
    "$$\\sigma = \\text{SD} = \\sqrt{\\frac{\\sum(Y_i - \\overline{Y})^2}{n-1}}$$\n",
    "\n",
    "As with variance, standard deviation has both population and sample versions, and the sample version is calculated by default. Conversion between the two takes the form\n",
    "\n",
    "$$\\sigma_p(x) = \\sqrt{\\frac{n-1}{n}}\\sigma_s(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard Error of the Mean (SEM)\n",
    "\n",
    "The **Standard Error of the Mean (SEM)** is a measure of how far our sample mean is likely to be from the true population mean. SEM is calculated as the standard deviation divided by the square root of the sample size:\n",
    "\n",
    "$$\\mathrm{SEM} = \\frac{\\mathrm{SD}}{\\sqrt{n}}$$\n",
    "\n",
    "The SEM gets smaller as our samples get larger. This implies that as we collect more data, the estimate of the mean gets more precise, which makes sense intuitively. Now, the Confidence Interval is a range of values, derived from the SEM, that predicts the probability of a parameter (like the mean) to fall within that range as $W = z^* \\times \\text{SEM}$, as illustrated in the next chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practical examples using `numpy` statistical functions\n",
    "\n",
    "Let's calculate the mean, standard deviation, and other statistics for a given dataset using numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the dataset is 70.86\n",
      "The sample (ddof=1) variance of the dataset is 105.013\n",
      "The sample (ddof=1) standard error of the mean is 4.583\n",
      "The corresponding 95% confidence interval is (58.136, 83.584)\n",
      "The population standard deviation (ddof=0) is 9.166\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# The dataset\n",
    "data = np.array([64.7, 65.4, 88.3, 64, 71.9])\n",
    "\n",
    "# Calculate the mean\n",
    "mean = np.mean(data)\n",
    "print(f\"The mean of the dataset is {mean:.2f}\")\n",
    "\n",
    "# Calculate the variance\n",
    "variance = np.var(data, ddof=1) # ddof=0 by default\n",
    "print(f\"The sample (ddof=1) variance of the dataset is {variance:.3f}\")\n",
    "\n",
    "# Calculate the standard error of the mean\n",
    "sem = np.std(data, ddof=1) / np.sqrt(len(data))\n",
    "print(f\"The sample (ddof=1) standard error of the mean is {sem:.3f}\")\n",
    "\n",
    "# Calculate the 95% Confidence Interval\n",
    "ci = stats.sem(data) * stats.t.ppf((1 + 0.95) / 2., len(data)-1)\n",
    "lower_bound = np.mean(data) - ci\n",
    "upper_bound = np.mean(data) + ci\n",
    "print(f\"The corresponding 95% confidence interval is ({lower_bound:.3f}, {upper_bound:.3f})\")\n",
    "\n",
    "# Calculate the standard deviation with ddof=0\n",
    "std_dev = np.std(data, ddof=0)\n",
    "print(f\"The population standard deviation (ddof=0) is {std_dev:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ddof` parameter stands for 'Delta Degrees of Freedom'. When calculating the standard deviation, the sum of squared deviations from the mean is usually divided by `n - ddof`, where `n` is the number of elements in the dataset.\n",
    "\n",
    "Setting `ddof=0` means that we divide by n, which gives us the _population standard deviation_. This assumes that the dataset represents the entire population.\n",
    "\n",
    "On the other hand, if we set `ddof=1`, we would be dividing by $n - 1$, giving the _sample standard deviation_. This is used when our data is a sample from a larger population.\n",
    "\n",
    "In most cases, if we're working with a sample, we would use `ddof=1` to get an unbiased estimate of the actual population standard deviation. However, if our data represents the entire population, we would use `ddof=0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The 3-sigma rule\n",
    "\n",
    "The **3-sigma** (3-$\\sigma$) rule, also known as the _empirical rule_ or _68-95-99.7 rule_, is a statistical rule which states that for a normal distribution:\n",
    "\n",
    "1. Approximately 68% of the data falls within one standard deviation of the mean.\n",
    "2. Approximately 95% of the data falls within two standard deviations of the mean.\n",
    "3. Approximately 99.7% of the data falls within three standard deviations of the mean.\n",
    "\n",
    "This rule is a quick way to get an understanding of where most of the values in the dataset are likely to fall, assuming the data follows a normal distribution.\n",
    "\n",
    "Here's how it relates to the standard deviation:\n",
    "\n",
    "- The standard deviation is a measure of the amount of variation or dispersion in a set of values. A low standard deviation means that the values tend to be close to the mean, while a high standard deviation means that the values are spread out over a wider range.\n",
    "- When we talk about the \"3-sigma\" or \"three standard deviations from the mean\", we're talking about the range of values that covers about 99.7% of the data (assuming a normal distribution). This means that the probability of a value being outside of this range is very low (0.3%).\n",
    "\n",
    "So, the 3-sigma rule gives us a way to understand the standard deviation in terms of the percentage of data that falls within these ranges. It's a useful rule of thumb for understanding the distribution of the data when working with statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometric Mean and Standard Deviation\n",
    "\n",
    "The **geometric mean** is a type of average that is calculated by multiplying all the numbers together and then taking the nth root of the product. It's particularly useful when we want to compare things with very different properties, or when dealing with values that increase exponentially.\n",
    "\n",
    "The geometric standard deviation is a measure of spread of a set of numbers where the geometric mean is preferred. It expounds the factor by which a set of numbers deviate from the geometric mean. It's used when dealing with values that increase exponentially.\n",
    "\n",
    "If a dataset follows a **log-normal** distribution, the logarithm of the values will have a normal distribution. In a log-normal distribution, the multiplicative spread of values (i.e., the ratio between different data points) is more important than the absolute difference between them. This is where the geometric mean and geometric standard deviation come into play.\n",
    "\n",
    "For a log-normal distribution, about 2/3 of the values fall between the geometric mean divided by the geometric standard deviation and the geometric mean multiplied by the geometric standard deviation. This is analogous to the rule for a normal distribution, where about 2/3 of the values fall within one standard deviation of the mean.\n",
    "\n",
    "When we take the logarithm of a log-normally distributed dataset, it becomes normally distributed (this is the \"log converts multiplicative scatter to additive\" part). The geometric mean corresponds to the arithmetic mean of the logged data, and the geometric standard deviation corresponds to the standard deviation of the logged data. So, the rule that 2/3 of the data falls within one standard deviation of the mean in a normal distribution translates to the rule we mentioned for a log-normal distribution.\n",
    "\n",
    "This property is very useful in many fields, including finance, biology, and engineering, where variables often exhibit exponential growth and one is interested in rates or multiplicative factors rather than raw differences. For example in biology, we can use geometric mean and standard deviation for potency of drug (EC50, IC50, Km, Ki etc.), blood serum concentrations of many natural or toxic compounds, etc.\n",
    "\n",
    "While `numpy` provides a lot of useful statistical functions, it does not directly provide a function to calculate the geometric mean or the geometric standard deviation. However, we can calculate these using numpy's other functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The geometric mean is 70.319\n",
      "The geometric standard deviation is 1.1451\n"
     ]
    }
   ],
   "source": [
    "# Calculate the geometric mean\n",
    "gmean = np.exp(np.mean(np.log(data)))\n",
    "print(f\"The geometric mean is {gmean:.3f}\")\n",
    "\n",
    "# Calculate the geometric standard deviation\n",
    "gstd = np.exp(np.std(np.log(data), ddof=1))\n",
    "print(f\"The geometric standard deviation is {gstd:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weighted statistics\n",
    "\n",
    "Weighted statistics are a type of statistics where some data points contribute more than others. The **\"weight\"** of a data point refers to the relative importance or frequency of that point in the dataset.\n",
    "\n",
    "For example, if we're calculating the average test score in a class, but some tests count for more of the final grade than others, we would use a weighted mean, where each test score is multiplied by the weight of that test before summing and dividing by the total weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The weighted mean is 74.01\n"
     ]
    }
   ],
   "source": [
    "# The dataset\n",
    "data = np.array([64.7, 65.4, 88.3, 64, 71.9])\n",
    "\n",
    "# Weights for each data point\n",
    "weights = np.array([0.1, 0.2, 0.3, 0.1, 0.3])\n",
    "\n",
    "# Calculate the weighted mean\n",
    "weighted_mean = np.average(data, weights=weights)\n",
    "print(f\"The weighted mean is {weighted_mean}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical analysis with other packages\n",
    "\n",
    "### The `scipy.stats` module\n",
    "\n",
    "This module contains a large number of probability distributions as well as a growing library of statistical functions:\n",
    "\n",
    "- _Interquartile Range (IQR)_ is a measure of statistical dispersion, being equal to the difference between the upper and lower quartiles. IQR is a measure of variability, based on dividing a data set into quartiles.\n",
    "- _Descriptive summary_ provides several descriptive statistics of the data. This includes measures such as the number of elements in the data set, the minimum and maximum value, the mean, variance, skewness, and kurtosis.\n",
    "- _Z-score_ is a statistical measurement that describes a value's relationship to the mean of a group of values. It is measured in terms of standard deviations from the mean.\n",
    "- _Coefficient of Variation (CV)_ is a statistical measure of the relative dispersion of data points in a data series around the mean, with $\\text{CV} = \\sigma / \\mu$. For example, if $\\text{CV} = 0.25$, we know that the SD is 25% of the mean\n",
    "- _Geometric mean_ is a kind of average of a set of numbers that is different from the arithmetic average. The geometric mean is calculated by multiplying all the numbers together and then taking the nth root of the product.\n",
    "- _Geometric standard deviation_ is a measure of spread of a set of numbers where the geometric mean is preferred, it expounds the factor by which a set of numbers deviate from the geometric mean.\n",
    "\n",
    "Here's how we can use it to calculate various statistical measures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Interquartile Range (IQR) is 7.20\n",
      "The Standard Error of the Mean (SEM) is 4.583\n",
      "The descriptive summary is:\n",
      "DescribeResult(nobs=5, minmax=(64.0, 88.3), mean=70.86000000000001, variance=105.01299999999995, skewness=1.1912013156755001, kurtosis=-0.24972334599204293)\n",
      "The z-scores are [-0.6720695  -0.59569796  1.90274222 -0.74844103  0.11346628]\n",
      "The Coefficient of Variation (CV) is 0.1293\n",
      "The geometric mean is 70.32\n",
      "The geometric standard deviation is 1.145\n"
     ]
    }
   ],
   "source": [
    "# Interquartile Range (IQR)\n",
    "iqr = stats.iqr(data)\n",
    "print(f\"The Interquartile Range (IQR) is {iqr:.2f}\")\n",
    "\n",
    "# Standard Error of the Mean (SEM)\n",
    "sem = stats.sem(data)\n",
    "print(f\"The Standard Error of the Mean (SEM) is {sem:.3f}\")\n",
    "\n",
    "# Descriptive summary\n",
    "summary = stats.describe(data)\n",
    "print(f\"The descriptive summary is:\\n{summary}\")\n",
    "\n",
    "# Z-score\n",
    "z_score = stats.zscore(data)\n",
    "print(f\"The z-scores are {z_score}\")\n",
    "\n",
    "# Coefficient of Variation (CV)\n",
    "cv = stats.variation(data)\n",
    "print(f\"The Coefficient of Variation (CV) is {cv:.4f}\")\n",
    "\n",
    "# Geometric mean\n",
    "gmean = stats.gmean(data)\n",
    "print(f\"The geometric mean is {gmean:.2f}\")\n",
    "\n",
    "# Geometric standard deviation\n",
    "gstd = stats.gstd(data)\n",
    "print(f\"The geometric standard deviation is {gstd:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `pandas` library\n",
    "\n",
    "We can perform similar statistical calculations using the `pandas` library, which provides a DataFrame structure to hold and manipulate data. These methods are very similar to the ones provided by `numpy` and `scipy.stats`, but they are methods of the DataFrame or Series object instead of standalone functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the dataset is 70.86\n",
      "The standard deviation of the dataset is 10.248\n",
      "The Standard Error of the Mean (SEM) is 4.583\n",
      "The Interquartile Range (IQR) is 7.20\n",
      "The descriptive summary is:\n",
      "count     5.000000\n",
      "mean     70.860000\n",
      "std      10.247585\n",
      "min      64.000000\n",
      "25%      64.700000\n",
      "50%      65.400000\n",
      "75%      71.900000\n",
      "max      88.300000\n",
      "Name: Data, dtype: float64\n",
      "The z-scores are:\n",
      "0   -0.601117\n",
      "1   -0.532808\n",
      "2    1.701864\n",
      "3   -0.669426\n",
      "4    0.101487\n",
      "Name: Z_score, dtype: float64\n",
      "The Coefficient of Variation (CV) is 0.1446\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(data, columns=['Data'])\n",
    "\n",
    "# Calculate the mean\n",
    "mean = df['Data'].mean()\n",
    "print(f\"The mean of the dataset is {mean:.2f}\")\n",
    "\n",
    "# Calculate the standard deviation\n",
    "std_dev = df['Data'].std(ddof=1)\n",
    "print(f\"The standard deviation of the dataset is {std_dev:.3f}\")\n",
    "\n",
    "# Calculate the Standard Error of the Mean (SEM)\n",
    "sem = df['Data'].sem()\n",
    "print(f\"The Standard Error of the Mean (SEM) is {sem:.3f}\")\n",
    "\n",
    "# Calculate the Interquartile Range (IQR)\n",
    "iqr = df['Data'].quantile(0.75) - df['Data'].quantile(0.25)\n",
    "print(f\"The Interquartile Range (IQR) is {iqr:.2f}\")\n",
    "\n",
    "# Descriptive summary\n",
    "summary = df['Data'].describe()\n",
    "print(f\"The descriptive summary is:\\n{summary}\")\n",
    "\n",
    "# Z-score\n",
    "df['Z_score'] = (df['Data'] - df['Data'].mean()) / df['Data'].std(ddof=1)\n",
    "print(f\"The z-scores are:\\n{df['Z_score']}\")\n",
    "\n",
    "# Coefficient of Variation (CV)\n",
    "cv = df['Data'].std(ddof=1) / df['Data'].mean()\n",
    "print(f\"The Coefficient of Variation (CV) is {cv:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `statsmodels` library\n",
    "\n",
    "`DescrStatsW` is a class in the `statsmodels` library in Python that provides [descriptive statistics and statistical tests with weights for case weights](https://www.statsmodels.org/dev/generated/statsmodels.stats.weightstats.DescrStatsW.html). It assumes that the data is 1D or 2D with observations in rows, variables in columns, and that the same weight applies to each column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The weighted mean is 74.01\n",
      "The variance is 96.1109\n",
      "The standard deviation is 9.8036\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.stats.weightstats import DescrStatsW\n",
    "\n",
    "# Create a DescrStatsW object\n",
    "d1 = DescrStatsW(data, weights=weights, ddof=0)\n",
    "\n",
    "# Calculate the weighted mean\n",
    "mean = d1.mean\n",
    "print(f\"The weighted mean is {mean}\")\n",
    "\n",
    "# Calculate the variance\n",
    "variance = d1.var\n",
    "print(f\"The variance is {variance:.4f}\")\n",
    "\n",
    "# Calculate the standard deviation\n",
    "std = d1.std\n",
    "print(f\"The standard deviation is {std:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interlude: understanding key statistical formulas\n",
    "\n",
    "In this section, we will explore some fundamental statistical formulas that are crucial in data analysis. These include the mean of a sum, the variance of a linear combination, and the standard error of the difference between means.\n",
    "\n",
    "### Mean of a sum\n",
    "\n",
    "The mean, or average, is a measure of central tendency. If we have two random variables, $X$ and $Y$, the mean of their sum is simply the sum of their means. Mathematically, this is expressed as $E(X+Y)=E(X)+E(Y)$. We demonstrate this relation as follows:\n",
    "\n",
    "Let $x = x_1, x_2, \\dots, x_n$ and $y = y_1, y_2, \\dots, y_m$ be samples of two random variables of length $n$ and $m$ respectively. If $m = n$ and $x + y$ is formed from the element-wise sum of $x$ and $y$, it is obvious that the mean of $x + y$ is equal to the sum of the mean of $x$ and the mean of $y$:\n",
    "\n",
    "$$\\frac{1}{n}\\sum_{i=1}^n x_i + \\frac{1}{n}\\sum_{i=1}^n y_i = \\frac{1}{n}\\sum_{i=1}^n (x_i + y_i)$$\n",
    "\n",
    "### Variance of a linear combination\n",
    "\n",
    "The variance measures how spread out a set of data is. If we have a linear combination of two random variables, $X$ and $Y$, with constants $a$ and $b$, the variance is given by $\\mathrm{Var}(aX+bY) = a^2\\mathrm{Var}(X) + b^2\\mathrm{Var}(Y) + 2ab\\mathrm{Cov}(X,Y)$.\n",
    "\n",
    "We can first demonstrate that if all values are scaled by a constant, the variance is scaled by the square of that constant $\\operatorname{Var}(aX) = a^2 \\operatorname{Var}(X)$:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\operatorname{Var}(aX) &= \\frac{1}{n}\\sum_{i = 1}^{n}(ax_i - \\overline{ax})^2 \\\\\n",
    "&= \\frac{1}{n}\\sum_{i = 1}^{n}(ax_i - a\\overline{x})^2 \\\\\n",
    "&= \\frac{1}{n}\\sum_{i = 1}^{n}a^2(x_i - \\overline{x})^2 \\\\\n",
    "&= a^2\\left(\\frac{1}{n}\\sum_{i = 1}^{n}(x_i - \\overline{x})^2\\right) \\\\\n",
    "&= a^2\\operatorname{Var}(x)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The variance of a sum of two random variables is given by $ \\operatorname {Var}(aX+bY) = a^2 \\operatorname{Var}(X) + b^2 \\operatorname{Var}(Y) + 2ab \\operatorname{Cov}(X,Y) $, and $ \\operatorname {Var}(aX-bY) = a^2 \\operatorname{Var}(X) + b^2 \\operatorname{Var}(Y) - 2ab \\operatorname{Cov}(X,Y) $ which can be simplified in the case of independent variables, in which case the covariate term drops out:\n",
    "\n",
    "$$\\operatorname{Var}(X \\pm Y) = \\operatorname{Var}(X) + \\operatorname{Var}(Y)$$\n",
    "\n",
    "### Standard Error of the difference between means\n",
    "\n",
    "The standard error measures the precision of an estimate. If we have two independent random variables, X and Y, the standard error of the difference between their means is given by $\\mathrm{SE}(\\hat{X} - \\hat{Y}) = \\sqrt{\\mathrm{SE}(\\hat{X})^2 + \\mathrm{SE}(\\hat{Y})^2}$.\n",
    "\n",
    "Let's first derive the formula for the standard error of the mean of $ x = x_1, x_2, \\dots, x_n $:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\operatorname{Var}(\\overline{x}) &= \\operatorname{Var}\\left(\\frac{1}{n}\\sum_{i = 1}^{n}x_i\\right) \\\\\n",
    "&= \\sum_{i = 1}^{n}\\operatorname{Var}\\left(\\frac{1}{n}x_i\\right) \\\\\n",
    "&= \\frac{1}{n^2}\\sum_{i = 1}^{n}\\operatorname{Var}\\left(x_i\\right) \\\\\n",
    "&= \\frac{1}{n^2} n \\operatorname{Var}\\left(x\\right) \\\\\n",
    "&= \\frac{1}{n}\\operatorname{Var}\\left(x\\right) \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Since each $x_i$ is independent and identically distributed, the standard error of the mean is the standard deviation of the mean, which is the square root of the variance of the mean. So taking the square root of the previous derivation:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "SE &= \\sqrt{\\operatorname{Var}(\\overline{x})} \\\\\n",
    "&= \\frac{\\sqrt{\\operatorname{Var}(x)}}{\\sqrt{n}} \\\\\n",
    "&= \\frac{\\sigma(x)}{\\sqrt{n}}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Finally, the standard error of the difference between the means of $x$ and $y$ (see chapter 30 - Comparing means) can be derived to\n",
    "\n",
    "$$\\mathrm{SE}_\\mathrm{mean\\_diff} = \\sqrt{\\frac{\\sigma^2(x)}{n} + \\frac{\\sigma^2(y)}{m}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In this chapter, we delved into the world of descriptive statistics and explored various measures of **central tendency** and **dispersion**. We started with the concept of scatter and how it is quantified using measures like range, interquartile range, variance, and standard deviation. We then discussed the importance of these measures in understanding the spread and variability of our data.\n",
    "\n",
    "We learned how to calculate these measures using Python and its powerful libraries like `numpy`, `scipy`, and `pandas`. We also touched upon the concept of weighted statistics and how some data points can contribute more than others based on their relative importance or frequency. We also introduced the `statsmodels` library and its `DescrStatsW` class for weighted descriptive statistics, and demonstrated how to use it to calculate various statistical measures.\n",
    "\n",
    "Finally, we took a mathematical interlude to understand key statistical formulas, including the mean of a sum, the variance of a linear combination, and the standard error of the difference between means. These formulas provide the mathematical foundation for many statistical tests and measures.\n",
    "\n",
    "Understanding these concepts and being able to calculate and interpret these measures are fundamental skills in data analysis and statistics. They allow us to summarize and understand our data, and form the basis for many advanced statistical tests and models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cheat sheet\n",
    "\n",
    "This cheat sheet provides a quick reference for essential code snippets used in this chapter.\n",
    "\n",
    "### Statistics with `numpy`\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "\n",
    "# Defining the sample dataset\n",
    "data = np.array([64.7, 65.4, 88.3, 64, 71.9])\n",
    "\n",
    "# Calculate basic statistics\n",
    "mean = np.mean(data)\n",
    "var = np.var(data, ddof=1) # ddof=0 by default\n",
    "std = np.std(data, ddof=1)\n",
    "sem = std / np.sqrt(len(data))\n",
    "\n",
    "# Geometric mean and std\n",
    "gmean = np.exp(np.mean(np.log(data)))\n",
    "gstd = np.exp(np.std(np.log(data), ddof=1))\n",
    "\n",
    "# Weighted statistics\n",
    "# Weights for each corresponding point in data\n",
    "weights = np.array([0.1, 0.2, 0.3, 0.1, 0.3])\n",
    "\n",
    "weighted_mean = np.average(data, weights=weights)\n",
    "```\n",
    "\n",
    "### Statistics with `scipy.stats`\n",
    "\n",
    "```python\n",
    "from scipy import stats\n",
    "\n",
    "# Basic statistics\n",
    "summary = stats.describe(data) # desciptive summary\n",
    "iqr = stats.iqr(data)\n",
    "sem = stats.sem(data)\n",
    "zscore = stats.zscore(data)\n",
    "cv = stats.variation(data)\n",
    "\n",
    "# Geometric mean and std\n",
    "gmean = stats.gmean(data)\n",
    "gstd = stats.gstd(data)\n",
    "```\n",
    "\n",
    "### Statistics with `pandas`\n",
    "\n",
    "```python\n",
    "import pandas as pd\n",
    "\n",
    "# Define the dataframe\n",
    "df = pd.DataFrame(data, columns=['values'])\n",
    "\n",
    "# Basic statistics\n",
    "summary = df['values'].describe() # desciptive summary\n",
    "mean = df['values'].mean()\n",
    "std = df['values'].std(ddof=1)\n",
    "sem = df['values'].sem()\n",
    "iqr = df['values'].quantile(0.75) - df['values'].quantile(0.25)\n",
    "df['zscore'] = (df['values'] - df['values'].mean()) / df['values'].std(ddof=1)\n",
    "cv = df['values'].std(ddof=1) / df['values'].mean()\n",
    "```\n",
    "\n",
    "### Statistics with `statsmodels`\n",
    "\n",
    "```python\n",
    "from statsmodels.stats.weightstats import DescrStatsW\n",
    "\n",
    "# Create a DescrStatsW object\n",
    "d = DescrStatsW(\n",
    "    data,\n",
    "    weights=weights, # 'None' if no weights\n",
    "    ddof=1)\n",
    "\n",
    "# Basic statistics\n",
    "mean = d.mean\n",
    "variance = d.var\n",
    "std = d.std\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Session Information\n",
    "\n",
    "The output below details all packages and version necessary to reproduce the results in this report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python 3.12.4\n",
      "-------------\n",
      "numpy: 1.26.4\n",
      "pandas: 2.2.2\n",
      "scipy: 1.13.1\n",
      "statsmodels: 0.14.2\n"
     ]
    }
   ],
   "source": [
    "!python --version\n",
    "print(\"-------------\")\n",
    "# List of packages we want to check the version\n",
    "packages = ['numpy', 'pandas', 'scipy', 'statsmodels']\n",
    "\n",
    "# Initialize an empty list to store the versions\n",
    "versions = []\n",
    "\n",
    "# Loop over the packages\n",
    "for package in packages:\n",
    "    # Get the version of the package\n",
    "    output = !pip show {package} | findstr \"Version\"\n",
    "    # If the output is not empty, get the version\n",
    "    if output:\n",
    "        version = output[0].split()[-1]\n",
    "    else:\n",
    "        version = 'Not installed'\n",
    "    # Append the version to the list\n",
    "    versions.append(version)\n",
    "\n",
    "# Print the versions\n",
    "for package, version in zip(packages, versions):\n",
    "    print(f'{package}: {version}')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "35e05afc3bbf8458d25376a416e873cf169d04b2fc4efc2b77204c35b02a6d38"
  },
  "kernelspec": {
   "display_name": "Python 3.9.6 64-bit",
   "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.12.5"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}