03_Missing_Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Missing Data and Imputation\n",
    "Compiled by [Morgan Williams](mailto:morgan.williams@csiro.au) for C3DIS 2018 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 'achilles heel' of compositional data analysis is the incompatibility with 'null' or zero components. In practice, most compositional data contains such values - as either true zeros (e.g. count data as found in surveys), below-detection values (e.g. geochemistry) or components which are simply missing (e.g. components are not measured, but may be of significant quantity).\n",
    "\n",
    "The simplest but ultimately least practical solution to this problem is to use subcompositions - a set of components free of missing/zero data. Beyond this approach, missing data may be able to be imputed using relationships between variables. In the absence of additional information, typical data imputation methods utilise some measure of central tendency (e.g. median, mode, mean) as a fill value for missing data; this preserves the expected values, and largely preserves the covariance stucture of the dataset. However, we can use information we already have to provide more accurate imputation methods - leading to more robust analysis workflows.\n",
    "\n",
    "For compositional data, missing values commonly represent a 'detection limit' - a specific value below which the variable cannot be quantified, nor certified to be above background/zero (e.g. commonly 2$\\sigma$ above zero). With regards to this form of missing data, a measure of central tendency is not likely an accurate representation - instead imputed values should lie below some threshold corresponding to the 'detction limit'. Here, the simplest form of imputation uses an arbitrarily small fill value. Imputation using nominal values may be marginally valid in some cases (i.e. using small values to represent values below decection, without altering the closure operation greatly), but overall this approach typically serves as a confounding factor (e.g. creating bimodal distributions and spurious clusters).\n",
    "\n",
    "Regardless of the method of imputation, when values are imputed from low density dataset, the output is strongly dependent on the data quality and 'representiveness' of the present values. For this reason, using imputed values for geological inference may be misleading - especially for rarely recorded parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.insert(0, './src')\n",
    "from geochem import *\n",
    "from compositions import *\n",
    "from EMimputation import EMCOMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simple_impute(arr, method='mean', threshold=0.01, limits=[]):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    if method == 'mean':\n",
    "        # Note this is not a compositional mean, which should be: inv_alr(np.nanmean(alr(arr), axis=0)); this would require 1 nan-free column\n",
    "        replace_vals = np.nanmean(arr, axis=0)\n",
    "    elif method == 'eps': # constant arbitrarily small value\n",
    "        replace_vals = 2*np.finfo(float).eps * np.ones(arr.shape)\n",
    "    elif method == 'threshold': # detection limits or some ratio of the minimum values\n",
    "        if not len(limits):\n",
    "            replace_vals = np.tile(np.nanmin(arr, axis=0), (arr.shape[0])).reshape(arr.shape)\n",
    "        else:\n",
    "            assert len(limits) == arr.shape[1]\n",
    "            replace_vals = np.tile(np.array(limits), (arr.shape[0])).reshape(arr.shape)\n",
    "    else:\n",
    "        raise NotImplementedError\n",
    "        \n",
    "    return np.where(~np.isfinite(arr), replace_vals, arr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>Ti</th>\n",
       "      <th>Ge</th>\n",
       "      <th>Tm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>45.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   SiO2  Al2O3   MgO   CaO    Ti    Ge    Tm\n",
       "0  45.0   10.0  20.0  25.0   NaN   NaN   NaN\n",
       "1  40.0   12.0  23.0  25.0   NaN   NaN   NaN\n",
       "2  10.0    NaN   NaN   NaN  10.0  20.0  25.0\n",
       "3  10.0    3.0   NaN   NaN  10.0  20.0  25.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r1 = pd.DataFrame(dict(SiO2=[45.0], Al2O3=[10.0], MgO=[20.0], CaO=[25.0]),)\n",
    "r2 = pd.DataFrame(dict(SiO2=[40.0], Al2O3=[12.0], MgO=[23.0], CaO=[25.0]),)\n",
    "r3 = pd.DataFrame(dict(SiO2=[10.0], Ti=[10.0], Ge=[20.0], Tm=[25.0]),)\n",
    "r4 = pd.DataFrame(dict(SiO2=[10.0], Al2O3=[3.0], Ti=[10.0], Ge=[20.0], Tm=[25.0]),)\n",
    "aggdf = r1\n",
    "for r in [r2, r3, r4]:\n",
    "    aggdf = aggdf.append(r, sort=False)\n",
    "    aggdf = aggdf.reset_index(drop=True)\n",
    "aggdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>Ti</th>\n",
       "      <th>Ge</th>\n",
       "      <th>Tm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>45.0</td>\n",
       "      <td>10.000000</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40.0</td>\n",
       "      <td>12.000000</td>\n",
       "      <td>23.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10.0</td>\n",
       "      <td>8.333333</td>\n",
       "      <td>21.5</td>\n",
       "      <td>25.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10.0</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>21.5</td>\n",
       "      <td>25.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   SiO2      Al2O3   MgO   CaO    Ti    Ge    Tm\n",
       "0  45.0  10.000000  20.0  25.0  10.0  20.0  25.0\n",
       "1  40.0  12.000000  23.0  25.0  10.0  20.0  25.0\n",
       "2  10.0   8.333333  21.5  25.0  10.0  20.0  25.0\n",
       "3  10.0   3.000000  21.5  25.0  10.0  20.0  25.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_imputed = aggdf.copy()\n",
    "mean_imputed.loc[:, :] = simple_impute(aggdf.values, method='mean')\n",
    "mean_imputed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>Ti</th>\n",
       "      <th>Ge</th>\n",
       "      <th>Tm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>45.0</td>\n",
       "      <td>1.000000e+01</td>\n",
       "      <td>2.000000e+01</td>\n",
       "      <td>2.500000e+01</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40.0</td>\n",
       "      <td>1.200000e+01</td>\n",
       "      <td>2.300000e+01</td>\n",
       "      <td>2.500000e+01</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10.0</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>1.000000e+01</td>\n",
       "      <td>2.000000e+01</td>\n",
       "      <td>2.500000e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10.0</td>\n",
       "      <td>3.000000e+00</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>4.440892e-16</td>\n",
       "      <td>1.000000e+01</td>\n",
       "      <td>2.000000e+01</td>\n",
       "      <td>2.500000e+01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   SiO2         Al2O3           MgO           CaO            Ti            Ge  \\\n",
       "0  45.0  1.000000e+01  2.000000e+01  2.500000e+01  4.440892e-16  4.440892e-16   \n",
       "1  40.0  1.200000e+01  2.300000e+01  2.500000e+01  4.440892e-16  4.440892e-16   \n",
       "2  10.0  4.440892e-16  4.440892e-16  4.440892e-16  1.000000e+01  2.000000e+01   \n",
       "3  10.0  3.000000e+00  4.440892e-16  4.440892e-16  1.000000e+01  2.000000e+01   \n",
       "\n",
       "             Tm  \n",
       "0  4.440892e-16  \n",
       "1  4.440892e-16  \n",
       "2  2.500000e+01  \n",
       "3  2.500000e+01  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eps_imputed = aggdf.copy()\n",
    "eps_imputed.loc[:, :] = simple_impute(aggdf.values, method='eps')\n",
    "eps_imputed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>Ti</th>\n",
       "      <th>Ge</th>\n",
       "      <th>Tm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>45.0</td>\n",
       "      <td>10.00</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.00</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40.0</td>\n",
       "      <td>12.00</td>\n",
       "      <td>23.0</td>\n",
       "      <td>25.00</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10.0</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.25</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10.0</td>\n",
       "      <td>3.00</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.25</td>\n",
       "      <td>10.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   SiO2  Al2O3   MgO    CaO    Ti    Ge     Tm\n",
       "0  45.0  10.00  20.0  25.00   0.1   0.2   0.25\n",
       "1  40.0  12.00  23.0  25.00   0.1   0.2   0.25\n",
       "2  10.0   0.03   0.2   0.25  10.0  20.0  25.00\n",
       "3  10.0   3.00   0.2   0.25  10.0  20.0  25.00"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "limit_imputed = aggdf.copy()\n",
    "limit_imputed.loc[:, :] = simple_impute(aggdf.values, method='threshold', limits=0.01*np.nanmin(aggdf.values, axis=0))\n",
    "limit_imputed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parametric Imputation using Multivariate Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parametric imputation attempts to preserve - and ideally *restore* (detection-limits effectively truncate distributions) - the distribution of multivariate compositional data. Missing values are imputed using regression against other variables. However, there remain difficulties with regards to imputation:\n",
    "* Omitted/not measured vs. 'below detection' values have different overall expectations: below detection values have upper bounds (and also upper error bounds), omitted values may be significant quantities, but are not well constrained\n",
    "* An iterative algorithm is needed (e.g. expectation-maximisation) as imputed values alter the closure operator and hence adjust other compostional components, if only slightly. Iteration continues until the imputed dataset resembles the original dataset within some specified tolerance (i.e. minimal perturbation to attain workable data)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use the EMCOMP algorithm of Palarea-Albaladejo and Martin-Fernandez (2008) to impute below-detection limit values. This requires at least one component free of zeros as a divisor. The method attempts to preserve the overall non-nan mean and covariance structure. This algorithm works well for relatively low-dimension data; future work will attempt to replace this algorithm with a bayesian imputation method along similar lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Limits at 0.1 percentile level.\n",
      "Proportion BDL: \t\t0.001\n",
      "Maximum Abs. Difference: \t0.25%\n",
      "\n",
      "Limits at 1 percentile level.\n",
      "Proportion BDL: \t\t0.006\n",
      "Maximum Abs. Difference: \t0.72%\n",
      "\n",
      "Limits at 10.0 percentile level.\n",
      "Proportion BDL: \t\t0.060\n",
      "Maximum Abs. Difference: \t9.28%\n",
      "\n",
      "Limits at 20.0 percentile level.\n",
      "Proportion BDL: \t\t0.120\n",
      "Maximum Abs. Difference: \t15.05%\n"
     ]
    }
   ],
   "source": [
    "def replace_below_detection(arr, limits=[0.1]):\n",
    "    \n",
    "    if len(limits) == 1:\n",
    "        arr[arr < lim] = np.nan\n",
    "    else:\n",
    "        for ix, lim in enumerate(limits):\n",
    "            arr[arr[:, ix] < lim, ix] = np.nan\n",
    "        \n",
    "    return arr\n",
    "\n",
    "\n",
    "for pctl in [0.1, 1, 10., 20.]:\n",
    "    print(f'\\nLimits at {pctl} percentile level.')\n",
    "    arrshape = (1000, 5)\n",
    "    arr = np.random.rand(*arrshape) * [1.0, 0.2, 0.01, 1.0, 0.5]\n",
    "    arr = np.divide(arr, np.nansum(arr, axis=1)[:, np.newaxis])\n",
    "\n",
    "    limits = np.percentile(arr, pctl, axis=0)\n",
    "\n",
    "    arr2 = arr.copy()\n",
    "    arr2[:, 2:] = replace_below_detection(arr2[:, 2:], limits=limits[2:])\n",
    "    imputed_arr, s, n = EMCOMP(arr2, limits)\n",
    "\n",
    "    print(f'Proportion BDL: \\t\\t{arr2[~np.isfinite(arr2)].size / arr2.size:2.3f}')\n",
    "    print(f'Maximum Abs. Difference: \\t{(imputed_arr-arr).max()*100:2.2f}%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}