{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Unbiased and Biased Distance Correlation (Dcorr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we explore\n", "\n", "- The theory behind the Dcorr test statistic and p-value and it's equivalency with HSIC\n", "- The features of the implementation\n", "- A fast implementation of Dcorr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Formulation of test statistic and p-value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following description is adapted from [[1]](https://arxiv.org/abs/1907.02088):\n", "\n", "Dcorr can be used to measure linear and nonlinear associations between two random vectors of arbitrary dimension. Let ${\\mathbf{D}}^{\\mathbf{x}}$ be the $n \\times n$ distance matrix of $\\mathbf{x}$ and ${\\mathbf{D}}^{\\mathbf{y}}$ be the $n \\times n$ distance matrix of $\\mathbf{y}$. The distance covariance (Dcov) is defined as\n", "\n", "$$\\mathrm{Dcov}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{1}{n^2} \\mathrm(tr) \\left( {{\\mathbf{D}}^{\\mathbf{x}} \\mathbf{H} {\\mathbf{D}}^{\\mathbf{y}} \\mathbf{H}} \\right).$$\n", "\n", "The normalized version of this covariance is Dcorr [[2]](https://projecteuclid.org/euclid.aos/1201012979) and can be calculated as\n", "\n", "$$\\mathrm{Dcorr}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{\\mathrm{Dcov}_n (\\mathbf{x}, \\mathbf{y})}{\\sqrt{\\mathrm{Dcov}_n ( \\mathbf{x}, \\mathbf{x}) \\mathrm{Dcov}_n (\\mathbf{y}, \\mathbf{y})}}.$$\n", "\n", "The original Dcorr is biased. A modified matrix ${\\mathbf{C}}^{\\mathbf{x}}$ based on $\\mathbf{H} {\\mathbf{D}}^{\\mathbf{x}} \\mathbf{H}$ can be calculated as\n", "\n", "$${\\mathbf{C}_{ij}}^{\\mathbf{x}} =\n", "\\begin{cases}\n", "{\\mathbf{D}_{ij}}^{\\mathbf{x}} - \\frac{1}{n - 2} \\sum_{t = 1}^n {\\mathbf{D}_{it}}^{\\mathbf{x}} - \\frac{1}{n - 2} \\sum_{s = 1}^n {\\mathbf{D}_{sj}}^{\\mathbf{x}} + \\frac{1}{\\left( n - 1 \\right) (n - 2)} \\sum_{s, t = 1}^n {\\mathbf{D}_{st}}^{\\mathbf{x}}& i \\neq j \\\\\n", "0 & \\mathrm{otherwise}\n", "\\end{cases},$$\n", "\n", "and similarly for ${\\mathbf{C}}^{\\mathbf{y}}$. Then\n", "the unbiased Dcov (UDcorr) [[3]](https://projecteuclid.org/euclid.aos/1413810731) is thus\n", "\n", "$$ \\mathrm{UDcov}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{1}{n (n - 3)} \\mathrm{tr} ({{\\mathbf{C}}^{\\mathbf{x}} {\\mathbf{C}}^{\\mathbf{y}}}).$$\n", "\n", "UDcorr can consequently be calculated using the same equation as the one used to Dcorr." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Equivalency of Dcorr and HSIC" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following description is adapted from [[1]](https://arxiv.org/abs/1907.02088):\n", "\n", "Hsic is a way to measure multivariate nonlinear associations given a specified kernel [[4]](https://arxiv.org/abs/0805.2368). Let ${\\mathbf{K}}^{\\mathbf{x}}$ be an $n \\times n$ kernel matrix of $\\mathbf{x}$ and ${\\mathbf{K}}^{\\mathbf{y}}$ be an $n \\times n$ kernel matrix of $\\mathbf{y}$. In addition, $\\mathbf{H} = \\mathbf{I} - (1 / n) \\mathbf{J}$ and is an $n \\times n$ centering matrix, where $\\mathbf{I}$ is the identify matrix of size $n \\times n$ and $\\mathbf{J}$ is a matrix of ones that is the same size. Then, Hsic is defined as\n", "\n", "$$\\mathrm{Hsic}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{1}{n^2} \\mathrm{tr} ({{\\mathbf{K}}^{\\mathbf{x}} \\mathbf{H} {\\mathbf{K}}^{\\mathbf{y}} \\mathbf{H}}).$$\n", "\n", "The default kernel choice is the Gaussian kernel using the median distance as the bandwidth, which is a characteristic kernel that guarantees Hsic being a consistent test [[4](https://arxiv.org/abs/0805.2368), [5](http://www.jmlr.org/papers/v11/gretton10a.html)]. Specifically, the Gaussian kernel transformation for $\\mathbf{x}$ of the form $\\mathcal{N} (0, \\mathrm{med}{ ( x_i - x_j) })$ where $\\mathrm{med}$ refers to the dimension-wise median of $\\mathbf{x}$. The transformation for $\\mathbf{y}$ is defined similarly.\n", "\n", "Hsic and Dcorr are indeed equivalent in the sense that every valid kernel has a corresponding valid semimetric to ensure their equivalence, and vice versa [[6](https://projecteuclid.org/euclid.aos/1383661264), [7](https://arxiv.org/abs/1806.05514)]." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using Dcorr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before delving straight into function calls, let's first import some useful functions, to ensure consistency in these examples, we set the seed:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt; plt.style.use('classic')\n", "import matplotlib.ticker as ticker\n", "import seaborn as sns; sns.set(style=\"white\")\n", "\n", "from mgcpy.independence_tests.dcorr import DCorr\n", "from mgcpy.benchmarks import simulations as sims\n", "\n", "np.random.seed(12345678)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, let's simulate some linear data:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAIlCAYAAABrdaqpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhkdXn3//fMwAzOzrArJLjE+0FQgaBBlGV8CD/RbBqXRAOiQYxRRAGNRhQwxLgAsgdRgkGIQfQxRmNcogyKCBgXRIRbRVBEBphpppm1Z+n+/XGqhpqa6u7q01VdVV3v13X1Vd2nTp2+qyimPv1dZ4yMjCBJkjRRMztdgCRJ6k2GCEmSVIohQpIklWKIkCRJpRgiJElSKYYISQIiYkana5B6zQ6dLkCa7iJiGbBTZh46xjlnAWcCT8jMDVNUWmkRsTPwd8CfAfsCQ8BdwDXA5Zm5uXLeUcANwLGZ+ZUpqOuTwIsyc88JPGYG8F5gA/DhyrGz6KH/HlKn2BIhdYdPAM+j+DDuahHxBODbwCuBi4AXA68GbgY+ShEkqn5A8by+O8VlTsQc4Gxgbs2xnvnvIXWSLRFSF8jM3wC/6XQdTXo5sD/w7Mz8cc3x/4qIVcDZEfHBzPxRZj4G3NKRKiehx/57SB1jiJC6QH3zeaVZfl/gCuAM4KnAfcAHMvNfax43u3L/ccATgV8Bl2XmBXXXPwF4E8WH/47ALyvnXVy5f1/gXuBU4PXAPsA5mXlug3KrXQWNWjKvBDYCqyvXPYqa7oxKHVcBh1K0WhwEPAScBXwFuBg4FngMuCoz31NX35sy8/LRXrf6YiJiJvB24Hjg94AZwN3ABzPzuprrApwZEWdm5oxG142Ig4B/AJ5L0WrxXeC9mXlLXY2vAv4UeAkwC/hv4K2ZubzB6yX1NLszpO51IEUz+z8BfwT8GvhkRDyz5pzrgdOAj1XO+SxwXkT8U/WEiHgj8C/A14E/pmhJuA+4qPIhX+uDlWsdD3xplLq+AmwGvhoR/xgRR1S6OMjMBzLzg5l5zzjP7XPApyr13EfRfbAM+BnwMuAbwN9HxMvHuc54/rHydTVFt8tfUYScayPiqcCDwJGVc6+k6MLYTkQcCdwKLAT+huL1mQvcGBGH153+MeBRitf5XZXneOkkn4fUlWyJkLrXIuDwzLwDICKSoqXhT4A7IuKFle9fn5lXVR7z9YjYALwvIi6tNMs/DbgwM8+oXjgivgOsBJZSfHhXfTEzLxmrqMy8IyL+HLgc+PvK16aIuBW4DrgiMzeO89zOy8x/rtSyqVLD92paHm4AXgE8nyIYlbUPcGZmnlc9EBH3At8HjsjMqyLitspdv6m2KjTwIeB+4Ojqc4uI/6Jo1TgX+IOac7+ZmW+pfP8/EfH7wF9FxIzMdLMiTSuGCKl7ra4GiIpqH/28yu0fVm7/MyJq/1/+AkULxv8F/jUz3wEQEQuBp1OEit+vnDu77nfeQRMy8z8j4r+BIyq/53CKD9IXAG+MiKMyc+UYl/hOzfcPVW63foBn5qaIGAR2bqaeMer8K4CI2IXHn/vSyt31z72hiJhH0YXx4dpwlJlDEfEZ4B0RMb/mId+pu8RvKLqQdqRoBZGmDUOE1L3W1f6QmcMRAY93Q+5auV0xyuOfBBARTwb+GTgG2ELRZVD9oKtfG+EhmpSZmyi6Hb5R+T0LgXdTNOH/HfDOMR7+WINja+t+nvRf7RFxMHAJj8+0uIvHg1Kz60IsrpzbaEzDg5X7FtYcW1d3znDl1u5jTTu+qaXetYriL9vnAs9p8PXJyhoIXwJ+l6KVYF5m7g+8rewvjYibI+K6+uOZ+VhmvpsipOxf9vqjqAaKWXXHF4z2gIhYQDF+Y5hifMm8zDyIYtzHRKyq/P5Ga088sXLfWK0u0rRliJB61zKKJvm5mfm/1S+Kv5zPAfYCdgOeAXwyM2+uaY5/ceW2zL8B9wB/GhHPqL8jIhYDuwM/3u5Rk1Ntudin7nj9oMZa+1E8/4sz8/bM3FI5Xv/ct2z3yBqZuRa4DXhFROxYPR4RcyjWyrg1M11PQn3J7gxpauwZEY3++n8wM7f7q75J/00xffL6iPhH4HaKD85zKJrZf1Lpt78X+JvK7SMUH7x/R/EX9LyGVx7b3wNHATdHxKUUC0+tq/zutwEDwHmjPrqEzHw0Im6ieB53AQ9QTEV98hgPuxsYBN4VEesrNb4YqA56nFe59qaIWAMcFhFHVJ5PvXcDXwO+EREfrRw7laLL6ITJPDepl9kSIU2N36VYF6H+6+1lL5iZwxTTOq+qXOdrFB/w1wNLa/46/hOKdSE+QTG18iXAiRQh5IgSv/d+ivUdPk6xHsL1FNNHTwe+DDw3M0cbpzEZr6UYy3EZ8O8UYeXvxqjzMYrnPgR8uvJ1MMVr9lMen9oJxfoPz6F4TepbO8jMG4AXUkxtvQb4JMUYjiMyc9mknpXUw2aMjDjjSJIkTZwtEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKmWHThcwnoj4GDArM08c45zrgZfXHf5GZh7d1uIkSepjXRsiImIGcDZwEnDlOKcfALwL+NeaY0NtKk2SJNGlISIinkIRHA4Afj3OubOBpwG3ZebyKShPkiTRvWMingf8EngmcO845+5HEYbuandRkiTpcV3ZEpGZ1wLXAkTEeKcfAGwEzo6IY4H1wPXAOZm5oZ11SpLUz7oyREzQ/sAMIIFLKFovzgf2AV7bzAUiYg7wHOBBYEt7ypQkqWvMAvYCvpeZpccQTocQcQZwbmYOVH6+IyK2AP8eEadm5somrvEc4Nttq1CSpO50OHBT2Qf3fIjIzGFgoO7wHZXbfYBmQsSDANdeey177rlnC6uTJKn7LF++nNe85jVQ+fwrq+dDRER8BtgxM19ac/gQiimev2jyMlsA9txzT/bee+8WVyhJUteaVBd+z4WIypTOJcBAZm4EPkul6wL4AnAQcC5FF8eazlUqSdL01q1TPMdyGEXzy2EAmfkZ4ATgdcBPgPOAC4H3dag+SZL6Qte3RGTmUXU/L6OYjVF77Grg6qmrSpIk9WJLhCRJ6gKGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklbJDpwsYT0R8DJiVmSeOcc4hwIXAQcADwD9k5tVTVKIkSX2pa1siImJGRLwfOGmc83YDvgr8ADgYuAi4MiKOaX+VkiT1r65siYiIpwBXAgcAvx7n9BOBQeCUzBwG7o6Ig4HTga+1tVBJkvpYt7ZEPA/4JfBM4N5xzj0c+FYlQFQtA54fEd36/CRJ6nld2RKRmdcC1wJExHin7w38sO7Yb4G5wBJgRavrkyRJ3dsSMRFzgQ11x4YqtztNcS2SJPWN6RAi1gNz6o5Vf147xbVIktQ3pkOIuB/Yq+7YE4E1FAMuJUlSG0yHEHETcEREzKg5thT4Tt1gS0nSNLBleISBwQ3cde9KBgY3MDw80umS+lZXDqwcS0TMphgwOZCZGymmgr4TuDwiLgCOBl4NvKhzVUqS2mVw9RAnn3cDj63dyMJ5s7n4tKUsWeQQuE7oxZaIw4AHK7dk5kMUgeEgilkabwGOz8xvdqxCSVLbPDSwlsfWbgTgsbUbeWjA4W+d0vUtEZl5VN3Py4AZdcduAZ47dVVJkjpljyXzWDhv9taWiD2WzOt0SX2r60OEJEm1Fi+Yw8WnLeWhgbXssWQeixfUT9DTVDFESJJ6ysyZM1iyaCfHQXSBXhwTIUnqQs6a6D+2REiSWsJZE/3HlghJUtPGam1w1kT/sSVCktS0sVobnDXRfwwRkqSmNWptqIYIZ030H0OEJKlpY7U2tHLWxJbhEQZXD20TSGbOnDH+AzWlDBGSpKZNVWuDgzR7gyFCktS0qVqjYaxuE3UPZ2dIkrpOtdsEcJBmF7MlQpLUdRyk2RsMEZKkruPS1r3B7gxJklSKIUKS1JB7YWg8dmdIkhpymqXGY0uEJKkh98LQeAwRkqSGnGap8didIUlqyGmWGo8hQpLUkNMsNR67MyRJUimGCEmSVIohQpIklWKIkCRNSLOLULlY1fTnwEpJElB86A+uHtpmNsbMmTO2O6/ZRahcrGr6M0RIkoDmP/QbLUI1mfPUu+zOkCQBza9Q2ewiVC5WNf3ZEiFJAmD3JXNZOG/21paIPZbMbXhes4tQuVjV9GeIkCQBMGvmDM5+w/N4+NF17L7z3IbjIaD5RahcrGr6M0RIkgB4cMVa3nXpTcyfO5s16zbywTe/gMULDAAanSFCkvpQo5kYeyyZx/y5td0ZjmHQ2AwRktSHGs3EcAyDJsoQIUl9aLTpl45h0EQ4xVOS+pDTL9UKtkRIUh/q5q6LZlfOVOcZIiSpD3Xz9MvB1UO85/Lv8NtH1jB/rstldzNDhCT1gH7563zL8Aibtgxz3LH7sdviJ3D55+9wuewuZoiQpB7QL5tZDa4e4tQLbtz6PM9+w/NYsnD6Pc/pwoGVktQDmt3XotfVP8+hTVu6aryGtmWIkKQe0OuzKbYMjzAwuIG77l3JwOAGhodHGp5X/zz32mXetOy2mS7szpCkHtDNsyma0Wx3TK8/z35jiJCkHtDNsymaMdriVvV6/Xn2G7szJElt1+vdMWrMlghJUtvZTTE9GSIkSW1nN8X0ZHeGJEkqxRAhSZJKMURIkqRSDBGSpKY0u2CU+ocDKyVJTemX/TvUPFsiJElN6Zf9O9Q8Q4QkdVC3dBE0U4cLRqme3RmS1EHd0kXQTB0uGKV6hghJ6qBm95TohjpcMEr17M6QpA5qZRfBZLpG7KpQGbZESFIHtaKLYMvwCIOrh9i0ZZhTL7ixVNeIXRUqwxAhSR3Uii6CwdVDvOfy73Dcsftt7ZJYs24jK1ata/q6dlWoDEOEJHWBamtCbUvAzJkzmnrsQwNr+e0ja9ht8RNYOG82e+4yjze97FkMbdrCwOCGCV1LmghDhCR1gcnM0thjyTzmz53N5Z+/g7Pf8DwWzJtdultDmghDhCR1gcnM0qgdz7DLop1YvrI7Znxo+jNESFIXqM6OqLYeTGR2RP14hpERSl9LmghDhCR1gWZnRzQzdsKZFpoqXRkiImIWcA5wArAA+Arw5sx8aJTzrwdeXnf4G5l5dDvrlKRWaXZ2RDNjJ5xpoanSrYtNnQW8FjgeOALYG/jcGOcfALwL2Kvm6xXtLVGSttfuvTDcBEvdpOtaIiJiNnAK8NbM/Hrl2F8A90bEYZl5c4PznwbclpnLp7xgSarR7r0wJjN2Qmq1rgsRwIEUXRjLqgcy876IuA84HLi57vz9KJ7HXVNTniSNrt17YTjeQd2kG0PE3pXbB+qO/xbYp8H5BwAbgbMj4lhgPXA9cE5mbmhblZLUQLtbChzvoG7SjSFiLjCcmZvqjg8Bjf6v2R+YASRwCfBM4HyKwPHaNtYpSdsZr6VgMitTSt2mG0PEemBmROyQmZtrjs8BGo0gOgM4NzMHKj/fERFbgH+PiFMzc2Wb65WkrcZrKWj3mAlpKnVjiLi/crtXzfcAT2T7Lg4ycxgYqDt8R+V2H8AQIanjtgyP8NiaIR54ZLWrSWra6MYpnrcDq4EjqwciYl9gX+Bb9SdHxGci4vN1hw+h6P74RduqlKQJqLZA7DR7BxbOmw3Q0jET7Z5aKjXSdS0RmTkUEZcB50bECuBh4DLgxsy8pTKlcwkwkJkbgc9S6boAvgAcBJxL0cWxpjPPQpK29dDAWgbXbOTyz9/BWSceyoaNm3nSbgsazq4oM27CbhJ1Qje2REAxzuFa4BrgBuBXPL4i5WHAg5VbMvMzFCtbvg74CXAecCHwvimtWJLGUJ218bNfP8rZV97C3rsvYMminRqGg2ogeOclN3HyeTewavXQuNd3ESp1Qte1RABUBlSeVvmqv28ZxWyM2mNXA1dPSXGSVEL9rI1F80df36HMWhMuQqVO6MoQIUm9arSuiIms71AmELgIlTrBECFJLdSKsQllAoGLUKkTDBGS1EKtWPZ6vEDgglXqFoYISWqhqRib4EwMdQtDhCS10GTHJjTTytDuTb6kZhkiJKmFJjs2oZlWBmdiqFsYIiSphRq1JIxA02MYxmtl2DI8wsyZcP7bjmTFqvXstctcZ2KoYwwRktRCjVoSgG2OXXTaUcxgRsNQMV4rQ/X6a9ZtZK9d5/OBNz3fQZXqGEOEJLVQfUvCY2uHWLth8zbHlq9cx8Wf+REPPLJmuy6L8cZU1F7/gUfWOB5CHdWty15LUk+qtiQALJo/m/lzZzNnx1nbbLq16+In8OCKYmuf+iWqq2Mq9nvyLg2Xxa69vuMh1Gm2REhSC9W2JOy5yzyWr1zLJ/7zTs468VAeWbWepzxpEUObNjN/brmBka5MqW5iiJCkFqqfnTEyAstXruX0i77FXrvO57RXH8zHPn8H57/tSAYG129tWRgY3NDUwEtXplQ3MURIUhtVWw4eXLmWOTvO4p//349ZvnItO86ayX5P3gWAgcENLh6lnmSIkKQarV5SutpysHjBHFatHuLEP9l/u24IF49SrzJESFKNdi0pPVY3hItHqVcZIiSpRidaBRwsqV5liJCkGp1oFXCwpHqVIUKSaky0VcBtudXPDBGSVGOirQJuy61+5oqVkvrGluERBgY3cNe9KxkY3MDw8Mikr9loDIXUL2yJkNQ32tFq4MwK9TNDhKS+0Y6ZF86sUD8zREjqG+1oNSg7s8IBmZoODBGS+kY3tRo4IFPTgSFCUt/o5HoM9S0Pq9ZscKlr9TxDhCRNgfqWhwtPPYpF82czuMYBmepdhghJmgL1gzofeXQdF5+2lOUrO9+1IpVliJCkKdBoUOfOC3di54V2Yah3GSIkaQp006BOqVUMEZI0BdxkS9PRqMteR8Sbp7IQSZLUW8baO+PCiPifiNhnyqqRJEk9Y6wQ8QJgD+AnEfHXU1SPJEnqEaOGiMy8BTgI+AhwSUR8OSKeNGWVSZKkrjbmwMrM3AycExGfBi4E7oiIy4B1ded9oH0lSlL3ce8LqfnZGb8BvgccA7weGKq5bwQwREjqK+59ITURIiLixcClwC7A2zPz0rZXJUldrh3biku9ZtQQERG7AxcBrwD+BzgpM381VYVJUjdrx7biUq8ZqyXi7srtiZl51VQUI0m9whUopbFDxI3AmzJz+VQVI0m9whUopTFCRGa+dCoLkaRu4cwLqTnunSFJdZx5ITVnrBUrJakvNZp5IWl7hghJqlOdeTFzBuy9+3xnXkijsDtDkuosXjCHS09fytDmYVasWg+MMDw84rgIqc5Y60R8s9mLZOYLW1OOJHXezJkzGB6BUy+40XER0hjG6s64p+brQeAoYB5wB/B9YAZwOPDT9pYoSVPPcRHS+Maa4vmG6vcR8SngQ5n57tpzIuIM4LntK0+SOsMVKaXxNTsm4qUU24LXuw74+9aVI0ndwRUppfE1GyIepui6+Hnd8WOB+1takSR1AVeklMbXbIg4D7gsIg4FfkAxHuIw4FXACe0pTZIkdbOmQkRmXhoRjwFvBl4NjAA/Al6RmV9oY32SJKlLNb1ORGZ+CvhUG2uRJEk9pOkQERFPA04D/g/wV8CfAj/NzGXtKU2S2sdNtqTJa2rZ64j4A+B24KkUYyHmAPsDX4+IP2pfeZLUHtVNtt55yU2cfN4NrFo91OmSpJ7T7N4ZHwQ+mJnHABsBMvPNwIeBs9pTmiS1j4tJSZPXbIg4GPh0g+NXAvu1rhxJmhrVxaQAF5OSSmp2TMRaYHfgF3XHnw4MtrQiSZoCLiYlTV6zIeLTwEcj4gSK6Z07RcQLgUuAz7apNklqGxeTkiav2e6MdwP3AXcC8yk24fo68L3KfZLUFbYMjzAwuIG77l3JwOAGhodHOl2SNG01u9jURuBVEfEe4ECKwZV3ZuY97SxOkiaqOuvCLbyl9msqRETEL4FDMvMX1IyLiIi9gNszc/c21SdJ2xlrjYdGsy4MEVJ7jBoiIuLFwCGVH/cF3hURa+pOe/pY1ygrImYB51Dsy7EA+Arw5sx8aJTzDwEupNhp9AHgHzLz6lbXJak9RgsFtcf32mUewyMjPDSwjt12nsuF1/2AH/1sxXatDW7hLU2dsQLAvcAFFJttAbwc2FJz/wiwGji5DXWdBbwWOB5YCVwGfA54Qf2JEbEb8FXg34C/Bv4QuDIilmfm19pQm6QWG60Lonp8zbqNnHfKkZz58e9uPeesEw/lxz//1natDc66kKbOqCEiM++iaGkgIm4AXpaZj7a7oIiYDZwCvDUzv1459hfAvRFxWGbeXPeQEymmmZ6SmcPA3RFxMHA6YIiQesBoXRDV44vmz+bhR9dtc86KwfXMn1us81Db2uCsC2nqNDU7IzOXNgoQETE7Ip7f4poOpOjCWFbz+++jmB1yeIPzDwe+VQkQVcuA50dEs7NPJHXQaAs/VY+vXruR3Xeeu805v7fPzpzxuudy8WlLbW2QOqTZgZW/D3wceCaNg8esFta0d+X2gbrjvwX2GeX8HzY4dy6wBFjRwtoktcFoXRC1x3ddvBMXn3YUDw2s23rOrouf0OHKpf7W7KDIC4H1wEkU4xNOAZ5cuT2+xTXNBYYzc1Pd8SGgUfvkXGBDg3MZ5XxJLdSK3TBH64JodHzJom2Dg7txSp3TbIg4CDgiM78fEScBmZlXRMRvgTfR2lUr1wMzI2KHzNxcc3wOxfLbjc6vb8us/uyOOlKbdXpdhk7/fqmfNTtmYAbwSOX7n1N0awB8EXh2i2u6v3K7V93xJ7J9F0f1/EbnrsF9PaS26/RumJ3+/VI/azZE/AR4ceX7nwLVwZR70NrxEAC3U0wdPbJ6ICL2pVir4lsNzr8JOCIiatsvlwLfqRtsKakNOr0bZqd/v9TPmu3O+BBwXURsodiM68yI+A+KVogbWllQZg5FxGXAuRGxAniYYhzGjZl5S2UK6BJgoLIc95XAO4HLI+IC4Gjg1cCLWlmXpMbasS7DRMY5uC6E1DnNTvH8HHAocFtm/oqiVWIj8GXgDW2o6wzgWuAaipDyK4rFrgAOAx6s3FJZxfJFFOM2fgi8BTg+M7/Zhrok1akOftzvybuwZNFOLRnUWB3n8M5LbuLk825g1eqhUc9tx++X1Jyml6zOzP+t+f4GWtwCUfe7NgOnVb7q71vG46toVo/dAjy3XfVImlrufyH1hrH2zmh6tcfMPKY15UiaCt0+LdL9L6TeMFZLRKOZEJKmgW6fFuk4B6k3jLV3xuumshBJU6fbuwvc/0LqDc0ue/3qse7PzH9rTTmSpoLdBZJaodmBldeMcnwD8BuKbbgl9Qi7CyS1QlMhIjO3mQoaEbMotgn/Z+BjbahLUht1orug2wdzSpq4pqd41srMLcBdEXEq8BmKBagkaVTdPphT0sQ1u+z1aDZT7FMhSWNyjwtp+pnMwMqFFFuD39rSiiRNSw7mlKafyQys3AR8F/jb1pUjabpyMKc0/ZQaWCmpe3XrAEbXfpCmn1IDKyV1LwcwSpoqzY6J+H3gEuAAYLs2yMyc3eK6JJXU7atRSpo+mm2J+ATF1t/vANa3rxxJk+UARklTpdkQEcBzMvPOdhYjafIcwChpqjQbIn4A/A5giJC63GgDGLt1wKWk3tVsiDgJ+HxEPAf4JTBce6cbcEndzwGXklqt2RDxMuD3gLMa3DeCG3BJXc8Bl5JardkQ8TbgDOCCzFzXxnoktYkDLiW1WrMhYhbwaQOE1LsccCmp1ZpdifITwN+0sxBJ7VUdcLnfk3dhyaKdHFQpadKabYlYBBwfEX8J3EOxb8ZWmXlMqwuTJEndrdkQsSPw6XYWIkmSekuzG3C9rt2FSJKk3jJqiIiIVwOfzcyNle9HM5KZtlJIbeZiUZK6zVgtEdcA/wM8XPl+NCPY1SG1XSsWizKISGqlUUNEZs5s9L2kzmjFYlGPrRnireffwOAaV62UNHkTDgcRsUNEHBwRT2pHQZIaqy4WBUx4sagtwyMMDG7gNw+v5sy/PpSn/87OW4OIJJU15sDKiDgOOAV4WWb+OiKeAXwZ2AcYiYhPAm/MzC1tr1Tqc5NZLKq+K+SsEw/l7CtvcdVKSZMyaktERLwS+CTwE6D658qngIXAi4DDgEMplsSW1GaTWSyqvitkw8bNXHzaUletlDQpY3VnvBU4IzNPyMyVEfFs4CDg4sz8embeBrwXcPqn1AWqXRZ33buSgcENDA+PbL2vvivkSbstYOeFrlopaXLG6s54FnBizc9HU8zE+GLNsR8DT21DXZImaKzZG7VdIXvtMo/hkRHuunelMzQkTcpYIWImsLHm5yOAQeD7NceeAGxoQ12SJmis2RvVrpAli3ZiYHADJ5+3bFJTRSUJxu7OuBN4PkBELAT+L/C1zBypOefPKcZMSOqwZmdvNAobklTGWC0RlwIXR8SzKMLEE4ALACJid+DVwLuAN7a7SEnja3b2RjVsVFsinKEhqayxFpu6OiJ2Ak4CtgCvysxbKnefSTFe4sOZeXX7y5Q0nuMPfXAAABh4SURBVNoui7FMZqqoJNUac52IzLwCuKLBXf8EvC8zV7alKkltUxs2tgyPsMplsCWV1OxW4NvIzN+0uhBJU68V+3FI6l+lQoSk6aF+kOWq1Ru2HrdlQtJ4DBFSH6sdZLlo/mwWzp9jy4SkphkipD5WO8hyz13msXzl5HcKldQ/DBFSH6uf0TEygtM/JTXNECFpK6d/SpoIQ4SkrZpda0KSYOxlryVNwFi7aHbTNSWpVWyJkFqkHWsuuI6DpG5mS4TUIu3Y2MrNsiR1M0OE1CLN7qLZ6WtKUqvYnSG1SDtmNjhbQlI3M0RILdKOmQ3OlpDUzezOkCbJGRSS+pUtEdIkOYNCUr+yJUKqUaZVwRkUkvqVLRFSjTKtCrU7YTqDQlI/MURINRq1KowXIpxBIalfGSKkGmVaFVo9g2LL8AiDq4e2CSUzZ85oybUlqZUMEVKNVrcqlAkEDtSU1CsMEVKNVrcqjBYIxgoXZbpUJKkTDBFSG40WCMZqbXCgpqReYYiQ2mi0QDBWa4MDNSX1CkOE1EYL583mwlOP4sEVa9hr1/lbN9Maq7XBpa4l9QpDhNRGj63dyNs+uoyREXjq3gt56ysP5pFH17HHkrlc+o6lPLjC1gZJvcsQIbXRQwNrGVyzkZkz4LgXPYNTzl+2zTiI/Z68S6dLlKTSui5ERMTuwCXAMcBG4CrgPZm5eYzHPAzsVnf4vZl5TtsKlZpQ7baYMQNWDK531oWkaaXrQgTwOWAEOBJ4EvBJYDPwnkYnR8QeFAHiCODnNXetbmuVUhNqB0nutvPcbcZBLFn0BAYGN7iYlKSe1VUhIiKeB7wAeEpm3gvcHhHvAC6OiPdn5lCDhx1AETJuzcyNU1iuBIy9oFTtIMnh4REuPm0pD65cy5wdZ/HhT/0vy1eudTEpST2rq0IEcDjwq0qAqFoGLAAOBG5t8JgDgHsMEOqUVas38Nbzasc6HMWSRU/Y7rxqoFixah2nXXgj1Q1C7daQ1Ku6bSvwvYEH6o79tnK7zyiPOQDYHBFfiojlEfH9iDiubRWqrzXaKnz5ynXbjHV4cOW6Ma+x6+K5zJ9bTPV0MSlJvWxKWyIiYl/g3lHuHgKuATbUHszMTRExAoz2p9r+wC7AeynGTRwLXBURO2TmVa2oW6pqtNLkroufsM1Yh10Xb98KUcvFpCRNF1PdnfEAsN8o9w0DJwPb/IsaETsCM4C1ozxuKTA7M6sDKW+PiN8FTqWY2SG1TKOVJp+42zzOfsPzePjRdey+81zm7Dh2A5+LSUmaLqY0RGTmJuDu0e6PiPuBF9cdfmLltr6bo3rNIYpWjFp3AH9ZskxpVI1Wmlwwdw5btsDw8DBLFu7Ewnm2LEjqD902sPIm4EMRsU9m3l85tpRiuuaP6k+OiB0oukfOz8yP1tx1CHBnu4tV/2m0jLUtC5L6VbeFiO8CtwDXRcRbgD2AD1GEhI0AETEfmJ+ZyzNzc0R8ETgjIu4Bfgr8GXAc8JKOPANNa7XLWM+YARed6vRMSf2rq2ZnZOYI8FLgIeDbFGMargTeX3Pa6cCDNT+/HbgcuIii9eE44JWZ+bWpqFn9pbqM9WNrNzK4phgTIUn9qttaIsjM5RRBYrT7zwLOqvl5iGJWRsMVLaVWGmv3TUnqN10XIqRuVp2euWLVOnZdPNfpmZL6Wld1Z0jdrrLIJFuqy01OQKOFqiSpl9kSIU1Ao8Wmmh1YOZnHSlI3siVCmoBGi01NxWMlqRsZIqQJqA6shInvezGZx0pSN7I7Q5qAyex74Z4ZkqYbQ4Q0AZNZndKVLSVNN3ZnSJKkUgwR6gutnl7pdE1JsjtDfaLV0yudrilJtkSoT7R6eqXTNSXJEKE+0erplU7XlCS7M9QnWj290umakmSIUJ9o9fRKp2tKkt0Z0pichSFJo7MlQhpD7SyMA5++K6e86mAeeXTd1i6MmTNndLpESeoYQ4Q0huosjJkz4Phjn8Ep5y9zWqckVdidIY2hOgtjwbzZrBhc77ROSaphS4Q0htpZGLvtPJeF82ZvbYlwWqekfmeIkMZQOwtjeHjEaZ2SVMMQITXJaZ2StC3HREiSpFIMEeoY12CQpN5md4Y6xp0wJam32RKhjnEnTEnqbYYIdYw7YUpSb7M7Qx3jTpiS1NsMEeoYp0xKUm+zO0OSJJViiJDG4DRUSRqd3RnqWVuGRxhcPbTNmIpWb83tNFRJGp0hQj1rKj7gG01DNURIUsHuDPWkLcMjPLiy/etMOA1VkkZnS4S60nhdFYOrh5iz46y2b83tNFRJGp0hQl1pvK6KhwbW8on/vJOzTjyUR1at56l7L27LB7zTUCVpdHZnqCuNtyT2HkvmsXzlWk6/6Ftc/eW72HHWzJYPqpQkjc2WCJXS7pkR1bEIo3VV2M0gSZ1niFAprZoZMVoYGS8k2M0gSZ1niFAprZr6OFoYMSRIUvdzTIRKadXUx/owsnxgLY8+5gqRktQLbIlQKa0ak1A/9mH3nedy8nk3MLjGFSIlqdsZIlTKZLobasdB7LXrPC4+7SgeGljHHkvmsWrNBgbXuEKkJPUCQ4SmXKNxEPs9eZet97d7ASlJUmsYIjTlxhqU6dRNSeodhghNubHWgHBWhiT1DkOEppytDZI0PRgiNOVsbZCk6cF1IiRJUimGCEmSVIohQh2zZXiEgUFXp5SkXuWYCHVMqzbxkiR1hi0R6phG60VIknqHIUId06pNvCRJnWF3hjrG9SIkqbcZItQxrhchSb3N7gxJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpXTtOhERMQe4DfhIZl4zzrmvAd4H/A5wO3ByZn6v/VVKktS/urIlIiIWAJ8HntXEuUcD/wKcBxwM3AF8LSJ2a2uRkiT1ua4LEZVQ8CNgjyYf8g7g05l5RWbeBbwRGADe0KYSJUkSXRgigBdTtCwcNt6JETETeD6wrHosM4eBbwGHt6k+NbBleISBwQ387FcDDAxuYHh4pNMlSZLarOvGRGTmqdXvI2K80xcD84AH6o7/FnhOaytT1ZbhEQZXD22zcdbqNUMMPLaBFYMbmDlzJjNnwuIF7okhSdPZlIaIiNgXuHeUu4cyc6KfOnMrtxvqrwX4CdYmg6uHOPm8G3hs7UYWzpvNxactZfOWYc78+He3Hjv/bUd2ukxJUptNdUvEA8B+o9w3XOJ66yu39XtIzwHWlriemvDQwFoeW7sRgD13mcfm4eFtjj22diMrVq1njyVzx7qMJKnHTWmIyMxNwN0tvOQARVjYq+74E9m+i0MtsseSeSycN5t16zfylpc/m9MuvJH3vf5QFs6bvbUlYq9dDBCSNN113ZiIicjMkYi4GTgS+BRsHWx5BPDxTtY2nS1eMIdLT1/Kxs1FC8T7Xn8o/3XzvZx14qFs2LiZJ+22gMUL6huHJEnTTc+FiIiYD8zPzOWVQ+cDX4yIHwLfBE4FFgGf6FCJ097MmTMYHoG3X3Dj1paHs048lPf/yy1cdOpSdl7ocBRJ6gfdOMVzPKcDD1Z/yMyvACcBpwE/AJ4BHJOZKzpTXn+oHwOxYeNmLjp1KYvm2wIhSf2iq1siMnNGg2NnAWfVHbsKuGpqqhI8Pi6i2hLxpN0W2AIhSX2mq0OEutfiBXO4+LSl26wVIUnqL4YIlTJz5gyWLNqJJYtsfZCkftWLYyIkSVIXMERIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJK2aHTBYwmIuYAtwEfycxrxjn3YWC3usPvzcxz2lWfJEn9ritDREQsAK4DntXEuXtQBIgjgJ/X3LW6PdVJkiTowhAREUcDHwNWNfmQA4DNwK2ZubFthUmSpG1045iIFwP/AhzW5PkHAPcYICRJmlpd1xKRmadWv4+IZh5yALA5Ir4EHAI8AFyQmZ9qT4WSJAmmOERExL7AvaPcPZSZO5W47P7ALsB7gfcAxwJXRcQOmXlVk9eYBbB8+fISv16SpN5S83k3azLXmeqWiAeA/Ua5b7jkNZcCszOzOpDy9oj4XeBUoNkQsRfAa17zmpIlSJLUk/YC7in74CkNEZm5Cbi7xdccAobqDt8B/OUELvM94HDgQWBLi0qTJKlbzaIIEN+bzEW6bkzERETEDhTdI+dn5kdr7joEuLPZ61SCyE0tLk+SpG5WugWiqudCRETMB+Zn5vLM3BwRXwTOiIh7gJ8CfwYcB7ykk3VKkjTd9VyIAE4HzgRmVH5+O/AocBFF08zdwCsz82udKU+SpP4wY2RkpNM1SJKkHtSNi01JkqQeYIiQJEmlGCIkSVIphghJklSKIUKSJJXSi1M8WyYi5gC3AR/JzGvGOfdhYLe6w+/NzHPaVV+3mODr9BrgfcDvALcDJ2fmpFZE6wURsTtwCXAMsJFiyfX3ZObmMR4z7d9TETELOAc4AVgAfAV4c2Y+NMr5hwAXAgdRLJP/D5l59dRU21klXqvrgZfXHf5GZh7dzjq7SUR8DJiVmSeOcU7fvqeqmnydSr2f+rYlIiIWAJ8HntXEuXtQ/GN/BMVaFNWvj471uOlggq/T0RTbuJ8HHEyx/PjXIqL+g3I6+hywJ3AkxYfA64CzRzu5j95TZwGvBY6neK57U7xW26m8T74K/IDi/XMRcGVEHDMllXbeWTT5WlUcALyLbd8/r2hvid0hImZExPuBk8Y5r6/fU82+ThWl3k992RJR+bD7GLCqyYccAGwGbs3MjW0rrMuUeJ3eAXw6M6+oPP6NwAuBNwAfaEuRXSAinge8AHhKZt5LsQncO4CLI+L9lWXV603791REzAZOAd6amV+vHPsL4N6IOCwzb657yInAIHBKZg4Dd0fEwRQLzE3rxeMm+lpVzn8acFtm9tX2wxHxFOBKiv+Hfj3O6f38nmr6dZrM+6lfWyJeTPEX82FNnn8AcM90/cd+DE2/ThExE3g+sKx6rPI/7bcoNjebzg4HflUJEFXLKJqkDxzlMf3wnjqQ4jVYVj2QmfcB99H4PXE48K3K+6ZqGfD8yvtrOpvoa7UfxR+Bd7W/tK7zPOCXwDMp9k4aSz+/pybyOpV+P/VlS0Rmnlr9PiKaecgBwOaI+BLF5l4PABdk5qfaU2F3mODrtBiYR/Ha1Pot8JzWVtZ19qbx8wbYB7i1wWP64T21d+W20Wuzzyjn/7DBuXOBJcCKllbXXSb6Wh1AMfbm7Ig4FlgPXA+ck5kb2lZlF8jMa4Froal/l/r2PTXB16n0+2nahYiI2JfRU9dQZu5U4rL7A7sA7wXeAxwLXBURO2TmVaUK7bA2vE5zK7f1b7ghoMxr3jXGe62Aa6h73pm5KSJGGP25T7v3VANzgeHM3FR3fLT3xFwav38Y5fzpZKKv1f4U+wclxYDeZwLnUwSO17axzl7Tz++piSj9fpp2IYIiye83yn3Doxwfz1Jgdmaurvx8e0T8LnAqxSj8XtTq12l95XZO3fE5wNoS1+sm471WJ1P3vCNiR4r/KUd77tPxPVVvPTCzEoxqZ6mM9p5YT+P3D6OcP51M9LU6Azg3MwcqP98REVuAf4+IUzNzZZvr7RX9/J6aiNLvp2kXIipJ/u4WX3OIx9Nr1R3AX7by90ylNrxOAxT/U+5Vd/yJbN9E21PGe60i4n6K8SO1nli5bfjcp+N7qoH7K7d71XwPo78n7qfx+2cNxeC46WxCr1Wlj3+g7vAdldt9AENEoZ/fU02bzPtpug8smbSI2CEi7o+It9fddQhwZydq6kaZOQLcTDHFEdg62PIIisGV09lNwFMiorbveimwGvhR/cl99J66neI1qH1P7AvsS+P3xE3AERExo+bYUuA7dQPjpqMJvVYR8ZmI+Hzd4UMogukv2lZl7+nn91TTJvN+mnYtEa0QEfOB+Zm5PDM3R8QXgTMi4h7gp8CfAccBL+lknZ1W+zpVDp0PfDEifgh8k6JpfhHwiQ6VOFW+C9wCXBcRbwH2AD4EnF+dfdGP76nMHIqIy4BzI2IF8DBwGXBjZt5SmVa2BBiovE5XAu8ELo+IC4CjgVcDL+rMM5g6JV6rz1Jpaga+QLGQ0rkUTdJrOvMsOs/3VHNa+X6yJaKx04EHa35+O3A5xUIld1L8Y//KzJzW84ybsM3rlJlfoVjU5DSKxV2eARyTmdN2BDRsbYV5KfAQ8G2KMQ1XAu+vOa1f31NnUIwQvwa4AfgVj6+KdxjFa3IYQGVlxhdR/AP2Q+AtwPGZ+c0prrlTJvJafYbHFzX7CcUCbxdSrBbbz3xPNadl76cZIyMjbatSkiRNX7ZESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFFSulPhERx1MstrM/xcZhPwYuyszras4ZAY7LzGvaVMMngb0z8+gmz38G8OTM/K9J/M5PAE/LzKMa3PdUiiWnr8nMv6m7748oVu/7/zLzf8r+fmk6syVC6gMRcRLFFr+XAc8G/gD4L+DTEVG71e9eFEvgdosvAM9p18Uz8x6K5dnfGBFbN1Gr7IPyr8CHDRDS6GyJkPrDG4GPZ+Yna479NCICOIXiA5OafVC6xYzxT5mczLwiIv4YuDIingmsAv4d+Dnw3nb/fqmXGSKk/rAFeH5ELMrM2i2QTwfmVX+o7c6odD1sotgy+cTKNS4A/h9wBXAwkMCJmfm/9Y9vdM36oiLiz4F3AQcAIxT7G7wtM78XEcuApwJnRsQJmblvRMwBPkCxidK8yvl/l5m31FzzzZXntUel1mZaXE+k2Pr4YuCuSj0HZubmJh4r9S27M6T+8BHgucBvI+I/I+L0iDgwMx/JzPvGeNzxldvfBz5KsanYfwD/VLneRuDSMgVFxHOAzwCfBPaj2AZ7BvDxyikvA+6j2Ayo2qVxNcX28q+k2Kr4m8ANEfH0yjWPo9hN9gMUmy79BvjL8WqpbNR0EvAXFK0Pb8zMe8s8L6mfGCKkPpCZ1wMvoBgHcQRFqPhhRPwgIvYf46GPAO+ojB34aOXYv2XmlzLzDoodSw8oWdYm4G8z89LMvC8zv0cRIJ5ZqXmAovVjTWY+EhFPowgPJ2TmtzPzZ5l5NnATxc6xUAwcvSYzP56FdwHfa7KeG4FHgc3AzSWfk9RX7M6Q+kRm3gzcHBGzKFoW/hg4GfjviHhaZm5s8LB7Kludk5lriyEU3FNz/3pgTsl6fhQRqyLi3RTbxv8ecCCj/3FzUOX21kodVXNqajiAyviOGrcAz2qipI8DKyn+Xbw6Il6YmcNNPE7qW4YIaZqrzDR4N/D+zFyemVuA24DbIuLbwFcpPmT/t8HDNzU41vQHa0SM+m9MRCwF/ptiBsZ3gH8Bng5cPspDqiHneRThpdZQ5XaE7QdjNgpH9bW8EXgp8EKKfxe/TtG68ZHxHiv1M7szpOlvPcXAwVc3uG8VxQfvwy36XZuAhTU//94Y5/4t8LXMfFVmXpSZNwD7AkRENQiM1Jx/Z+V2j8z8RfULeDvwp5X7fgQcVvd7Dhmr4IjYj6Kr5tzMvDEzv0ExwPKciHj2WI+V+p0tEdI0l5krIuLDwD9FxELgcxTB4pnAOcC/ZuavW/TrvgucFBHfAWZRfDgPjXLuI8BLIuJQ4CHgj4C3Ve6bA2wAVgNPj4gnZuYvIuI64IrKDIyfAa8H/gY4pvK4c4HPRsRtFK0cr6QYC3JTowIiYieK6ZzJttM531W55jURcUhmjvYcpL5mS4TUBzLzDIq1Iv6Q4gP1TuAfKcYPnNTCX/Umiimht1IsWnUFxQyJRt5HMUXzq8D3KWZjVBe+qs7GOB84FvhxRMykaFH5MsWAzp9U7ntZpfWAzPwP4ATgzRQrch5aqWE0H6HoQnlN7ZiQzFwPHAf8H4qZKJIamDEyMjL+WZIkSXVsiZAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSV8v8D5+zlv4QXwTsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x, y = sims.linear_sim(num_samp=100, num_dim=1, noise=0.1)\n", "\n", "fig = plt.figure(figsize=(8,8))\n", "fig.suptitle(\"Linear Simulation\", fontsize=17)\n", "ax = sns.scatterplot(x=x[:,0], y=y[:,0])\n", "ax.set_xlabel('Simulated X', fontsize=15)\n", "ax.set_ylabel('Simulated Y', fontsize=15) \n", "plt.axis('equal')\n", "plt.xticks(fontsize=15)\n", "plt.yticks(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The test statistic and p-value can be called by creating the `DCorr` object and simply calling the corresponding test statistic and p-value methods. The value of `which_test` can be modified to 3 different values: `biased` for biased Dcorr, `unbiased` for unbiased Dcorr, and `mantel` for Mantel." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dcorr test statistic: 0.9724847196117954\n", "P Value: 0.001\n" ] } ], "source": [ "dcorr = DCorr(which_test='biased')\n", "dcorr_statistic, independence_test_metadata = dcorr.test_statistic(x, y)\n", "p_value, _ = dcorr.p_value(x, y)\n", "\n", "print(\"Dcorr test statistic:\", dcorr_statistic)\n", "print(\"P Value:\", p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Due to the similarities between Dcorr and MGC, many of the properties of the values of the test statistics are similar (value of 1 for linear, etc.). The p-value is bounded by the number of repetitions (in this case 1000). This is because since we are estimating the null distribution via permutation, this is the lowest value that we can be sufficiently sure is the p-value. It is worth noting that as in most of the other tests that use permutation to approximate the p-value, the `replication_factor` parameter can be set to the desired number.\n", "\n", "Unlike in the case of MGC, Dcorr struggles to consistently identify some nonlinear relationships." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAIlCAYAAADbvPFqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZiddX3//2cmmIQEQghhs1jxa+UtbmCrFrCKWEurtdpalxZc0B+CSEUWF6oW4atVaxURxLrXBbVatS1+616WuuG+G99WhKoIgSQwZCEJZOb3x30PnhxmJmcmZ7nP+Twf15XrzNxnzjnvOTOZ1/1Z7wWTk5NIkqRyjA26AEmS1F+GvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EsaGRGxYNA1SMNgt0EXIA2riLgv8DLg0cB+wK3AN4CLMvNT83i+SeAfMvPsLtR2LXBVZv7VTr7uOOAk4HBgCXAd8FngtZn5y7k+XzdExKOAy4HHZuZn5vC4o4GzgcfWnx8MXAOckplv636l0vCy5S/NQ0TcH/g6cDDwYuCPgJOBbcB/RsTz5vG0RwIXd6vGnYmIc4D3A98CjqcKzTfWt9+MiN9u+fK/AF7Rr9rm6RTg0JbPr6d6Tz8xmHKk5rLlL83PmcAm4DGZuW3qYER8AvgC8JqIeEdmTnT6hJl5VffLnF5ELKJqJV+YmWe13HV5RHwG+CnV93h6Xdt3+lVbt2TmVqBv76k0TAx/aX4OqG8Xth7MzMm6RX0ksAjYEhHnAs8DngmcD9ybKlxflZkfm3psa7d/S5f1mcBzgHsAr87MN0TEw4GXA78PLAfWAP8GvDQzN3dY/17A7kzT+5eZP4+IF1ANAUzVdi0t3f51racBDwSeUr8P/wY8v/73AmBv4MvASZn5i/pxVwBLMvOIluee+l5n7J7f2fdcP+/RLbU9G7ii/XkjYjlVD8ZfUL2n1wJvBy7IzMmWGq8Fvlt/j78FrAZeNp/hHKmJDH9pfj4JPA74ekS8B7gM+EFmTmTml6lCr9VewEeA1wI/oBpn/2hEPCkz/32W13kdcBbwC+CnEfFAqvHwT1J11d9R13EGsA44t5PiM/OmiPgacFpEHEDVNf7FzLy+vv+fOniaV1MF8F8CfwCcBzwMuJGqC/4g4M3AW4HHd1LXdDr8np9fv9YDqIL9amBZ2/MsAb5IFfqvBBL4Y+ANwCF1zVP+nOrE5iXA5vp7/URE3CMzb5rv9yI1heEvzUNmvi0i9gX+lqo1D3Br3Wp8d2Ze2vaQJcDpmfl2gLpr/btUwTVb+H8yM98y9UlEPJ2qRfu0zLyjPvyFiHgMcAwdhn/tScB7gafV/4iInwP/SdUS/vlOHn9NZp5Qf3xZRDybKvCPzMxb6ud7WP06u+IwdvI9Z+aPI2IdsHVq+CQilrU9zwnAg6iGav6rPva5iNgGnB0RF2Xmj+vjuwN/nJlr6+faCFwJHAt8cBe/H2ngnPAnzVNmvoqq+/+pVK3bXwBPAP4jIj7Ytuxskmpy3dRjJ4GPA4dFxIpZXuYHba95SWYeCyyMiAdExBMj4hXA/lTDDHOp/9f1cwXwwrqePam67H8cEX+2k6f4Stvna4CfTAV/bS0w2/fXSZ3d+p6PAda0BP+UqZ/Lo3Z82Sr4a7+qb9tPKKShZMtf2gWZeSvwr/W/qfHri4HjqFqIU2PE6zLztraH31jf7g3cwvTWtH4SEYupurefRdWb8AuqVQebgXmtcc/Mn1LNQbiwPmH5E+B9wDsi4rdmmbR46zTHNrV9vsuXDe3i97wSuGGa49fXt60nKe1zJ6beAxtMGgn+IktzFBG/FRHX1ZPidpCZ1wIn1p/ev+WuvSNiYduX708VjnMZQ76Aatz7WcDyzLxnZj5ljs9BRLwwItZFxMq2+icz89PAm6h6NfaZy/N2YJK2SZJUvQ2z6cr3DKznNxM1W929vl07zX3SSDL8pbm7AbgdOCUi9pjm/kPq2++3HFtINSQA3LkT3ZOBr2bmxjm89tHAlzLzo5m5oX6ug6gmp83l//MPqVrCZ8xw/yFULeJuB+KtwG+1DYk8cieP6fR73r6T57kc2D8i/rDt+DPq2yt38nhpZNjtL81RZm6PiJOBS4HvRMRbqIJ+DDiKKlAvzczPtj30nRGxP/C/VEv/7gs8Zo4vfxVwXET8DdV8gKDaZXAxcxiPzsz/ioj3Aa+oNyz6KNXSvlVUS/eeRjXBbpe77dtcSnUSdHFEfJRqZ8EzmT24O/2eb6YK98dSTaZs9z7gVKpVFucCP6HanOks4F2ZmbvwfUlDxfCX5iEzPxsRD6Ha3e90qu7k7VSBch7T79T3XOD1VOvGvwMcm5lzbW2eRdWLcA6wlGr8+z31a58bEQdk5nTj2tN5NtUSxWcBF1EtR7yFaiLfI+sli932XuBeVHsXPIdq7P6JwFdneUyn3/M7qU6m/qP+2n9pfZLMvK3eAvg1VHsG7A38rH7+C7vz7UnDYcHkZLdP7CW1qluZrwR2z8wtAy5HkhzzlySpNIa/JEmFsdtfkqTC2PKXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwuw26AJmEhEnAi8B7gH8GHhxZl422KokSRp+jWz5R8SzgIuB1wEPBK4ELo2IgwdZlyRJo6Bx4R8RC4DzgH/IzPdk5s+AFwE/A44aaHGSJI2AJnb7B3BP4CNTBzJzAjh8YBVJkjRCmhj+h9S3KyLiMuABwE+AszPzK4MrS5Kk0dDE8F9e374POIcq+E8ELouIB2fm6tkeHBGLgYcC1wPbe1moJEkNsRA4EPhGZm7d2Rc3Mfxvr2//PjM/BBARpwKPAE4BTtvJ4x8KfLF35UmS1FiPAL60sy9qYvhfV9/+YOpAZk5GxGrgXh08/nqAD37wgxxwwAE9KE+SpGa54YYbOP7446HOwJ1pYvh/G9hE1YL/Jty5AuB+wBc6ePx2gAMOOICDDjqoVzVKktREHQ13Ny78M3NzRLwJ+PuIWEPVA/B84N7AXw60OEmSRkDjwr92DrAZuADYD/gucGxm5kCrkiRpBDQy/DNzEnht/U+SJHVR43b4kyRJvWX4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwhj+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwhj+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S0Nk+8Qk68e3sPqadawf38LExOSgS5I0hHYbdAGSOje+YSsveOPl3LppG8uXLeKis45h5V5LBl2WpCFjy18aImvWb+LWTdsAuHXTNtas3zTgigbP3hBp7mz5S0Nk/5XLWL5s0Z0t//1XLht0SQNnb4g0d4a/NERW7LmYi846hjXrN7H/ymWs2HPxoEsauOl6Qwx/aXaGvzRExsYWsHKvJYZbi270hmyfmGR8w9YdTqrGxhb0oFqpGQx/SUOtG70hDh2oNIa/pKHWjd4Qhw5UGmf7SxpJc1kFMDV0ADiRUkWw5S9pJM2lK9+JlCqN4S9pJM2lK9+JlCqN3f6SRpJd+dLMbPlLGkl25UszM/ylLnPNeDPYlS/NrNHhHxFHAF8CHpOZVwy4HKkjrhmX1HSNHfOPiGXAB4CFg65FmgsvviOp6Rob/sD5wK8GXYQ0V040k9R0jez2j4jHAX8KPBb4/oDLkebEiWaSmq5x4R8Rq4B3Ac8Bbh5wOdKcOdFMUtM1sdv/7cAnM/Mzgy5EkqRR1KiWf0Q8C3gw8KBB1yJJ0qhqWsv/BOAg4IaI2AhkffzTEfG2gVUlSdIIaVTLH3g6sHvL5wcAXwROBD4/kIokSRoxjQr/zLyu9fOI2FJ/eF1m3jiAkiRJGjlN6/aXJEk91qiWf7vM/BXgpujqCvfcl6RKo8Nf6ib33Jekit3+KoZ77ktSxfBXMdxzX5IqdvurGO65L0kVw1/FcM99SarY7a++2z4xyfrxLay+Zh3rx7cwMTE56JIkqSi2/NV3zrqXpMGy5a++c9a9JA2W4a++c9a91B0OoWm+7PZX3znrXuoOh9A0X4a/+s5Z91J3TDeE5v8rdcJuf0kaUg6hab5s+UvSkHIITfNl+EvSkHIITfNlt78kSYUx/NU1LjuSpOFgt7+6xmVHkjQcbPmra9y5T5KGg+GvrnHZkdRfDrVpvuz2V9e47EjqL4faNF+Gv7rGZUdSf7nDn+bLbn9JGlIOtWm+bPlL0pByqE3zZfhL0pByqE3zZbe/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDv3DbJyZZP76F1desY/34FiYmJgddkiSpx7ykb+HGN2zlBW+8nFs3bWP5skVcdNYxXh5UkkacLf/CrVm/iVs3bQPg1k3bWLN+04ArkiT1muFfuP1XLmP5skUALF+2iP1XLhtwRZKkXrPbv3Ar9lzMRWcdw5r1m9h/5TJW7Ll40CVJknrM8C/c2NgCVu61xHF+SSqI3f6SJBXG8JckqTCGvyRJhTH8JUkqjOEvaai4K6W065ztLw3Y9olJxjds3WG55djYgkGX1VjuSintOsNfGjDDbG6m25XS90uaG7v9pQFzi+W5cVdKadfZ8pcGbCrMplr+oxhm3RzacFdKadcZ/tKAlRBm3RzacFdKadcZ/tKAlRBmjtNLzeKYv6Sec5xeahZb/pJ6roShDWmYGP6Seq6EoQ1pmNjtL0lSYQx/qQHcslZSP9ntLzWAu/xJ6idb/lIDuMufpH4y/KUGcCmcpH5qXLd/ROwPvB44Ftgd+BpwVmb+cKCFST3kUjhJ/dSoln9EjAH/BhwCPBE4ChgH/isi9hlkbVIvTS2FO/Re+7ByryVe0ldSTzWt5X8YcCRwv8xcDRARzwDWA38KvH+AtUmSNBIa1fIHfgE8HsiWYxPAAmDvgVQkSdKIaVTLPzPXAf/Zdvg0YAnwuf5XJGlUdfMyw9o53+9maVT4t4uIJwCvBc6fGgaQpG5wb4X+8v1ulqZ1+98pIk4APg58BHjJYKuRNGrcW6G/fL+bpZHhHxEvB/4ZeBvwzMycGHBJkkaMeyv0l+93szSu2z8iXgK8GjgnM1816HokjSb3Vugv3+9maVT4R8SDgNcA7wHeGREHtNy9ITPtJ5LUFV5muL98v5ulad3+fwUsBJ4DXN/274wB1iVJ0shoVMs/M18GvGzQdUitXKJUBn/OKkmjwl9qIpcolcGfs0rStG5/qXFcolQGf84qieEv7YRLlMrgz1klsdtf2gmXKJXBn7NKYvhLO+ESpTL4c1ZJ7PaXJKkwhr9U2z4xyfrxLay+Zh3rx7cwMTE56JIkqSfs9pdqLvWSVApb/lLNpV6SSmH4SzWXekkqhd3+Us2lXpJKYfhLNZd6SSqF3f6SNEeuDNGws+Wv4nj1Nu0qV4Zo2Bn+Ko5/uLWrplsZ4u+Qhond/iqOS/q0q1wZomFny18ja6bu/ak/3FMtf/9wa65cGaJhZ/hrZM3Uve8fbu0qV4Zo2Bn+DeEktO6baVzWP9ySSmf4N4ST0LrP7n1Jmp7h3xDOHu4+u/claXqGf0PYSu0+u/claXqGf0PYSpUk9Yvh3xC2UqXR4ORdDQPDX5K6yMm7Ggbu8CdJXeQOkhoGhr8kdZFb/2oY2O0vSV3k5F0NA8NfkrrIybsaBoa/hlLrjOr9Vi5l4dgCrl/r7GpJ6sSM4R8Rp2bmxf0sRupU+4zq8557JGdf/CX2WOrsaknamdkm/L05Ir4QEffoWzVSh9pnVN9482b2WLrI2dWS1IHZwv8PgP2BH0bE/9eneqSOtM+o3m/vpWzc7NbIktSJGbv9M/OqiHgwcDbwloj4S+C5mXld36qTZrDjjOqljI0t4HWn/oGzqyWpA7NO+MvMO4BXR8SHgTcDP4iItwKb277uNb0rUbqr6WZUr9jTcX5J6kSns/1/BXwDOBZ4DrC15b5JwPCXJHWN10jorZ2Gf0Q8DrgY2Ac4wxUAkqReu2XDFk574xUt10h4FCv32n3QZY2M2Zb67QdcCDwF+AJwUmb+b78KkySV64Z1m3dY0XP9us2GfxfN1vL/SX17Ymb+cz+KkSQJYNWK3Vm+bNGdLf9VKwz+bpot/K8ETsnMG/pVjCRJAIvvNsZ5zz2SG2/ezH57L2Xx3bwOXTfNttTvL/pZiKT+mW4y1SQ4wUqNsXzZYiYmYAGT7L18CcuXuYS3m9zbXypQ+/bIF511DMBdjrlNsgZlsr7ddsfEQOsYVYa/VKD27ZHX3rKZ7ROTOxxbs36T4a+Bme4E1d/H7nEQpY+2T0yyfnwLq69Zx/rxLUxMTO78QVIPtG+PvGrF0rscc5tkDVL7CarX7OguW/595JmsmmLH7ZF/syXy1LH9Vi5lbAxWX7PO8X8NxNTJ6NTfS09Gu2u2df6Xdfokmfno7pQz2qY7kzX8NQjTbY8MsNeei1mwAH69dhOL77aQd136I25Yt2nGE1V3YVOvzHSCqu6YreV/dcvHS4G/ptri9ypgG/AQqiv/vb1n1Y2YqTPZjZu3cfd99/BMVo3T3jt17olH8KIL/3vGE1V7s9QrM52gqjtmW+r33KmPI+IDwD9k5t+2fk1EvAJ4WO/KGy0r9lzMxS86hq13TLD2ltuASSYmJm0pqSfm0ypv75266ZbbOHDVzCeq9mZJw6nTMf+/AB48zfGPAC/rXjmjbWxsAROTcOYFV9pSUs/Np1XePs5674NW8JpTHj5jl6vjstJw6jT8bwQeAfxP2/HHAr/sakUjzpaS+mU+v2vTjbPO1lvguKw0nDoN/zcCb42II4BvAwuAo4CnASf0prTRZEtJ/TKf37W5jrM6LisNp47CPzMvjohbgVOB46g2X/ou8JTM/I8e1jdybCmpX3bld81Z/NJo63idf2Z+APhAD2spgi0l9ct0v2udhrqz+KXR1nH4R8TvAGcB9wWeDjwR+HFmXtGb0iR1W6eh7twUabR1tL1vRPw+8D3g3lRj/YuB+wOfj4jH9648lchtkOdmLu9Xp1umutWvNNo6bfm/DnhdZr4qIjYAZOapEXELcC7w/3pUnwpkl/PczOX96nQSoHNTpNHWafj/LvDcaY6/Gzi9e+VIdjnP1Vzer05D3bkp0mjr9Kp+m4D9pjl+CDDevXIku5znai7v11SoH3qvfVi51xJn8EuF6rTl/2HgTRFxAtUyvyUR8WjgLcDHelSbCmWX89w0+f3aPjHJho1b79zS+sB9lrJiT086XErZGd+n3uk0/P+Wapnfj+rPf1DffrS+T+oau5znpsnv1/iGray/dQuvfOdXncPRwnktnfF96p1ON/nZBjwtIl4OHE51Vb8fZebVsz9SUq8MQ6to7S2bWTu+xTkcbZzX0hnfp97pKPwj4ufAQzLzZ8DPWo4fCHwvM6ebDyCph4ahVbRqxVLGxsbc0rqN23x3xvepd2YM/4h4HPCQ+tODgbMjYmPblx0y23PMV0QsBF5Ndd2APYHPAKdm5ppuv5Y0rIahVVT1RsD5px/dMubfnDkJg9LkeRpN4vvUO7MF9zXABVQX8QF4MrC95f5JYAPwgh7UdS7wLOCZwDrgrcDHgT/owWtJQ2f7xCT77r20ka2imYYj9l+5dNClNUaT52k0ie9T78wY/pm5mqplT0RcDjwpM2/udUERsQh4IXBaZn6+PvZXwDURcVRmfqXXNUhNd8uGLbz5I9/m3BOPYO34bdznHns3plXUPhxx8YuOYWKSRs9NkErT6YS/Y6Y7Xgf1QzPzy12s6XCqrv4rWl7/2oi4FngEYPiPqGGYwNYUN6zbzHd/upbv/89/s8fSRbzshIexasXugy4L2HE4YuPmbWy9Y4IzL7iy0XMTpNJ0OuHv94B3Ag9k+o2BFnaxpoPq2+vajv8auEcXX0cNMwwT2Jpi1Yrd7+zyn/q8KVonad193z1Ye8ttjZ+bIJWm08l6bwZuA06iGn9/IXCv+vaZXa5pKTCRmbe3Hd8K+BdjhN0wBBPYBm2qd2S33RZw3nOP5MabN7Pf3ktZfLdON+vsvfZJWjDZyLkJUsk6Df8HA4/MzG9FxElAZuY7IuLXwCl0d5e/24CxiNgtM+9oOb6Yapthjaj92iaw7bu3E8TabdhYbZqzbnwL97r7chYuXMBeyxazfFkzxvvhrpO0JiYmnbEtNUyn4b8AuKn++H+ouv+vBD4JvLLLNf2yvj2w5WOAu3PXoQCNkFs3beXcE4/gpltuY98Vu3Prpq2N6s5ugq13TOywW96bTj+68b0jztiWmqfTvsIfAo+rP/4x8PD64/3p7ng/wPeolhAePXUgIg6m2mvgv7v8WmqQFXss4bx3X8XFH/se5737KlbsYVi0ax8/v+mW2wZc0fS2T0yyfnwLq69Zx/rxLUxMTA66JEktOm35/wPwkYjYTnWRn1dGxL8DhwGXd7OgzNwaEW8F3hARa4EbqeYZXJmZV3XztdQsK/ZczIVn2j08mwP32XFo5MB9mjk04uRNqdk6Xer38Yg4Arg9M/+33v3vFOBTwDk9qOsVwN2AS+rbzwCn9uB11CB2D+/cij2XDMX4+TDsPiiVrOOteTPzmy0fX06XW/xtr3UHcFb9T1JtWE6Q3JNdarbZ9vb/XKdPkpnHdqcclcxNfkaHe7JLzTZby9+Z9eorx4lHx7D0UEilmm1v/2f3sxDJcWJJ6o9Ot/c9brb7M/ND3SlHJXOceOemhkbW3rKZVSuWOjQiaV46nfB3yQzHtwC/Agx/7TLHiXduaoe/teNbGBsbY2ysWgGg4eC8FjVFp0v9dtgMKCIWUl3u95+At/egLhXIceKda9/h7/zTj975g9QY4xu28vK3fZlf37SRPZY6r0WDM6+rgWTm9sxcDZwJvKq7JUmaSfsOf2sbusOf7mr7xCS3b5/gGY89lDec9kgO2GcZa9Z7uRINRsfr/GdwB9We+5L6YFh2+NNdjW/YypkXXHnnz+685x7JyuW2+jUYuzLhbznVJX6/1tWKVCTHQjszLDv86a7aV7NsvX27Pz8NzK5M+Lsd+Crw/O6Vo1K5xr8zzosYXu2rWQ7cZ5knuBqYeU34k7rNNf4ada5mUZPs6pi/1BWu8deos9dGTdLpmP/vAW8BHgDc5XQ1Mxd1uS4VxlaRJPVPpy3/dwHbgBcDri1S19kqkqT+6TT8A3hoZv6ol8VIkqTe63Qi37eB3+5lIZIkqT86bfmfBPxbRDwU+Dkw0XqnF/aR+st9ESTtik7D/0nAfYBzp7lvEi/sI/WV+yJI2hWdhv/pwCuACzJzcw/rkdQB90WQtCs6HfNfCHzY4JeaYWpfBMB9ESTN2VyW+j0PeGkPa5HUIfdFkLQrOg3/vYBnRsRfA1dT7et/p8w8ttuFSZqZ+yJI2hWdhv/dgA/3shBJktQfnV7Y59m9LkSSJPXHjOEfEccBH8vMbfXHM5nMTHsFJBXH/RY0rGZr+V8CfAG4sf54JpM4JCCpQO63oGE1Y/hn5th0H0tSybZPTHLrxq3cdPNmbt8+yYIFMLbA/RY0XDqd8HeniNgNeBCwJjOv635JktRM2ycmufnWLfzPL2/m7vssY8Xyu3Hqkw9j1V678/5P/9j9FjQ0Zg3/iHgG8ELgSZn5i4i4H/Ap4B7AZES8Fzg5M7f3vFJJGrDxDVt54flXsHHzNt74wqN52ZuuvLPL/81nPsr9FjQ0ZuzOj4inAu8Ffghsqg9/AFgO/AlwFHAE1da/kjTyprZV3nPZIm68efMOWyzfdPNmJ/tpaMw2ln8a8IrMPCEz10XEYcCDgYsy8/OZ+XXg7wCXAUoNtH1ikvXjW1h9zTrWj29hYmJy0CUNvaltlTds2sZ+ey91i2UNrdm6/R8EnNjy+WOoZvZ/suXY94F796AuSbvImejd17qt8qoVS7jorEexZv1mt1ieh9vvmGB841auX7uRA1ftwYo9FrPbbs4t75fZwn8M2Nby+SOBceBbLcd2B7b0oC5Ju8gr/3XfdNsqr9xr9wFWNLzGN1bzJ1rnTKxa4XvZL7OdZv0IeDhARCwH/hD4XGa29h3+JdWcAEkN45X/1GTXr924w8np9Ws3DriisszW8r8YuCgiHkR1ErA7cAFAROwHHAecDZzc6yIlzZ1X/lOTHbhqD5YvW3Rny//AVXsMuqSizLbJz/sjYglwErAdeFpmXlXf/Uqq+QCvz8z3975MSXPllf/UZCv2WMybz3zUDmP+6p9Z1/ln5juAd0xz12uBczJzXU+qktRT7kmvQdtttzFWrdjdcf4BmfMOfwCZ+atuFyKpf0pfCTC1Re8N6zz5UZlcVyEVaLqVAKWY2qL3J/+7nt0WjvGmf/kWt2zYOuiypL6aV8tf0nCbWgkw1fIvaSXA1Ba9U9/7uSce4TJIFcfwlwrQPsY/1dVf4kqA9l6PteO3cd97rhxwVVJ/Gf5SAWYa4y+xtdve63Gfe+zNXs40V2EMf6kA7vb3G9Ptf+BkP5XG8FfjuSxt15U8xt/O/Q8kw19DoPRlad3gbn+SWhn+ajy7rHedrV1JrVznr8bzAjWS1F22/NV4dlkPD+dnSMPB8Ffj2WU9PJyfIQ0Hu/0ldU2n2wZvn5hk/fgWVl+zjvXjW7jjjokdPp+YmOxn2eqS9p+rP8fmsuUvqWs6XVLY3kPw5jMfxelvuoLxjfYYDDN7foaH4a+h4phys3U6P6O9h+D6tRuZrBuJrugYXq7MGR6Gv4aKLYtm63R+RnsPwYGr9mBBfYSpcAQAABMnSURBVA7nio7h5WZSw8Pw11CxZTEa2nsIli9bxIVnuqJj2LkyZ3gY/hoqrS2Lww9Zxb57L2X1NescAhgyrT0ErUM5B+yzjL328Oc4rFyZMzwMfw2V1pbFvnsv3eG67A4BDCeHcqT+c6mfhspUy+LQe+3DTTdv7mhZmQank6VfnS4PlNQ9tvw1tJxc1HydtOr9OUr9Z/hraDm5qPk6maDpz1HqP8NfQ8vJRc3XSaven6PUf4a/pJ6xVS81k+EvqWds1UvN5Gx/SVLHvHjPaGhcyz8ifhd4PfAQYDPwKeAlmbl+oIVJktyXYUQ0quUfEXcHvgBcAxwJPAV4GPDRQdYlSaq4L8NoaFT4A08DtgDPy8zVmfll4FTgDyPitwdbmiRpagUHeBGmYda0bv9LgW9m5vaWYxP17d7AL/pfkiRpiis4RkOjwj8zrwaubjv8UuA64If9r0iS1MoVHKOhr+EfEQdTjedPZ2tm7vDbFBGvAx4P/Hlbb4AkSZqnfrf8rwMOneG+qe59ImIh8BbgZOCUzLy0D7VJklSEvoZ/Zt4O/GS2r4mIJVSz+/8EeHpmfqgftUmSVIpGjflHxBjwr8CjgT/LzM8OuCRJkkZOo8IfOIVqjP9E4HsRcUDLfevqngNpzrZPTDK+YesOM5THxhYMuixJGoimhf/x9e27prnvEcCX+liLRoi7kknSbzQq/DPzqEHXoNHUyXXlJdlLVopGhb/UK51cV16SvWSlMPxVBHclkzpjL1kZmra3v9QTUxcdXWj3pTQr9+4vgy1/FcGuTKkz9pKVwfBXEdq7Mq9ft8mJTNI03Lu/DHb7qwjtXZmL77aQWzZsHXBVkjQYtvxVhBV7Lub804/m6l/dwr4rduefPvF9TnzC/W3dSCqS4a8ijI0t4G4Lx3j/p1Zz/dqN7LHUiUySymX4qxgr9lzMa055uBOZJBXP8FcxnMgkSRUn/EmSVBjDX5Kkwhj+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/NEfbJyZZP76F1desY/34FiYmJgddkiTNiVf1k+ZofMNWXvDGy7l10zaWL1vERWcd45UCJQ0VW/7SHK1Zv4lbN20D4NZN21izftOAK5KkuTH8pTnaf+Uyli9bBMDyZYvYf+WyAVckSXNjt780Ryv2XMxFZx3DmvWb2H/lMlbsuXjQJUnSnBj+0hyNjS1g5V5LHOeXNLTs9pckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5L6yGtDqAlc5y9JfeS1IdQEtvwlqY+8NoSawPCXpFo/uuS9NoSawG5/Sap1o0t++8Qk4xu27nDth7GxBXfe77Uh1ASGvyTVpuuSn2v47+wEwmtDqAns9pekWje65B3T1zCw5S9JtW50yU+dQEy1/B3TVxMZ/pJU60aXvGP6GgaGvyR1kWP6GgaO+Uty17ku8r3UMLDlLw2RnS0jm69Sdp3r1fvXqpT3UsPN8JeGSK+CpRtL3IZBP4K5lPdSw81uf2mI9GoZWTd3nWvv9r7jjonGdIP3YxmeO/hpGNjyl4ZIr5aRdXOGemvreq89FnHBGY/ihedf0Yhu8H4sw3O2v4aB4S8NkV4FSzdnqLe2rgGuX7uxZ93gcx3D70cwO9tfw8Dwl4bIMARLa+sa4MBVe/SstT3XMfxheP+kfjD8JXVVe+t6KpR70dp2cp00P4a/pK6arnXdq9a2W+lK82P4SyOsH+vaB8nJddL8GP7SCBv1DWccw5fmx3X+0gjz8rKSpmP4SyPMDWckTcduf2mEOSYuaTqGvzTCHBOXNB27/SVJKozhL0lSYQx/SZIK0+jwj4gXR8Tgrv8pSdIIamz4R8QDgVcNug5JkkZNI8M/IhYBlwBfHXQtkiSNmkaGP/Bq4Drg3YMuRJKkUdO48I+IRwDPBk4cdC2SJI2ivm7yExEHA9fMcPdWYF/g/cBpmfnriOhXaZIkFaPfO/xdBxw6w30TwIXANzPzw/0rSZKksvQ1/DPzduAnM90fEScAWyJiY31ot/r4RuDkzPxgz4uUJGnENW1v//u0ff5E4A3A4cCa/pcjSdLoaVT4Z+bPWj+PiDXTHZckSfPXuNn+kiSptxrV8m+XmZdQbfYjSZK6xJa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwhj+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwhj+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpMIY/pIkFcbwlySpMIa/JEmFMfwlSSqM4S9JUmEMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhTH8JUkqjOEvSVJhDH9Jkgpj+EuSVBjDX5Kkwhj+kiQVZrdBF9AuIhYBrwWOB5YCXwT+JjOvGWhhkiSNiCa2/N8OPBU4DjgK2B24NCIWDLQqSZJGRKPCPyL+D3AC8KzMvCwzfwicAiwH7j3I2iRJGhVN6/Y/FrgpMy+bOpCZCdxzcCVJkjRamhb+hwA/j4jjgJcC+wJfBs7IzF8NtDJJkkZEX8M/Ig4GZpq4txW4BLgvcBZwRn3stcB/RcRhmbmlg5dZCHDDDTfscr2SJA2Dlsxb2MnX97vlfx1w6Az3TVAF/l7Ak6dm90fEk4HrgccBn+jgNQ4EOP7443e5WEmShsyBwNU7+6K+hn9m3g78ZKb7I+I6YFPrsr7MvDEi1gH36vBlvgE8guqEYfsulCtJ0rBYSBX83+jkixdMTk72tpw5iIijgSuA+2Xm6vrYAcCvgSdl5r8PsDxJkkZC08J/AXAlsCfwfGATcAHV2cxhmbltgOVJkjQSGrXOPzMngScA3wL+k2qm/zjwRwa/JEnd0aiWvyRJ6r1GtfwlSVLvGf6SJBXG8JckqTCGvyRJhTH8JUkqTNMu7NNVEbGI6toAxwNLgS8Cf9O6g6B2FBEvBl6fmQsGXUvTRMTvAq8HHgJsBj4FvCQz1w+0sAaIiIXAq6kuyb0n8Bng1MxcM8i6miYi9qf6HToW2B34GnBWfflyTSMijgC+BDwmM68YcDmNFBEnAi8B7gH8GHhx69VxpzPqLf+3A08FjgOOovrPdmm9mZDaRMQDgVcNuo4mioi7A1+gujDVkcBTgIcBHx1kXQ1yLvAs4JnAI4GDgI8PsqCmiYgx4N+orl76RKq/SeNUFy7bZ5C1NVVELAM+QIcXqylRRDwLuBh4HfBAqo3yLq0vpDejkQ3/iPg/VK2QZ2XmZfWZ9SnAcuDeg6ytiepekkuArw66loZ6GrAFeF5mrs7MLwOnAn8YEb892NIGq/7deSHwssz8fGZ+G/gr4OERcdRgq2uUw6hOHJ+TmV/PzB8DzwD2AP50oJU11/mAl3OfQd2QPQ/4h8x8T2b+DHgR8DOqk8sZjXK3/7HATa1dH5mZwD0HV1KjvZrqqosfAh412FIa6VLgm5nZerGoifp2b+AX/S+pMQ6n6uq/YupAZl4bEddSXWTrKwOpqnl+ATweyJZjE8ACqt8htYiIx1GdFD0W+P6Ay2mqoMq0j0wdyMwJqv+Tsxrl8D8E+HlEHAe8FNiXarvgMzLTM8kWEfEI4NlULZNHD7icRsrMq7nrZTJfSnXCVPp47UH17XVtx39NNQYpIDPXUW1b3uo0YAnwuf5X1FwRsQp4F/Ac4OYBl9Nkh9S3KyLiMuABVFfOPTszZz3pHtrwr8czZpq4t5WqC/u+wFnAGfWx11KNrx2WmVv6UeegdfA+7Qu8HzgtM38dEf0qrVF29j5l5pK2r38dVSvuz9t6A0q0FJioL9ndaitVsGkaEfEEqr9J509dxVR3ejvwycz8TEQctNOvLtfy+vZ9wDlUwX8icFlEPHi236uhDX+qVsahM9w3QRX4ewFPnprdHxFPBq4HHgd8oh9FNsDO3qcLqbqzP9y/khppZ+8TcOes9rcAJwOnZOalfait6W4DxiJit8y8o+X4Yqorc6pNRJwAvBP4F6pZ2qrVE9geDDxo0LUMgakT7r/PzA8BRMSpVMNtp1D1LE1raMO/bmX8ZKb7I+I6YFPrsr7MvDEi1gH36kOJjdDB+3QCsCUiNtaHdquPbwROzswP9rzIBtjZ+wQQEUuoZvf/CfD0qf9s4pf17YEtHwPcnbsOBRQvIl5ONcfmLVQ9bl5dbUcnUA0l3VD3RE6tzvp0RLwvM583qMIaaOr/1w+mDmTmZESsZic5N7Kz/anW9C+LiDtbcxFxALCKu47dluw+VMtDDq//vbw+fjjVJDdx5zKtfwX+EPgzg38H3wM2AEdPHaiHUQ4G/nswJTVTRLyEKvjPycwXGPzTejpwP37zN+mP6+MnUnVt6ze+TdW79tCpA/UKgPuxk5wb2Uv61m/AlVSzkJ9P9QZdQNU6OSwztw2wvMaKiKcDH3CTnx3VXWlvofoD1D5pa900491FqedAnFD/uxF4K7AlMx81uKqaJSIeRPXH+n385iR7yobMdIhkGvWY/y+BY9zk564i4lVUy45PpOoBeD7wPODweoXbtEa25V+fUT8B+BbVH+svU22o8UcGv+bh+Pr2XVTzRlr//f6gimqQVwAfpJpoeznwv8CTB1pR8/wV1WY1z+Guv0NnDLAuDbdzgH+katz+gGoviWNnC34Y4Za/JEma3si2/CVJ0vQMf0mSCmP4S5JUGMNfkqTCGP6SJBXG8JckqTCGvyRJhRnavf2lkkXEM4G/Ae5PdeGh7wMXZuZHWr5mEnhGZl7SoxreCxyUmY/p8OvvB9wrM9t3SJzLa74L+J3pdg6MiHtTbTV8Sfv+7xHxeOA/gD/OzC/M9/WlUWHLXxoyEXES1VbDbwUOo9ph8D+BD9dXRJtyIPCx/lc4o/+gZQ/ybsvMq4EzgZMj4nFTxyPiHlRb6r7e4Jcqtvyl4XMy8M7MfG/LsR9HdQm0F1IFHZl5wwBqm03PrxeRme+IiD8D3h0RDwRuobps7v8Af9fr15eGheEvDZ/twMMjYq/MHG85/iJg2dQnrd3+dRf97VTXtzixfo4LgE8A7wB+F0jgxMz8Zvvjp3vO9qIi4i+Bs4EHAJPAd4DTM/MbEXEFcG/glRFxQmYeHBGLgdcAx9V1fwd4aWZe1fKcp9bf1/51rZ30Vk5d4OQiYHVdz+GZeUcHj5WKYLe/NHz+EXgY8OuIuDQiXhQRh2fmTZl57SyPe2Z9+3vAm4D/C/w78Nr6+bYBF8+noIh4KPBR4L3AoVSX910AvLP+kicB1wJv5Ddd/+8HHgk8FXgIcBlweUQcUj/nM4DzqU4QHgz8CvjrndWSmWuAk6gupPN3wMmZec18vi9pVBn+0pDJzH8F/oBqnP+RVCcD34mIb0fE/Wd56E3Ai+ux8TfVxz6Umf8vM38A/DNVK3k+bgeen5kXZ+a1mfkNquB/YF3zeqreho2ZeVNE/A5V6J+QmV/MzJ9m5nnAl4Cz6uf8G6rJe+/MytnANzqs50rgZuAO4Cvz/J6kkWW3vzSEMvMrwFciYiFVS/7PgBcAn46I35nhstVX15e6JjM3VVMEuLrl/tuAxfOs57sRcUtE/C1wP+A+wOHM3MB4cH37tbqOKYtbangA9fyFFlcBD+qgpHcC66j+xr0/Ih6dmRMdPE4qguEvDZF65vrfAv83M2/IzO3A14GvR8QXgc9SheM3p3n47dMc6zgQI2LGvxcRcQzwaaoZ/V8G3gMcArxthodMnZwcSXXS0WprfTvJXScJTndS017LycBfAI+m+hv3earehH/c2WOlUtjtLw2X26gmtB03zX23UAXmjV16rduB5S2f32eWr30+8LnMfFpmXpiZlwMHA0TEVIBPtnz9j+rb/TPzZ1P/gDOAJ9b3fRc4qu11HjJbwRFxKNWQxhsy88rM/C+qiX+vjojDZnusVBJb/tIQycy1EfF64LURsRz4ONUJwQOBVwPvy8xfdOnlvgqcFBFfBhZSherWGb72JuBPI+IIYA3weOD0+r7FwBZgA3BIRNw9M38WER8B3lHP6P8p8BzgecCx9ePeAHwsIr5O1avwVKq5Dl+aroCIWEK1rC/ZcVnf2fVzXhIRD8nMmb4HqRi2/KUhk5mvoFrr/0dUQfgj4O+pxsdP6uJLnUK1NPBrVJsFvYNqxv10zqFaqvdZ4FtUs/unNhyamt1/PvBY4PsRMUbVg/EpqomGP6zve1LdWicz/x04ATiVagfDI+oaZvKPVEMNx7fOecjM24BnAPelWtkgFW/B5OTkzr9KkiSNDFv+kiQVxvCXJKkwhr8kSYUx/CVJKozhL0lSYQx/SZIKY/hLklQYw1+SpML8/yWLaP1lyq72AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x, y = sims.spiral_sim(num_samp=100, num_dim=1, noise=0.1)\n", "\n", "fig = plt.figure(figsize=(8,8))\n", "fig.suptitle(\"Spiral Simulation\", fontsize=17)\n", "ax = sns.scatterplot(x=x[:,0], y=y[:,0])\n", "ax.set_xlabel('Simulated X', fontsize=15)\n", "ax.set_ylabel('Simulated Y', fontsize=15) \n", "plt.axis('equal')\n", "plt.xticks(fontsize=15)\n", "plt.yticks(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Applying Dcorr on this data," ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dcorr test statistic: 0.012101670566301154\n", "P Value: 0.19970250056635008\n" ] } ], "source": [ "dcorr = DCorr(which_test='unbiased')\n", "dcorr_statistic, independence_test_metadata = dcorr.test_statistic(x, y)\n", "p_value, _ = dcorr.p_value(x, y)\n", "\n", "print(\"Dcorr test statistic:\", dcorr_statistic)\n", "print(\"P Value:\", p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, Dcorr has a p-value greater than the alpha value of 0.05, which is not desirable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A faster version of Dcorr (Fast Dcorr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dcorr operates in $\\mathcal{O} (n^2)$ time. To make things faster, Fast Dcorr can be used. This version of Dcorr takes bootstrapped subsamples of the input data matrices and estimates the null distribution using the subsamples. Fast Dcorr, thus, operates in $\\mathcal{O} (n)$ time.\n", "\n", "Calling Fast Dcorr is similar to calling Dcorr except that the `is_fast` flag should be set to `True` in either the `test_statistic` or `p_value` methods. The number of subsamples used can either be defined or can be chosen optimally from the provided data. Using the some linear data as before," ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Biased P Value: 0.001\n", "Unbiased P Value: 0.001\n" ] } ], "source": [ "x, y = sims.linear_sim(100000, 1)\n", "\n", "b_dcorr = DCorr(which_test=\"biased\")\n", "u_dcorr = DCorr(which_test=\"unbiased\")\n", "\n", "b_p_value, _ = b_dcorr.p_value(x, y, is_fast=True)\n", "u_p_value, _ = u_dcorr.p_value(x, y, is_fast=True)\n", "\n", "print(\"Biased P Value:\", b_p_value)\n", "print(\"Unbiased P Value:\", u_p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is able to calculate many more samples in a time faster than Dcorr!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mantel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we explore\n", "\n", "- The theory behind the Mantel test statistic and p-value\n", "- The features of the implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mantel is one of the earliest distance-based statistic for independence [[8]](https://cancerres.aacrjournals.org/content/27/2_Part_1/209), which achieves some success in testing without consistency testing (Mantel is later shown to be not consistent, see [[4]](https://projecteuclid.org/euclid.aop/1378991840)). Given $\\mathbf{x}$ and $\\mathbf{y}$, suppose that an appropriate distance measure is used to calculate two $n \\times n$ distance matrices ${\\mathbf{D}}^{\\mathbf{x}}$ and ${\\mathbf{D}}^{\\mathbf{y}}$ respectively, e.g., Euclidean distance. Let the centered matrix be calculated as,\n", "\n", "$$ {\\mathbf{C}_{ij}}^{\\mathbf{x}} = {\\mathbf{D}_{ij}}^{\\mathbf{x}} - \\frac{1}{n (n - 1)} \\sum_{i, j = 1}^n {\\mathbf{D}_{ij}}^{\\mathbf{x}},$$\n", "\n", "similarly for ${\\mathbf{C}}^{\\mathbf{y}}$. Consider a function of $\\mathbf{x}$ and $\\mathbf{y}$, $M$, such that,\n", "\n", "$$M_n (\\mathbf{x}, \\mathbf{y}) = \\mathrm{trace}({{\\mathbf{C}}^{\\mathbf{x}} {\\mathbf{C}}^{\\mathbf{y}}}).$$\n", "\n", "Then, the Mantel coefficient is calculated as,\n", "\n", "$$ \\mathrm{Mantel}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{M_n (\\mathbf{x}, \\mathbf{y})}{\\sqrt{M_n (\\mathbf{x}, \\mathbf{x}) M_n (\\mathbf{y}, \\mathbf{y})}}.$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using Mantel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, let's simulate some linear data:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAIlCAYAAABrdaqpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhkdXn3/3c3MAOzMQwgYCTiej8IGlA0ioLgQ/iJZtMY14iaIMa4Amo0ooASNQrIJkGUoAhRQB9jNIlLlEERWRIVUfFWcVBEBphpZpiF6Rmm+/fHqYKaorq7+nTt9X5d11zVfeqc6m8VRden7+82Mjk5iSRJ0myNdrsBkiSpPxkiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkCYiIkW63Qeo323e7AdKgi4jlwI6Z+fRpzjkZOAnYKTM3dahppUXELsDfA38O7AOMAzcDlwDnZ+b9lfMOA64EjsrMr3agXZ8CnpuZe87imhHgPcAm4MOVYyfTR/89pG6xEiH1hk8Cz6D4MO5pEbET8B3gxcDZwPOAlwPXAB+lCBJV36d4Xt/rcDNnYz5wCrCg5ljf/PeQuslKhNQDMvO3wG+73Y4mvQjYD/iDzPxRzfH/iIg1wCkR8aHM/GFm3gtc25VWzkGf/feQusYQIfWA+vJ5pSy/D3ABcCLwGOBW4AOZ+ema6+ZV7n8l8HDg18B5mXlm3eO/Gng9xYf/DsCvKuedU7l/H2AFcDzw18DewKmZeVqD5la7ChpVMi8ENgPrKo97GDXdGZV2XAQ8naJqcSBwJ3Ay8FXgHOAo4F7gosx8d137Xp+Z50/1utU3JiJGgeOAo4HHASPAz4APZeZlNY8LcFJEnJSZI40eNyIOBN4PPI2iavE94D2ZeW1dG18C/BnwfGA74L+AN2fmygavl9TX7M6QetcBFGX2DwJ/DPwG+FREPLHmnCuAE4CPV875PHB6RHywekJEvA74F+AbwJ9QVBJuBc6ufMjX+lDlsY4GvjJFu74K3A98LSL+MSIOrXRxkJm3Z+aHMvOWGZ7bF4DPVNpzK0X3wXLg58ALgW8C/xARL5rhcWbyj5V/F1N0u/wVRci5NCIeA9wBPLty7oUUXRgPERHPBq4DlgB/S/H6LACuiohD6k7/OHAPxev8zspz/Ngcn4fUk6xESL1rZ+CQzLwJICKSotLwp8BNEfGcytd/nZkXVa75RkRsAt4bER+rlOUfC5yVmSdWHzgivgusBg6n+PCu+nJmnjtdozLzpoj4C+B84B8q/7ZExHXAZcAFmbl5hud2emb+c6UtWyptuKGm8nAl8JfAMymCUVl7Aydl5unVAxGxAvhf4NDMvCgirq/c9dtqVaGBfwJuA46oPreI+A+KqsZpwB/WnPutzHxj5ev/joinAH8VESOZ6WZFGiiGCKl3rasGiIpqH/3Cyu0fVW7/PSJq/1/+EkUF4/8Cn87MtwNExBLg8RSh4imVc+fV/cybaEJm/ntE/BdwaOXnHELxQfos4HURcVhmrp7mIb5b8/WdldsHPsAzc0tErAV2aaY907TzrwAiYlcefO6HV+6uf+4NRcRCii6MD9eGo8wcj4jLgbdHxKKaS75b9xC/pehC2oGiCiINDEOE1Ls21n6TmRMRAQ92Q+5WuV01xfW/BxARjwL+GTgS2ErRZVD9oKtfG+FOmpSZWyi6Hb5Z+TlLgHdRlPD/HnjHNJff2+DYhrrv5/xXe0Q8GTiXB2da3MyDQanZdSGWVs5tNKbhjsp9S2qObaw7Z6Jya/exBo5vaql/raH4y/ZpwFMb/PtUZQ2ErwCPpKgSLMzM/YC3lv2hEXFNRFxWfzwz783Md1GElP3KPv4UqoFiu7rji6e6ICIWU4zfmKAYX7IwMw+kGPcxG2sqP7/R2hMPr9w3XdVFGliGCKl/LacoyS/IzP+p/qP4y/lUYC9gd+AJwKcy85qacvzzKrdlfgfcAvxZRDyh/o6IWAo8DPjRQ66am2rlYu+64/WDGmvtS/H8z8nMGzNza+V4/XPf+pAra2TmBuB64C8jYofq8YiYT7FWxnWZ6XoSGkp2Z0idsWdENPrr/47MfMhf9U36L4rpk1dExD8CN1J8cJ5KUWb/caXffgXwt5Xbuyk+eP+e4i/ohQ0feXr/ABwGXBMRH6NYeGpj5We/FRgDTp/y6hIy856IuJriedwM3E4xFfVR01z2M2At8M6IuK/SxucB1UGPCyuPvSUi1gMHR8ShledT713A14FvRsRHK8eOp+gyevVcnpvUz6xESJ3xSIp1Eer/HVf2ATNzgmJa50WVx/k6xQf8FcDhNX8d/ynFuhCfpJha+XzgGIoQcmiJn3sbxfoOn6BYD+EKiumjbwP+E3haZk41TmMuXkUxluM84HMUYeXvp2nnvRTPfRz4bOXfkyles5/y4NROKNZ/eCrFa1Jf7SAzrwSeQzG19RLgUxRjOA7NzOVzelZSHxuZnHTGkSRJmj0rEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKmX7bjdgJhHxcWC7zDxmmnOuAF5Ud/ibmXlEWxsnSdIQ69kQEREjwCnAscCFM5y+P/BO4NM1x8bb1DRJkkSPhoiIeDRFcNgf+M0M584DHgtcn5krO9A8SZJE746JeAbwK+CJwIoZzt2XIgzd3O5GSZKkB/VkJSIzLwUuBYiImU7fH9gMnBIRRwH3AVcAp2bmpna2U5KkYdaTIWKW9gNGgATOpahenAHsDbyqmQeIiPnAU4E7gK3taaYkST1jO2Av4IbMLD2GcBBCxInAaZk5Vvn+pojYCnwuIo7PzNVNPMZTge+0rYWSJPWmQ4Cry17c9yEiMyeAsbrDN1Vu9waaCRF3AFx66aXsueeeLWydJEm9Z+XKlbziFa+AyudfWX0fIiLicmCHzHxBzeGDKKZ4/rLJh9kKsOeee/KIRzyixS2UJKlnzakLv+9CRGVK5zJgLDM3A5+n0nUBfAk4EDiNootjffdaKknSYOvVKZ7TOZii/HIwQGZeDrwaeA3wY+B04CzgvV1qnyRJQ6HnKxGZeVjd98spZmPUHrsYuLhzrZIkSf1YiZAkST3AECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkqS+tHVikrG1m7h5xWrG1m5iYmKy200aOtt3uwGSJJWxdt04bzr9Su7dsJklC+dxzgmHs2znHbvdrKFiJUKS1JfuHNvAvRs2A3Dvhs3cObahyy0aPoYISVJf2mPZQpYsnAfAkoXz2GPZwi63aPjYnSFJ6itbJyZZu26cNes3cdbxh3H3PRvZY9lCli6e3+2mDR1DhCSpr9SOhdh5UTEWYpcljoXoBrszJEl9pXYsxNr1m1m52rEQ3WKIkCT1FcdC9I6e786IiI8D22XmMdOccxBwFnAgcDvw/sy8uENNlCR10NLF8znnhMO5c2yDYyG6rGcrERExEhHvA46d4bzdga8B3weeDJwNXBgRR7a/lZKkThsdHWHZzjuy76N2ZdnOOzI6OtLtJg2tnqxERMSjgQuB/YHfzHD6McBa4C2ZOQH8LCKeDLwN+HpbGypJ0hDr1UrEM4BfAU8EVsxw7iHAtysBomo58MyI6NXnJ0lS3+vJSkRmXgpcChARM53+COAHdcd+BywAlgGrWt0+SZLUu5WI2VgAbKo7Nl65deKwJEltMggh4j6gfmhu9XsnD0uS1CaDECJuA/aqO/ZwYD3FgEtJktQGgxAirgYOjYjaOT6HA9+tG2wpSZJaqCcHVk4nIuZRDJgcy8zNFFNB3wGcHxFnAkcALwee271WSpI0+PqxEnEwcEfllsy8kyIwHEgxS+ONwNGZ+a2utVCSpCHQ85WIzDys7vvlwEjdsWuBp3WuVZIkqR8rEZIkqQcYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSS2xdWKSsbWbuHnFasbWbmJiYrLbTVKb9fyy15Kk/rBm3SbefPpy7t2wmSUL53HOCYexbOedut0stZGVCElSS6xcvZF7N2wG4N4Nm7lj9cZZXW8lo/9YiZAktcRuS3diycJ5D1Qidls6uyrE2nXjvOn0K2sqGYezbOcd29RatYIhQpLUEvN3GOWU1z6Du+7ZyMN2WcD8HWZX7L5zbENdJWMDSxfPZ3T0wY2bt05MsnbdOHeObWCPZQsfcr86yxAhSWqJJQvnMzEBExMTLFuyI0sWzp/V9XssW7hNJWP+DtuxZt34NtUIqxW9xRAhSWqJ0dERlu28Y+kP9aWL53PGW5/NLb9dw+5Ld+Kf/9+POOZP99vm8eqrFXeObTBEdJEhQpLUE0ZHR9hhu1Eu/s+buWPVehYtmMceyxZuc059taL+fnWWIUKS1DOWLp7PB17/zG3GPNTff84Jh095vzrLECFJKq3VAx1n6hKZa5eJWssQIUkqzYGOw83FpiRJpTUa6KjhYYiQJJVWHegIONBxCNmdIUkqzYGOw80QIUkqzYGOw83uDEmSVIohQpIklWKIkKQB5Lba6gTHREjSAHL9BnWClQhJGkD9sH6D1ZL+ZyVCkgZQP2xUZbWk/xkiJGlA1O5jsdduCznnhMO4c2xj6fUbWr0vRj239e5/hghJGhCN/rLf91G7tvTxWvkh3w/VEk3PECFJA6LVf9m3u1Lgapf9zxAhSQOi1X/Zt7tS4GqX/c8QIUkDotV/2Vsp0EwMEZI0IFr9l72VAs3EdSIkSdNyPQdNxUqEJA2gVk7PdD0HTcUQIUkDaKoP/jLhwvUcNBVDhCQNoKk++MtUFVzPQVMxREjSAJrqg79MVcFZGpqKIUKSBtBUH/xlqgrO0tBUDBGS1CHt3oui1lQf/FYV1EqGCEnqkF6Y5WBVQa3kOhGS1CGNxiNI/cwQIUkdUh2PADjLQQPB7gxJ6hDHI2jQGCIkqUMcj6BBY3eGJA0Y97pQp1iJkKQB0wuzQDQcrERI0oCZyywQqxiaDSsRkjRg5rLXhVUMzYYhQpIGzFxmgbhjp2bDECFJA2Yus0AetmxBXRVjQRtaqEFhiJAkPWC70RFOee0zuOuejTxslwVt29tDg8EQIUl6wB2rNvDOj13NogXzWL9xMx96w7NYutjuDDXm7AxJ0gP2WLaQRQuK7oxFC1yaW9OzEiFJeoBLc2s2DBGSpAe4NLdmw+4MSeohLvakfmIlQpJ6iIs9qZ9YiZCkNptNdaHMktVWL9QtPVmJiIjtgFOBVwOLga8Cb8jMO6c4/wrgRXWHv5mZR7SznZLUjNlUF8osWW31Qt3SkyECOBl4FXA0sBo4D/gC8Kwpzt8feCfw6Zpj421snyQ1bTZLSZeZHeFS1eqWngsRETEPeAvw5sz8RuXYS4EVEXFwZl7T4PzHAtdn5sqON1iSZjCb6kKZ2RFz2XBLmoueCxHAARRdGMurBzLz1oi4FTgEuKbu/H0pnsfNnWmeJM1Ou9decG0HdUsvhohHVG5vrzv+O2DvBufvD2wGTomIo4D7gCuAUzNzU9taKUlNavfaC67toG7pxRCxAJjIzC11x8eBRv+H7AeMAAmcCzwROIMicLyqje2UJGmo9eIUz/uA0YioDzjzgUZznU4E9szMMzLzpsz8V4oxFUdHxK5tbqskSUOrFysRt1Vu96r5GuDhPLSLg8ycAMbqDt9Uud2bYnaHJElqsV6sRNwIrAOeXT0QEfsA+wDfrj85Ii6PiC/WHT6Iovvjl21rpSS1gAtFqZ/1XCUiM8cj4jzgtIhYBdxFsU7EVZl5bWVK5zJgLDM3A58HPhcRxwNfAg4ETgNOy8z13XkWktQcF4pSP+vFSgQU4xwuBS4BrgR+zYMrUh4M3FG5JTMvp1jZ8jXAj4HTgbOA93a0xZI0jakqDmWWuZZ6Rc9VIgAy837ghMq/+vuWU8zGqD12MXBxRxonSSVMVXFwoSj1s54MEZI0aKZamnq2C0VtnZhk7brxbc4fHR2Z9hqpXQwRktQBU1UcZrtQlGMo1EsMEZLUAa1amtrNttRLDBGS1AGtWpraMRTqJYYISeojbralXmKIkKQ+4mZb6iW9uk6EJEnqcYYISepzLp2tbrE7Q5L6nNM+1S1WIiSpz7l0trrFECFpILWqxN8PXQXVaZ+A0z7VUXZnSBpIrSrx90NXgdM+1S2GCEkDqVUrO/bDCpFO+1S32J0haSC1qsRvV4E0NSsRkgZSq0r8SxfP5+wTDmPl6o3stnQnRkdhYmLSnTMlDBGSBlSjEn+ZbbRHR0cYYYRzLv8hd6xaz6IFvTkuQuoGQ4SkoVF2kOSdYxu4/e71QO+Oi5C6wTERkoZG2fUUWjkuoh+mjErNshIhaWiU3Ua70fiKMl0jWycmuefeTbzljOU9PWVUapYhQtLQKDvYstH4irG1m2bdNXLv+nF+cds9PT9lVGqWIULS0Gjlegpl1o9YuXoDu+28U6lqiNSLDBGSVEKZrpE9li3ko5/7X04+5umsWnsfj9t7F1eXVF8zREhSCWW6RpYuns9xL30Ka9Zt4vG/vwt3jW1kdGSkqfEUUi8yREjSFKYbPFmma6R6DdDz+3FIzTBESNIU2rX5Vj/sxyE1w3UiJGkKZdeVmIn7cWhQWImQpCmUXVdiJm7drUFhiJCkKbTrw96tuzUoDBGSNAU/7KXpOSZCkiSVMmWIiIg3dLIhkiSpv0xXiTgrIv47IvbuWGskqcXcNVNqn+lCxLOAPYAfR8TfdKg9ktRS1bUe3nHu1bzp9CtZs268202SBsaUISIzrwUOBD4CnBsR/xkRv9exlklSC7RrrQdJM8zOyMz7gVMj4rPAWcBNEXEesLHuvA+0r4mSVF671nqQ1PwUz98CNwBHAn8N1NYDJwFDhKSe5MJOUvvMGCIi4nnAx4BdgeMy82Ntb5Uktch0az1Mt8GWpJlNGSIi4mHA2cBfAv8NHJuZv+5UwySp3dq1wZY0LKarRPyscntMZl7UicZIUie5m6Y0N9NN8bwKeIIBQtKgcjdNaW6mrERk5gs62RBJ6jQHXUpz4wZckoaWG2xJc+MGXJI6zqWopcFgJUJSxzkrQhoMViIkdVy7lqK2wiF11nTrRHyr2QfJzOe0pjmShkG7lqK2wiF11nTdGbfUfL0AeBnF0tfXApuBgyh2+vx421onaSC1a1aE6z5InTXdFM/XVr+OiM8A/5SZ76o9JyJOBJ7WvuZJGkTtmhXhZltSZzU7sPIFFNuC17sM+IfWNUeSynPdB6mzmg0RdwGHAL+oO34UcFtLWyRJJbnug9RZzYaI04HzIuLpwPeBEeBg4CXAq9vTNEmS1MuaChGZ+bGIuBd4A/ByYBL4IfCXmfmlNrZPkiT1qKYXm8rMzwCfaWNbJElSH2k6RETEY4ETgP8D/BXwZ8BPM3N5e5omSQ/aOjHJ2nXj2wyaHB0d6XazpKHW1IqVEfGHwI3AYyjGQswH9gO+ERF/3L7mSVKhupDUO869mjedfiVr1o13u0nS0Gt22esPAR/KzCMpFpoiM98AfBg4uT1Nk6QHtWupbEnlNRsingx8tsHxC4F9W9ccSWqsupAU0NRCUu6jIbVfs2MiNgAPA35Zd/zxwNqWtkiSGpjtQlLuoyG1X7OViM8CH42IfSmmd+4YEc8BzgU+367GSepvZasBja6rLiS176N2ZefF81mzbnzax7X7Q2q/ZisR76KY3vmTyvc3VW4vr9wnSQ9Rthpw7/px3nzGlaxd3/i6Zh7XfTSk9mt2sanNwEsi4t3AARSDK3+SmbdMf6WkYTbbXTWr0zhvv3sdJ/3N0zn/izfx89/c85Drmnlc99GQ2q+pEBERvwIOysxfUjMuIiL2Am7MzIe1qX2S+thsqwH1FYaTj3k6p1x47UOua+Zx3UdDar8pQ0REPA84qPLtPsA7I2J93WmPn+4xyoqI7YBTKfblWAx8FXhDZt45xfkHAWdR7DR6O/D+zLy41e2SNDuzrQbUVxg2bb6fc044nJ0XbXudVQapN0wXAFYAZ1JstgXwImBrzf2TwDrgTW1o18nAq4CjgdXAecAXgGfVnxgRuwNfA/4V+Bvgj4ALI2JlZn69DW2T1KTZVgPqKwy/t/tidlny0GutMki9YcoQkZk3U1QaiIgrgRdm5j3tblBEzAPeArw5M79ROfZSYEVEHJyZ19RdcgzFNNO3ZOYE8LOIeDLwNsAQIfURKwxSf2lqimdmHt4oQETEvIh4ZovbdABFF8bymp9/K3ArcEiD8w8Bvl0JEFXLgWdGRLNTWCX1gNppnMt23tG9MaQe1+zAyqcAnwCeSOPgsV0L2/SIyu3tdcd/B+w9xfk/aHDuAmAZsKqFbZMkSRXN/qV+FnAfcCzF9M7XU+ybMQ68pMVtWgBMZOaWuuPjQKMO0AXApgbnMsX5knpMKxelktQ5zc6sOBA4NDP/NyKOBTIzL4iI31EEilauWnkfMBoR22fm/TXH51Msv93o/PqO0+r3LlEn9ZhGW3o3uyhV/bWjo7i0tdRFzYaIEeDuyte/oOjWuAr4MnBSi9t0W+V2r5qvAR7OQ7s4qufvVXfs4cB63NdD6jn1geHctx3e9KJU9dee8dZns35j84tZSWqtZrszfgw8r/L1T4HqYMo9aO14CIAbKaaOPrt6ICL2oVir4tsNzr8aODQiakdgHQ58t26wpaQeUB8YVq7esM0OnQc8fjd232VBwy6K+mtXrbmPvXZbBDS3s6ek1mq2EvFPwGURsZViM66TIuLfgD8ArmxlgzJzPCLOA06LiFXAXRTrRFyVmddWpoAuA8Yqy3FfCLwDOD8izgSOAF4OPLeV7ZLUGo1Wm6yd2rn7Lgt4yxnLG3ZR1F+7164L+MDrn+mUUKlLmt074wsR8XRgS2b+urKa5euB/wTe24Z2nQjsAFxSuf0q8IbKfQdTBJfDgeWZeWdEPBc4m2KWxq+BozPzW21ol6QmNBr3UJ2u2WgtiNrFo25esXrKro3prpXUeU0vWZ2Z/1Pz9ZW0uAJR97PuB06o/Ku/bzkPrqJZPXYt8LR2tUfS7Ew3UHK61Sa3Tkyy+y4LptwXw5Uqpd4y3d4ZTa/2mJlHtqY5kgbBbHfvrFq7bpyzLvs+Jx/zdFatvY/H7b2LXRRSD5uuEtFoJoQkzWi2u3dW3Tm2gR/+fBU/+sW3WbRgHie+5mnstnSnNrdWUlnT7Z3xmk42RNLgKLsHRm34qH5fa7qxFpI6r9llr18+3f2Z+a+taY6kQVB27MJM4aPZRakkdUazAysvmeL4JuC3FNtwSxoA3fxrf6bwUXashaT2aHaK5zaLUkXEdhTbhP8z8PE2tEtSl/TyX/tlx1pIao+mp3jWysytwM0RcTxwOcUCVJIGQC//tV92rIWk9igVImrcT7FPhaQB0Yq/9tvVJeI6EVJvmcvAyiUUW4Nf19IWSeqqufy1Xw0PW7ZOcPyZV/Vkl4ik1pnLwMotwPeAv2tdcyR121z+2l+7bpx3n/9dXnnUvj3bJSKpdUoNrJSkRu4c28Dv7l7P7kt3cgCkNATmOiZC0gAqO6Zhj2ULWbRgHud/8SZOee0zGN+ylb12dQCkNKiaHRPxFOBcYH/gIb8NMnNei9slqYvKTvOsHU+xbMmOrigpDbhmKxGfBDYDbwfua19zJPWCstM8nT0hDZdmQ0QAT83Mn7SzMZJ6Q6Npnu5bIalesyHi+8DvA4YIaQg0mua5podXspTUHc2GiGOBL0bEU4FfARO1d7oBlzRYGnVL9PJKlpK6o9kQ8ULgccDJDe6bxA24pIHXrX0r7EaRelezIeKtwInAmZm5sY3tkdSj6rs4liycx9jaTW3/cO/lDcGkYddsiNgO+KwBQhpe9V0cY2s3deTD3W4UqXc1uxLlJ4G/bWdDJPWGrROTjK3dxM0rVjO2dhMTE5MNz2v04d4O1aoH4OqXUo9pthKxM3B0RLwMuIVi34wHZOaRrW6YpIfqxPiAZrsPOjVGwu2/pd7VbIjYAfhsOxsiaWadGB/QbPdBpz7cXcBK6l3NbsD1mnY3RNLMOjE+oNkKgx/ukqYMERHxcuDzmbm58vVUJjPTKoXUAZ3oQrD7QFKzpqtEXAL8N3BX5eupTGJXh9QRnfiAt8IgqVlThojMHG30taTu8QNeUi9pdmDlAyJie+BJwJ2ZeXvrmyRJkvrBtBWGiHhlRPxPRPx+5fsnAL8EbgB+HRGfjIjtOtBOSZLUY6YMERHxYuBTwI+B6ioynwGWAM8FDgaeTrEktiRJGjLTdWe8GTgxMz8IEBF/ABwIvD8zv1E59h7g/cDp7W6opP7nZlrSYJkuRDwJOKbm+yMoZmJ8uebYj4DHtKFdkgaQm2lJg2W6MRGjwOaa7w8F1gL/W3NsJ2BTG9olaQB1ar8NSZ0xXYj4CfBMgIhYAvxf4OuZWbsbz19QjJmQpBm5mZY0WKbrzvgYcE5EPIkiTOwEnAkQEQ8DXg68E3hduxspaTC4GqY0WKZbbOriiNgROBbYCrwkM6+t3H0SxXiJD2fmxe1vpqRuatWASBfLkgbLtItNZeYFwAUN7vog8N7MXN2WVknqKQ6IlNTIrFesBMjM37a6IZ1nTiMAABf+SURBVJKm183pkZ3YPVRS/ykVIiR1Xn014OwTDmOEkY6Eik7sHiqp/xgipDZpdeWgvhqwcvVGzrn8h9x+9/q2dzE4IFJSI4YIqU1aPY6gvhqw29KduGPVeqD9XQwOiJTUiCFCapPZjCNopmpRXw0YHYVFC+xikNQ9hgipTWYzjqCZqkV9NWBiYtIuBkldZYiQ2mQ24wjKzH6wi0FStxkipDaZzYe8sx8k9SNDhNQDaqsWe+26kInJSW5esdrtsiX1NEOE1IR2L/RUW7UYW7uJN52+3NUhJfU8Q4TUhE4u++zqkJL6xXRbgUuqaPTB3i5uly2pX1iJkJrQiYGP1S6TNes3cdbxh3H3PRuduimppxkipCZ0Ytnn2i6TnRcVXSa7LLEbQ1LvMkRITZhpumYrBl7WdpmsXb+Zlas3GCIk9TRDhNQCrRh46VoRkvqNIUJqgVbMqHCnTEn9xhAhtUArqgguYy2p3xgipBawiiBpGBkipBawiiBpGLnYlCRJKsUQIUmSSjFESJKkUgwRUo/ZOjHJ2NpN3LxiNWNrNzExMdntJklSQw6slHpMJ3cMlaS5sBIh9ZhO7hgqSXPRc5WIiHgYcC5wJLAZuAh4d2beP801dwG71x1+T2ae2raGSm3i8teS+kXPhQjgC8Ak8Gzg94BPAfcD7250ckTsQREgDgV+UXPXura2UmoTF66S1C96KkRExDOAZwGPzswVwI0R8XbgnIh4X2aON7hsf4qQcV1mbu5gc6W2cOEqSf2i18ZEHAL8uhIgqpYDi4EDprhmf+AWA4QkSZ3VU5UI4BHA7XXHfle53Ru4rsE1+wP3R8RXgIMq15+ZmZ9pWyslSVJnQ0RE7AOsmOLuceASYFPtwczcEhGTwFS13f2AXYH3UIybOAq4KCK2z8yLWtFuqV22Tkyydt34NuMfRkdHut0sSWpKpysRtwP7TnHfBPAmYJtRZBGxAzACTDXP7XBgXmZWB1LeGBGPBI6nmNkh9SzXhJDUzzoaIjJzC/Czqe6PiNuA59Udfnjltr6bo/qY4xRVjFo3AS8r2UypYxqtCWGIkNQvem1g5dXAoyNi75pjh1NM1/xh/ckRsX1E3BYRx9XddRDwk/Y1U2qN6poQgGtCSOo7vTaw8nvAtcBlEfFGYA/gn4AzqrMvImIRsCgzV2bm/RHxZeDEiLgF+Cnw58Arged35RlooLV6DINrQkjqZz0VIjJzMiJeAPwz8B2KCsSFwPtqTnsbcBLFOAmA44B7gLOBvSi6S16cmV/vVLs1PNauG+fd53+X3929nkUL5j6GwTUhJPWzngoRAJm5EnjBNPefDJxc8/04xayMhitaSq2ydWKSLVsneOVR+7L70p04/4s3OYZB0lDruRAh9aq168Y5/syrHphJccprn8GyJQYIScOr1wZWSj2rfibF+JatjmGQNNQMEVKT6mdS7LXrQheGkjTU7M6QmuRMCknaliFCapIzKSRpW3ZnSG20dWKSsbWbuHnFasbWbmJiYrLbTZKklrESIbWRe2NIGmRWIqQ2arQ3hiQNCkOE1EbujSFpkNmdoaHX6v0wajmjQ9IgM0Ro6NWPWzj7hMMYYaQlocIZHZIGmSFCQ69+3MLK1Rs55/Ifcvvd6x0MKUnTcEyEhl79uIXdlu7EHavWAw6GlKTpWInQ0KsftzA6CosWzHuge8PBkJLUmCFCQ69+3MLExKSDISWpCYYIqY6DISWpOY6JUE8ouzy0y0pLUvdYiVBPKLs8tMtKS1L3WIlQTyi7PHTtdXvuupAtWyesSkhSh1iJUE+oTrOc7YyI6nXrN27m9S98EsefeVXpqkQ7V66UpEFkiFBPKLs8dPW6VWs2Mr5l60OqGbMJEXaNSNLsGCLUE8rOiKi9bmztplLVjKpGXSqGCEmamiFCfalR18NcN7sq26UiScPKEKG+NFXXw1zWd1i6eD5nn3AYK1dvZLelOzE6Wiw85bgISWrMEKG+1I6uh9HREUYY4ZzLf8gdq9azaIHjIiRpOk7xVF+q3zSrVV0Pd45t4Pa71zMx6eZbkjQTKxHqS3Md/zAVx0VIUvMMEepL7drfol3hRJIGkSFCquHmW5LUPMdESJKkUgwRkiSpFEOEJEkqxRAhSZJKMUSo67ZOTDK2dpNbeEtSn3F2hrrO3TMlqT9ZiVDXNVrCWpLU+wwR6rp2LWEtSWovuzPUda4SKUn9yRChrpvLKpFbJyZZu258mwDi1t2S1BmGCPU1B2VKUvc4JkJ9zUGZktQ9hgj1NQdlSlL32J2hvuagTEnqHkOE+ppbd0tS99idIUmSSjFESJKkUgwRkiSpFEOEJEkqxRChOXEbb0kaXs7O0Jy4YqQkDS8rEZqT2a4YaeVCkgaHlQjNSXXFyGolYqYVI61cSNLgMERoTma7YmSjyoUhQpL6kyFCs7J1YpJ168cZv3+CVWvuY69dF7B0cfMrRs62ciFJ6l2GCM3K2nXjjN27iZM+8b1SXRLudSFJg8MQoVlZtWYjq9ZuKt0l4V4XkjQ4DBGald2WLmB0dHSbLondd1nAxMQko6Mj3W6eJKmDDBGalaWL5zM6Cme89dncdc9GFu64A2dd9n2Oe+lT2lJd2Doxydp149t0fxhWJKk3GCI0K6OjIyxdvCM///UYH774Bu7dsJmJSdo2y8IpoZLUu1xsSqXstnQBE5MwMUlbZ1nMdjErSVLnWIlQKa2aZTFTd8Ueyxay86J523wvSeoNhgiV0qpZFjN1VyxZOI8zjzuMO1atZ6/dFrFk4bxpHk2S1Ek9GyIiYj5wPfCRzLxkhnNfAbwX+H3gRuBNmXlD+1upuZppBct7N2zmLWcsd0yEJPWgnhwTERGLgS8CT2ri3COAfwFOB54M3AR8PSJ2b2sj1RLVFSyh8dgKx0RIUu/quUpEJRR8HFjT5CVvBz6bmRdUrn8d8BzgtcAH2tJItcxMYytcJluSelfPhQjgeRSVhdOATdOdGBGjwDOBN1aPZeZERHwbOKSdjVRr1I+tqG4VXg0V1S4Ml8mWpN7TcyEiM4+vfh0RM52+FFgI3F53/HfAU1vbMnXCVAMtHQchSb2noyEiIvYBVkxx93hmzvaTYkHltr5iMQ74qdNFZVeadKtwSeofna5E3A7sO8V9EyUe777KbX2Nez7gCLwuKrvSpGMgJKl/dDREZOYW4GctfMgxirCwV93xh/PQLg51UNmKgluFS1L/6Mkpns3KzEngGuDZ1WOVwZaHAt/uVrs089TNqVQHWu77qF1ZtvOObrYlST2s5wZWziQiFgGLMnNl5dAZwJcj4gfAt4DjgZ2BT3apicKKgiQNg36sRLwNuKP6TWZ+FTgWOAH4PvAE4MjMXNWd5g226hTMm1esZmztJiYmJhueZ0VBkgZfT1ciMvMhnzyZeTJwct2xi4CLOtOq4ebW3JKkqn6sRKiLXIZaklRliNCslB0wKUkaPD3dnaHe44BJSVKVIUKzUr/XhSRpeNmdIUmSSjFESJKkUgwRkiSpFEOEJEkqxRChpjS7UqUkaXg4O0NNcaVKSVI9KxFqiitVSpLqGSLUFFeqlCTVsztDTXGlSklSPUOEmlK/UmV1oGVtqHC7b0kaLoYIleJAS0mSYyJUigMtJUmGCJXiQEtJkt0ZKsWBlpIkQ4RKcUtwSZLdGZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRStu92A6YSEfOB64GPZOYlM5x7F7B73eH3ZOap7WqfJEnDridDREQsBi4DntTEuXtQBIhDgV/U3LWuPa2TJEnQgyEiIo4APg6safKS/YH7gesyc3PbGiZJkrbRi2Mingf8C3Bwk+fvD9xigJAkqbN6rhKRmcdXv46IZi7ZH7g/Ir4CHATcDpyZmZ9pTwslSRJ0OERExD7AiinuHs/MHUs87H7ArsB7gHcDRwEXRcT2mXlRk4+xHcDKlStL/HhJkvpLzefddnN5nE5XIm4H9p3ivomSj3k4MC8zqwMpb4yIRwLHA82GiL0AXvGKV5RsgiRJfWkv4JayF3c0RGTmFuBnLX7McWC87vBNwMtm8TA3AIcAdwBbW9Q0SZJ61XYUAeKGuTxIz42JmI2I2J6ie+SMzPxozV0HAT9p9nEqQeTqFjdPkqReVroCUdV3ISIiFgGLMnNlZt4fEV8GToyIW4CfAn8OvBJ4fjfbKUnSoOu7EAG8DTgJGKl8fxxwD3A2RWnmZ8CLM/Pr3WmeJEnDYWRycrLbbZAkSX2oFxebkiRJfcAQIUmSSjFESJKkUgwRkiSpFEOEJEkqpR+neLZMRMwHrgc+kpmXzHDuXcDudYffk5mntqt9vWKWr9MrgPcCvw/cCLwpM+e0Ilo/iIiHAecCRwKbKZZcf3dm3j/NNQP/noqI7YBTgVcDi4GvAm/IzDunOP8g4CzgQIpl8t+fmRd3prXdVeK1ugJ4Ud3hb2bmEe1sZy+JiI8D22XmMdOcM7TvqaomX6dS76ehrURExGLgi8CTmjh3D4pf9odSrEVR/ffR6a4bBLN8nY6g2Mb9dODJFMuPfz0i6j8oB9EXgD2BZ1N8CLwGOGWqk4foPXUy8CrgaIrn+giK1+ohKu+TrwHfp3j/nA1cGBFHdqSl3XcyTb5WFfsD72Tb989ftreJvSEiRiLifcCxM5w31O+pZl+nilLvp6GsRFQ+7D4OrGnykv2B+4HrMnNz2xrWY0q8Tm8HPpuZF1Sufx3wHOC1wAfa0sgeEBHPAJ4FPDozV1BsAvd24JyIeF9lWfV6A/+eioh5wFuAN2fmNyrHXgqsiIiDM/OaukuOAdYCb8nMCeBnEfFkigXmBnrxuNm+VpXzHwtcn5lDtf1wRDwauJDi/6HfzHD6ML+nmn6d5vJ+GtZKxPMo/mI+uMnz9wduGdRf9tNo+nWKiFHgmcDy6rHK/7TfptjcbJAdAvy6EiCqllOUpA+Y4ppheE8dQPEaLK8eyMxbgVtp/J44BPh25X1TtRx4ZuX9Nchm+1rtS/FH4M3tb1rPeQbwK+CJFHsnTWeY31OzeZ1Kv5+GshKRmcdXv46IZi7ZH7g/Ir5CsbnX7cCZmfmZ9rSwN8zydVoKLKR4bWr9Dnhqa1vWcx5B4+cNsDdwXYNrhuE99YjKbaPXZu8pzv9Bg3MXAMuAVS1tXW+Z7Wu1P8XYm1Mi4ijgPuAK4NTM3NS2VvaAzLwUuBSa+r00tO+pWb5Opd9PAxciImIfpk5d45m5Y4mH3Q/YFXgP8G7gKOCiiNg+My8q1dAua8PrtKByW/+GGwfKvOY9Y6bXCriEuuedmVsiYpKpn/vAvacaWABMZOaWuuNTvScW0Pj9wxTnD5LZvlb7UewflBQDep8InEEROF7Vxnb2m2F+T81G6ffTwIUIiiS/7xT3TUxxfCaHA/Myc13l+xsj4pHA8RSj8PtRq1+n+yq38+uOzwc2lHi8XjLTa/Um6p53ROxA8T/lVM99EN9T9e4DRivBqHaWylTvifto/P5hivMHyWxfqxOB0zJzrPL9TRGxFfhcRByfmavb3N5+Mczvqdko/X4auBBRSfI/a/FjjvNgeq26CXhZK39OJ7XhdRqj+J9yr7rjD+ehJdq+MtNrFRG3UYwfqfXwym3D5z6I76kGbqvc7lXzNUz9nriNxu+f9RSD4wbZrF6rSh//WN3hmyq3ewOGiMIwv6eaNpf306APLJmziNg+Im6LiOPq7joI+Ek32tSLMnMSuIZiiiPwwGDLQykGVw6yq4FHR0Rt3/XhwDrgh/UnD9F76kaK16D2PbEPsA+N3xNXA4dGxEjNscOB79YNjBtEs3qtIuLyiPhi3eGDKILpL9vWyv4zzO+pps3l/TRwlYhWiIhFwKLMXJmZ90fEl4ETI+IW4KfAnwOvBJ7fzXZ2W+3rVDl0BvDliPgB8C2K0vzOwCe71MRO+R5wLXBZRLwR2AP4J+CM6uyLYXxPZeZ4RJwHnBYRq4C7gPOAqzLz2sq0smXAWOV1uhB4B3B+RJwJHAG8HHhud55B55R4rT5PpdQMfIliIaXTKErS67vzLLrP91RzWvl+shLR2NuAO2q+Pw44n2Khkp9Q/LJ/cWYO9DzjJmzzOmXmVykWNTmBYnGXJwBHZubAjoCGB6owLwDuBL5DMabhQuB9NacN63vqRIoR4pcAVwK/5sFV8Q6meE0OBqiszPhcil9gPwDeCBydmd/qcJu7ZTav1eU8uKjZjykWeDuLYrXYYeZ7qjktez+NTE5Otq2VkiRpcFmJkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKK1ZKQyIijqZYbGc/io3DfgScnZmX1ZwzCbwyMy9pUxs+BTwiM49o8vwnAI/KzP+Yw8/8JPDYzDyswX2PoVhy+pLM/Nu6+/6YYvW+/y8z/7vsz5cGmZUIaQhExLEUW/yeB/wB8IfAfwCfjYjarX73olgCt1d8CXhqux48M2+hWJ79dRHxwCZqlX1QPg182AAhTc1KhDQcXgd8IjM/VXPspxERwFsoPjCp2QelV4zMfMrcZOYFEfEnwIUR8URgDfA54BfAe9r986V+ZoiQhsNW4JkRsXNm1m6B/DZgYfWb2u6MStfDFootk4+pPMaZwP8DLgCeDCRwTGb+T/31jR6zvlER8RfAO4H9gUmK/Q3empk3RMRy4DHASRHx6szcJyLmAx+g2ERpYeX8v8/Ma2se8w2V57VHpa3NVFyPodj6+Bzg5kp7DsjM+5u4VhpadmdIw+EjwNOA30XEv0fE2yLigMy8OzNvnea6oyu3TwE+SrGp2L8BH6w83mbgY2UaFBFPBS4HPgXsS7EN9gjwicopLwRupdgMqNqlcTHF9vIvptiq+FvAlRHx+MpjvpJiN9kPUGy69FvgZTO1pbJR07HASymqD6/LzBVlnpc0TAwR0hDIzCuAZ1GMgziUIlT8ICK+HxH7TXPp3cDbK2MHPlo59q+Z+ZXMvIlix9L9SzZrC/B3mfmxzLw1M2+gCBBPrLR5jKL6sT4z746Ix1KEh1dn5ncy8+eZeQpwNcXOsVAMHL0kMz+RhXcCNzTZnquAe4D7gWtKPidpqNidIQ2JzLwGuCYitqOoLPwJ8CbgvyLisZm5ucFlt1S2OiczNxRDKLil5v77gPkl2/PDiFgTEe+i2Db+ccABTP3HzYGV2+sq7aiaX9OG/amM76hxLfCkJpr0CWA1xe/FiyPiOZk50cR10tAyREgDrjLT4F3A+zJzZWZuBa4Hro+I7wBfo/iQ/Z8Gl29pcKzpD9aImPJ3TEQcDvwXxQyM7wL/AjweOH+KS6oh5xkU4aXWeOV2kocOxmwUjurb8jrgBcBzKH4vfoOiuvGRma6VhpndGdLgu49i4ODLG9y3huKD964W/awtwJKa7x83zbl/B3w9M1+SmWdn5pXAPgARUQ0CkzXn/6Ryu0dm/rL6DzgO+LPKfT8EDq77OQdN1+CI2Jeiq+a0zLwqM79JMcDy1Ij4g+mulYadlQhpwGXmqoj4MPDBiFgCfIEiWDwROBX4dGb+pkU/7nvAsRHxXWA7ig/n8SnOvRt4fkQ8HbgT+GPgrZX75gObgHXA4yPi4Zn5y4i4DLigMgPj58BfA38LHFm57jTg8xFxPUWV48UUY0GubtSAiNiRYjpnsu10zndWHvOSiDgoM6d6DtJQsxIhDYHMPJFirYg/ovhA/QnwjxTjB45t4Y96PcWU0OsoFq26gGKGRCPvpZii+TXgfylmY1QXvqrOxjgDOAr4UUSMUlRU/pNiQOePK/e9sFI9IDP/DXg18AaKFTmfXmnDVD5C0YXyitoxIZl5H/BK4P9QzESR1MDI5OTkzGdJkiTVsRIhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKuX/B6Crc36eRUyVAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x, y = sims.linear_sim(num_samp=100, num_dim=1, noise=0.1)\n", "\n", "fig = plt.figure(figsize=(8,8))\n", "fig.suptitle(\"Linear Simulation\", fontsize=17)\n", "ax = sns.scatterplot(x=x[:,0], y=y[:,0])\n", "ax.set_xlabel('Simulated X', fontsize=15)\n", "ax.set_ylabel('Simulated Y', fontsize=15) \n", "plt.axis('equal')\n", "plt.xticks(fontsize=15)\n", "plt.yticks(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The test statistic and p-value can be called by creating the `DCorr` object and simply calling the corresponding test statistic and p-value methods. The value of `which_test` can be modified to 3 different values: `biased` for biased Dcorr, `unbiased` for unbiased Dcorr, and `mantel` for Mantel." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dcorr test statistic: 0.9591913132744958\n", "P Value: 0.001\n" ] } ], "source": [ "dcorr = DCorr(which_test='mantel')\n", "dcorr_statistic, independence_test_metadata = dcorr.test_statistic(x, y)\n", "p_value, _ = dcorr.p_value(x, y)\n", "\n", "print(\"Dcorr test statistic:\", dcorr_statistic)\n", "print(\"P Value:\", p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The p-value is bounded by the number of repetitions (in this case 1000). This is because since we are estimating the null distribution via permutation, this is the lowest value that we can be sufficiently sure is the p-value. It is worth noting that as in most of the other tests that use permutation to approximate the p-value, the `replication_factor` parameter can be set to the desired number." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }