Fig2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8b748901",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  from pandas import Int64Index as NumericIndex\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "import numpy as np\n",
    "import scipy as sc\n",
    "from itertools import product\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import PIL\n",
    "from numpy import log10\n",
    "import random\n",
    "from math import factorial\n",
    "from scipy.stats import linregress, gaussian_kde, skew\n",
    "from scipy import stats\n",
    "from scipy.spatial import distance\n",
    "import warnings\n",
    "import pandas as pd\n",
    "import re\n",
    "import os\n",
    "import math\n",
    "from collections import Counter\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from statsmodels.stats.outliers_influence import summary_table\n",
    "import statsmodels.api as sm\n",
    "from scipy.optimize import linear_sum_assignment\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "%config InlineBackend.print_figure_kwargs={'facecolor' : \"w\"}\n",
    "\n",
    "pd.set_option('display.max_columns', None)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e696cc30",
   "metadata": {},
   "source": [
    "# Figure 2.\n",
    "\n",
    "## Analyzing species rarity\n",
    " \n",
    "The analysis of shift is highly relevant to the study of species-abundance distributions (SADs) (Fig 8a) [10, 12, 14, 21]. These histograms of species abundance underpin thousands of ecological studies spanning all domains of life and major habitats [10, 12, 14, 21]. SADs are often compared to theoretical predictions, to SADs sampled from other communities of similar taxa, and are the basis for many measures of species dominance, diversity, evenness, and rarity [9, 10, 12, 14, 21]. Consequently, the field of ecology is replete with techniques for analyzing SADs and for quantifying the aspects of biodiversity they contain.\n",
    "\n",
    "\n",
    "In considering that RDS provides a means of quantifying the magnitude and direction by which one SAD is concentrated to lesser or greater abundances relative to another SAD,  \n",
    "\n",
    "Beyond providing yet another means of comparing SADs, DS and RDS provide an intuitive and easily interpretable means of quantifying species rarity, i.e., the concentration of species at low abundances [14, 21]. \n",
    "\n",
    "However, \n",
    "\n",
    "As the concentration of frequencies away from the discrete class having the greatest value, DS is highly similar in concept to species rarity. Consequently, given its bounded values and intuitive interpretation, we asked whether DS provides a preferrable measure of species rarity. in particular, because skewness based measures are not bounded an\n",
    "\n",
    "Unlike skewness-based measures, RDS is bounded, intuitive, and directly reflects the concentration of species towards the class of lowest abundance. Likewise, RDS represents a better comparative metric of species rarity for the same reasons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32ce696d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_pts_within_radius(x, y, radius, scale=0):\n",
    "    \"\"\"Count the number of points within a fixed radius in 2D space\"\"\"\n",
    "    \n",
    "    raw_data = np.array([x, y])\n",
    "    x = np.array(x)\n",
    "    y = np.array(y)\n",
    "    raw_data = raw_data.transpose()\n",
    "    \n",
    "    # Get unique data points by adding each pair of points to a set\n",
    "    unique_points = set()\n",
    "    for xval, yval in raw_data:\n",
    "        unique_points.add((xval, yval))\n",
    "    \n",
    "    count_data = []\n",
    "    for a, b in unique_points:\n",
    "        if scale == 'sqrt':\n",
    "            num_neighbors = len(x[((sqrt(x) - sqrt(a)) ** 2 +\n",
    "                                   (sqrt(y) - sqrt(b)) ** 2) <= sqrt(radius) ** 2])\n",
    "        else:        \n",
    "            num_neighbors = len(x[((x - a) ** 2 + (y - b) ** 2) <= radius ** 2])\n",
    "        count_data.append((a, b, num_neighbors))\n",
    "    return count_data\n",
    "\n",
    "\n",
    "\n",
    "def plot_color_by_pt_dens(x, y, radius, loglog=0, scale=0, plot_obj=None, point_size=10):\n",
    "    \n",
    "    \"\"\"Plot bivariate relationships with large n using color for point density\n",
    "\n",
    "    Inputs:\n",
    "    x & y -- variables to be plotted\n",
    "    radius -- the linear distance within which to count points as neighbors\n",
    "    scale -- a flag to indicate the use of a scale plot (scale = 1)\n",
    "\n",
    "    The color of each point in the plot is determined by the logarithm (base 10)\n",
    "    of the number of points that occur with a given radius of the focal point,\n",
    "    with hotter colors indicating more points. The number of neighboring points\n",
    "    is determined in linear space regardless of whether a scale plot is\n",
    "    presented.\n",
    "    \"\"\"\n",
    "    \n",
    "    plot_data = count_pts_within_radius(x, y, radius, scale)\n",
    "    sorted_plot_data = np.array(sorted(plot_data, key=lambda point: point[2]))\n",
    "\n",
    "    if plot_obj == None:\n",
    "        plot_obj = plt.axes()\n",
    "        \n",
    "    plot_obj.scatter(sorted_plot_data[:, 0],\n",
    "            sorted_plot_data[:, 1],\n",
    "            facecolors='none',\n",
    "            s = point_size, \n",
    "            edgecolors='0.1', \n",
    "            linewidths=1., \n",
    "            #cmap='Greys'\n",
    "            )\n",
    "    \n",
    "    # plot points\n",
    "    c = np.array(sorted_plot_data[:, 2])**0.25\n",
    "    c = np.max(c) - c\n",
    "    plot_obj.scatter(sorted_plot_data[:, 0],\n",
    "                    sorted_plot_data[:, 1],\n",
    "                    c = c,\n",
    "                    s = point_size, \n",
    "                    edgecolors='k', \n",
    "                    linewidths=0.0, \n",
    "                    cmap='Greys_r',\n",
    "                    #alpha = 0.5,\n",
    "                    )\n",
    "        \n",
    "    return plot_obj\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f13e90fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def DS(p):\n",
    "    p_bins = len(p)\n",
    "    p_obs = sum(p)\n",
    "    \n",
    "    z_p = (p_bins + 1)/p_bins\n",
    "    p = [sum(p[:ii+1])**(z_p) for ii in range(len(p))]\n",
    "    Sp = np.sum(np.array(p)/(p_obs**z_p)) - 1\n",
    "    ds = Sp/(p_bins - 1)\n",
    "    return ds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8deb13da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Rlogskew(sad):\n",
    "    '''\n",
    "    Calculation of rarity used in:\n",
    "        A.E. Magurran and B.J. McGill, eds. Biological diversity: frontiers in measurement and assessment. \n",
    "        OUP Oxford, 2010.\n",
    "\n",
    "    '''\n",
    "    \n",
    "    S = len(sad)\n",
    "\n",
    "    if S <= 2.0:\n",
    "        print('S < 2, cannot compute log-skew')\n",
    "        sys.exit()\n",
    "\n",
    "    sad = np.log10(sad)\n",
    "    mu = np.mean(sad)\n",
    "\n",
    "    num = 0\n",
    "    denom = 0\n",
    "    for ni in sad:\n",
    "        num += ((ni - mu)**3.0)/S\n",
    "        denom += ((ni - mu)**2.0)/S\n",
    "\n",
    "    t1 = num/(denom**(3.0/2.0))\n",
    "    t2 = (S/(S - 2.0)) * np.sqrt((S - 1.0)/S)\n",
    "\n",
    "    return t1 * t2\n",
    "\n",
    "\n",
    "def get_RADs(path, name, closedref=True):\n",
    "\n",
    "    # Get rank-abundance distributions, i.e, abundance vectors\n",
    "    \n",
    "    RADdict = {}\n",
    "    DATA = path + name + '-data.txt'\n",
    "\n",
    "    with open(DATA) as f:\n",
    "\n",
    "        for d in f:\n",
    "\n",
    "            if d.strip():\n",
    "                d = d.split()\n",
    "                length = len(d)\n",
    "\n",
    "                if name == 'GENTRY':\n",
    "                    site = d[0]\n",
    "                    #species = d[1] # Dataset name plus species identifier\n",
    "                    abundance = float(d[-1])\n",
    "\n",
    "                else:\n",
    "                    site = d[0]\n",
    "                    #year = d[1]\n",
    "\n",
    "                    if closedref == True:\n",
    "                        for i in d:\n",
    "                            if 'unclassified' in i:\n",
    "                                #print('unclassified')\n",
    "                                continue\n",
    "                            elif 'unidentified' in i:\n",
    "                                #print('unidentified')\n",
    "                                continue\n",
    "\n",
    "                    abundance = float(d[-1])\n",
    "\n",
    "\n",
    "                if abundance > 0:\n",
    "                    if site in RADdict:\n",
    "                        RADdict[site].append(abundance)\n",
    "                    else:\n",
    "                        RADdict[site] = [abundance]\n",
    "\n",
    "    RADs = RADdict.values()\n",
    "    filteredRADs = []\n",
    "    for rad in RADs:\n",
    "        if len(rad) >= 10:\n",
    "            filteredRADs.append(rad)\n",
    "\n",
    "    return filteredRADs\n",
    "\n",
    "\n",
    "\n",
    "def EMP_RADs(path, name):\n",
    "\n",
    "    minS = 10\n",
    "\n",
    "    IN = path + '/EMPclosed-SADs.txt'\n",
    "    rads = []\n",
    "    \n",
    "    with open(IN) as f:\n",
    "\n",
    "        for rad in f:\n",
    "            rad = eval(rad)\n",
    "            if len(rad) >= minS:\n",
    "                rads.append(rad)\n",
    "\n",
    "    return rads\n",
    "\n",
    "def Louca_RADs(path, name):\n",
    "    \n",
    "    # Get rank-abundance distributions, i.e, abundance vectors\n",
    "    \n",
    "    RADdict = {}\n",
    "    DATA = path + 'SSADdata.txt'\n",
    "\n",
    "    with open(DATA) as f:\n",
    "\n",
    "        for d in f:\n",
    "\n",
    "            if d.strip():\n",
    "                d = d.split()\n",
    "                length = len(d)\n",
    "\n",
    "                site = d[1]\n",
    "                for i in d:\n",
    "                    abundance = float(d[-1])\n",
    "\n",
    "\n",
    "                if abundance > 0:\n",
    "                    if site in RADdict:\n",
    "                        RADdict[site].append(abundance)\n",
    "                    else:\n",
    "                        RADdict[site] = [abundance]\n",
    "\n",
    "    RADs = RADdict.values()\n",
    "    filteredRADs = []\n",
    "    for rad in RADs:\n",
    "        if len(rad) >= 10:\n",
    "            filteredRADs.append(rad)\n",
    "\n",
    "    return filteredRADs\n",
    "\n",
    "\n",
    "def NSECF(p):\n",
    "    p_bins = len(p)\n",
    "    p_obs = sum(p)\n",
    "    \n",
    "    z_p = (p_bins + 1)/p_bins\n",
    "    p = [sum(p[:ii+1])**(z_p) for ii in range(len(p))]\n",
    "    return np.sum(np.array(p)/(p_obs**z_p)) - 1\n",
    "    \n",
    "                \n",
    "def getMetrics():\n",
    "    \n",
    "    name_ls = []\n",
    "    kind_ls = []\n",
    "    N_ls = []\n",
    "    S_ls = []\n",
    "    skew_ls = []\n",
    "    logskew_ls = []\n",
    "    log_mod_skew_ls = []\n",
    "    log_mod_skew_log_ls = []\n",
    "    ds_ls = []\n",
    "    nsecf_ls = []\n",
    "    \n",
    "    datasets = []\n",
    "    for name in os.listdir('data/ecological/micro'):\n",
    "            datasets.append([name, 'micro'])\n",
    "    for name in os.listdir('data/ecological/macro'):\n",
    "            datasets.append([name, 'macro'])\n",
    "\n",
    "    for dataset in datasets:\n",
    "        \n",
    "        name = dataset[0] # name of dataset\n",
    "        kind = dataset[1] # micro or macro\n",
    "\n",
    "        if name == '.DS_Store' or name == 'MGRAST': \n",
    "            continue\n",
    "        \n",
    "        OUT = open('data/ecological/'+kind+'/'+name+'/'+name+'-SADMetricData.txt','w+')\n",
    "        RADs = []\n",
    "\n",
    "        if kind == 'macro':\n",
    "            RADs = get_RADs('data/ecological/'+kind+'/'+name+'/', name)\n",
    "            print('macro', name, len(RADs))\n",
    "\n",
    "\n",
    "        if kind == 'micro':\n",
    "            if name == 'EMPclosed' or name == 'EMPopen':\n",
    "                RADs = EMP_RADs('data/ecological/'+kind+'/'+name+'/', name)\n",
    "            \n",
    "            elif name == 'Louca':\n",
    "                RADs = Louca_RADs('data/ecological/'+kind+'/'+name+'/', name)\n",
    "            \n",
    "            else:\n",
    "                RADs = get_RADs('data/ecological/'+kind+'/'+name+'/', name)\n",
    "\n",
    "            print('micro', name, len(RADs))\n",
    "\n",
    "        ct = 0\n",
    "        numRADs = len(RADs)\n",
    "        for RAD in RADs:\n",
    "\n",
    "            if kind == 'micro':\n",
    "                RAD = list([x for x in RAD if x > 0])\n",
    "\n",
    "            elif kind == 'macro':\n",
    "                RAD = list([x for x in RAD if x > 0])\n",
    "\n",
    "\n",
    "            N = sum(RAD)\n",
    "            S = len(RAD)\n",
    "\n",
    "            if S < 10: \n",
    "                continue\n",
    "            if max(RAD) == min(RAD): \n",
    "                continue\n",
    "\n",
    "            # Measures of Rarity\n",
    "            \n",
    "            # 1. skewness of abundances\n",
    "            skewness = skew(RAD)\n",
    "            \n",
    "            # 2. log-modulo transformation of skewnness\n",
    "            lms = np.log10(np.abs(float(skewness)) + 1)\n",
    "            if skewness < 0: \n",
    "                lms = lms * -1\n",
    "            log_mod_skew = float(lms)\n",
    "            \n",
    "            # 3. skewness of log-transformed abundances\n",
    "            logskew = Rlogskew(RAD)\n",
    "            \n",
    "            # 4. log-modulo transformation of skewnness of log-transformed abundances\n",
    "            lms = np.log10(np.abs(float(logskew)) + 1)\n",
    "            if skewness < 0: \n",
    "                lms = lms * -1\n",
    "            log_mod_skew_log = float(lms)\n",
    "            \n",
    "            # 5. Distributional shift (DS)\n",
    "            # Convert the abundances to logarithmic scale (base 2)\n",
    "            abundances = np.log2(RAD)\n",
    "\n",
    "            # Define the bins for the histogram\n",
    "            min_abundance = 0 #np.floor(min(abundances))\n",
    "            max_abundance = np.ceil(max(abundances))\n",
    "            bins = np.arange(min_abundance, max_abundance + 1, 1)\n",
    "\n",
    "            # Compute the histogram\n",
    "            hist, bin_edges = np.histogram(abundances, bins=bins)\n",
    "\n",
    "            # Use the right side of the bin edges as bin values\n",
    "            bin_values = bin_edges[1:]\n",
    "\n",
    "            # Convert histogram to list\n",
    "            bin_heights = hist.tolist()\n",
    "            ds = DS(bin_heights)\n",
    "            \n",
    "            # 6. Normalized sums of exponentiated cumulative frequencies\n",
    "            nsecf = NSECF(bin_heights)\n",
    "\n",
    "            ct+=1\n",
    "\n",
    "            print(name, kind, N, S, skewness, logskew, log_mod_skew, log_mod_skew_log, ds, nsecf, file=OUT)\n",
    "            \n",
    "            name_ls.append(name)\n",
    "            kind_ls.append(kind)\n",
    "            N_ls.append(N)\n",
    "            S_ls.append(S)\n",
    "            skew_ls.append(skewness)\n",
    "            logskew_ls.append(logskew)\n",
    "            log_mod_skew_ls.append(log_mod_skew)\n",
    "            log_mod_skew_log_ls.append(log_mod_skew_log)\n",
    "            ds_ls.append(ds)\n",
    "            nsecf_ls.append(nsecf)\n",
    "            \n",
    "        OUT.close()\n",
    "        \n",
    "    return name_ls, kind_ls, N_ls, S_ls, skew_ls, logskew_ls, log_mod_skew_ls, log_mod_skew_log_ls, ds_ls, nsecf_ls\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7874d9a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "micro CHU 44\n",
      "micro CATLIN 130\n",
      "micro BOVINE 16\n",
      "micro LAUB 84\n",
      "micro SED 34\n",
      "micro BIGN 13\n",
      "micro CHINA 186\n",
      "micro TARA 139\n",
      "micro HMP 4303\n",
      "micro FUNGI 128\n",
      "micro Louca 8308\n",
      "micro HUMAN 525\n",
      "micro HYDRO 123\n",
      "micro EMPclosed 14631\n",
      "macro FIA 10355\n",
      "macro CBC 1999\n",
      "macro MCDB 103\n",
      "macro BBS 2769\n",
      "macro GENTRY 222\n",
      "(44090, 10)\n",
      "['CHU', 'CATLIN', 'BOVINE', 'LAUB', 'SED', 'BIGN', 'CHINA', 'TARA', 'HMP', 'FUNGI', 'Louca', 'HUMAN', 'HYDRO', 'EMPclosed', 'FIA', 'CBC', 'MCDB', 'BBS', 'GENTRY']\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>kind</th>\n",
       "      <th>N</th>\n",
       "      <th>S</th>\n",
       "      <th>skew</th>\n",
       "      <th>logskew</th>\n",
       "      <th>log_mod_skew</th>\n",
       "      <th>log_mod_skew_log</th>\n",
       "      <th>DS</th>\n",
       "      <th>NSECF</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHU</td>\n",
       "      <td>micro</td>\n",
       "      <td>2487.0</td>\n",
       "      <td>350</td>\n",
       "      <td>11.376136</td>\n",
       "      <td>1.415601</td>\n",
       "      <td>1.092585</td>\n",
       "      <td>0.383025</td>\n",
       "      <td>0.839037</td>\n",
       "      <td>6.712296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CHU</td>\n",
       "      <td>micro</td>\n",
       "      <td>2032.0</td>\n",
       "      <td>306</td>\n",
       "      <td>9.868336</td>\n",
       "      <td>1.355590</td>\n",
       "      <td>1.036163</td>\n",
       "      <td>0.372100</td>\n",
       "      <td>0.832462</td>\n",
       "      <td>6.659695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CHU</td>\n",
       "      <td>micro</td>\n",
       "      <td>872.0</td>\n",
       "      <td>181</td>\n",
       "      <td>4.654126</td>\n",
       "      <td>1.179171</td>\n",
       "      <td>0.752365</td>\n",
       "      <td>0.338291</td>\n",
       "      <td>0.801513</td>\n",
       "      <td>4.809080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CHU</td>\n",
       "      <td>micro</td>\n",
       "      <td>1241.0</td>\n",
       "      <td>202</td>\n",
       "      <td>9.932251</td>\n",
       "      <td>1.676968</td>\n",
       "      <td>1.038710</td>\n",
       "      <td>0.427643</td>\n",
       "      <td>0.861797</td>\n",
       "      <td>6.894373</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CHU</td>\n",
       "      <td>micro</td>\n",
       "      <td>785.0</td>\n",
       "      <td>174</td>\n",
       "      <td>7.223195</td>\n",
       "      <td>1.498230</td>\n",
       "      <td>0.915041</td>\n",
       "      <td>0.397632</td>\n",
       "      <td>0.815210</td>\n",
       "      <td>4.891259</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  name   kind       N    S       skew   logskew  log_mod_skew  \\\n",
       "0  CHU  micro  2487.0  350  11.376136  1.415601      1.092585   \n",
       "1  CHU  micro  2032.0  306   9.868336  1.355590      1.036163   \n",
       "2  CHU  micro   872.0  181   4.654126  1.179171      0.752365   \n",
       "3  CHU  micro  1241.0  202   9.932251  1.676968      1.038710   \n",
       "4  CHU  micro   785.0  174   7.223195  1.498230      0.915041   \n",
       "\n",
       "   log_mod_skew_log        DS     NSECF  \n",
       "0          0.383025  0.839037  6.712296  \n",
       "1          0.372100  0.832462  6.659695  \n",
       "2          0.338291  0.801513  4.809080  \n",
       "3          0.427643  0.861797  6.894373  \n",
       "4          0.397632  0.815210  4.891259  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_ls, kind_ls, N_ls, S_ls, skew_ls, logskew_ls, log_mod_skew_ls, log_mod_skew_log_ls, ds_ls, nsecf_ls = getMetrics()\n",
    "\n",
    "main_df = pd.DataFrame({\n",
    "    'name': name_ls,\n",
    "    'kind': kind_ls,\n",
    "    'N': N_ls,\n",
    "    'S': S_ls,\n",
    "    'skew': skew_ls,\n",
    "    'logskew': logskew_ls,\n",
    "    'log_mod_skew': log_mod_skew_ls,\n",
    "    'log_mod_skew_log': log_mod_skew_log_ls,\n",
    "    'DS': ds_ls,\n",
    "    'NSECF': nsecf_ls,\n",
    "})\n",
    "\n",
    "print(main_df.shape)\n",
    "\n",
    "print(main_df['name'].unique().tolist())\n",
    "main_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d251120e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(44090, 10)\n"
     ]
    }
   ],
   "source": [
    "print(main_df.shape)\n",
    "\n",
    "its = 20\n",
    "tdf = None\n",
    "names = ['CHU', 'CATLIN', 'BOVINE', 'LAUB', 'SED', 'BIGN', \n",
    "         'CHINA', 'TARA', 'HMP', 'FUNGI', 'Louca', 'HUMAN', \n",
    "         'HYDRO', 'EMPclosed', 'FIA', 'CBC', 'MCDB', 'BBS', \n",
    "         'GENTRY']\n",
    "\n",
    "for n in range(its):\n",
    "\n",
    "    for i, name in enumerate(names):\n",
    "                \n",
    "        tdf_nm = main_df[main_df['name'] == name]\n",
    "        kind = tdf_nm['kind'].iloc[0]\n",
    "        numlines = tdf_nm.shape[0]\n",
    "                \n",
    "        small = ['BIGN', 'BOVINE', 'CHU', 'LAUB', 'SED']\n",
    "        big = ['HUMAN', 'CHINA', 'CATLIN', 'FUNGI', 'HYDRO']\n",
    "\n",
    "        if name == 'Louca':\n",
    "            tdf_nm = tdf_nm.sample(1000, replace=True)\n",
    "\n",
    "        elif kind == 'macro':\n",
    "            tdf_nm = tdf_nm.sample(100, replace=True)\n",
    "        elif name in small:\n",
    "            tdf_nm = tdf_nm.sample(20, replace=True)\n",
    "        elif name in big:\n",
    "            tdf_nm = tdf_nm.sample(50, replace=True)\n",
    "        elif name == 'TARA':\n",
    "            tdf_nm = tdf_nm.sample(50, replace=True)\n",
    "        else:\n",
    "            tdf_nm = tdf_nm.sample(50, replace=True)\n",
    "            \n",
    "        if n == 0 and i == 0:\n",
    "            tdf = tdf_nm.copy(deep=True)\n",
    "        else:\n",
    "            tdf = pd.concat([tdf, tdf_nm])\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b6c182bc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regression for log_mod_skew:\n",
      "Slope: 0.18601026595889422\n",
      "Intercept: 0.13896086999748047\n",
      "r2: 0.5922371138946283 \n",
      "\n",
      "Regression for DS:\n",
      "Slope: 0.027741989515227494\n",
      "Intercept: 0.602322091588741\n",
      "r2: 0.06779993014824855 \n",
      "\n",
      "Regression for NSECF:\n",
      "Slope: 2.405479534031236\n",
      "Intercept: -1.8782641482141953\n",
      "r2: 0.8477463486565046\n"
     ]
    }
   ],
   "source": [
    "\n",
    "tdf.dropna(subset=['DS'], how='any', inplace=True)\n",
    "\n",
    "# Log scale the 'N' values\n",
    "log_N = np.log10(tdf['N'])\n",
    "\n",
    "# log_mod_skew values\n",
    "log_mod_skew = tdf['log_mod_skew']\n",
    "\n",
    "# DS values\n",
    "DS_vals = tdf['DS']\n",
    "\n",
    "# NSECF\n",
    "NSECF_vals = tdf['NSECF']\n",
    "\n",
    "# Reshape the data to fit the model\n",
    "X = log_N.values\n",
    "y1 = log_mod_skew.values\n",
    "y2 = DS_vals.values\n",
    "y3 = NSECF_vals.values\n",
    "\n",
    "# Fit the linear regression model for log_mod_skew\n",
    "slope1, intercept1, r_value1, p1, se1 = linregress(X, y1)\n",
    "print(\"Regression for log_mod_skew:\")\n",
    "print(\"Slope:\", slope1)\n",
    "print(\"Intercept:\", intercept1)\n",
    "print(\"r2:\", r_value1 ** 2, '\\n')\n",
    "\n",
    "# Fit the linear regression model for DS\n",
    "slope2, intercept2, r_value2, p2, se2 = linregress(X, y2)\n",
    "print(\"Regression for DS:\")\n",
    "print(\"Slope:\", slope2)\n",
    "print(\"Intercept:\", intercept2)\n",
    "print(\"r2:\", r_value2 ** 2, '\\n')\n",
    "\n",
    "# Fit the linear regression model for NSECF\n",
    "slope3, intercept3, r_value3, p3, se3 = linregress(X, y3)\n",
    "print(\"Regression for NSECF:\")\n",
    "print(\"Slope:\", slope3)\n",
    "print(\"Intercept:\", intercept3)\n",
    "print(\"r2:\", r_value3 ** 2)\n",
    "\n",
    "# Plotting\n",
    "fig = plt.figure(figsize=(10, 2.5))\n",
    "x_lab = 'log(N)'\n",
    "\n",
    "# Scatter plot and regression line for log_mod_skew\n",
    "plot_color_by_pt_dens(X, y1, radius=0.05, loglog=0, plot_obj=plt.subplot(1, 3, 1), point_size=10)\n",
    "slope, intercept, r_val, p_val, std_err = linregress(X, y1)\n",
    "fitted_vals = slope * np.array(X) + intercept\n",
    "s = r'$r_{2}$' + ' = ' + str(round(r_val**2, 2))\n",
    "plt.plot(X, fitted_vals, color='k', linewidth=2, label=s)\n",
    "plt.xlabel(x_lab, fontweight='bold')\n",
    "plt.ylabel('log-modulo skewness', fontweight='bold')\n",
    "plt.title(r'Rarity = '+str(round(10**intercept,2))+'*'+r'$N$'+'$^{'+str(round(slope,2))+'}$')\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "# Scatter plot and regression line for DS (non-standardized)\n",
    "plot_color_by_pt_dens(X, y3, radius=0.05, loglog=0, plot_obj=plt.subplot(1, 3, 2), point_size=10)\n",
    "slope, intercept, r_val, p_val, std_err = linregress(X, y3)\n",
    "fitted_vals = slope * np.array(X) + intercept\n",
    "s = r'$r_{2}$' + ' = ' + str(round(r_val**2, 2))\n",
    "plt.plot(X, fitted_vals, color='k', linewidth=2, label=s)\n",
    "plt.xlabel(x_lab, fontweight='bold')\n",
    "plt.ylabel('DS, non-standardized', fontweight='bold')\n",
    "plt.title(r'Rarity = '+str(round(slope,2))+'*'+r'$N$'+ ' - ' + str(abs(round(intercept,2))))\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "\n",
    "# Scatter plot and regression line for DS\n",
    "plot_color_by_pt_dens(X, y2, radius=0.05, loglog=0, plot_obj=plt.subplot(1, 3, 3), point_size=10)\n",
    "slope, intercept, r_val, p_val, std_err = linregress(X, y2)\n",
    "fitted_vals = slope * np.array(X) + intercept\n",
    "s = r'$r_{2}$' + ' = ' + str(round(r_val**2, 2))\n",
    "plt.plot(X, fitted_vals, color='k', linewidth=2, label=s)\n",
    "plt.xlabel(x_lab, fontweight='bold')\n",
    "plt.ylabel('DS, standardized', fontweight='bold')\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "\n",
    "fig.patch.set_facecolor('white')\n",
    "plt.subplots_adjust(wspace=0.45, hspace=0.35)\n",
    "plt.savefig('Final_Figs/manuscript/Fig2.pdf', bbox_inches='tight', format='pdf', dpi=600)\n",
    "plt.savefig('Final_Figs/manuscript/Fig2.jpg', bbox_inches='tight', format='jpg', dpi=600)\n",
    "plt.close()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}