<jupyter_start><jupyter_text>3. Missing Data and ImputationCompiled by [<NAME>](mailto:<EMAIL>) for C3DIS 2018  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).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.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).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.<jupyter_code>%matplotlib inline
%load_ext autoreload
%autoreload 2
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

sys.path.insert(0, './src')
from geochem import *
from compositions import *
from EMimputation import EMCOMP
def simple_impute(arr, method='mean', threshold=0.01, limits=[]):
    """
    """
    if method == 'mean':
        # 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
        replace_vals = np.nanmean(arr, axis=0)
    elif method == 'eps': # constant arbitrarily small value
        replace_vals = 2*np.finfo(float).eps * np.ones(arr.shape)
    elif method == 'threshold': # detection limits or some ratio of the minimum values
        if not len(limits):
            replace_vals = np.tile(np.nanmin(arr, axis=0), (arr.shape[0])).reshape(arr.shape)
        else:
            assert len(limits) == arr.shape[1]
            replace_vals = np.tile(np.array(limits), (arr.shape[0])).reshape(arr.shape)
    else:
        raise NotImplementedError
        
    return np.where(~np.isfinite(arr), replace_vals, arr)

r1 = pd.DataFrame(dict(SiO2=[45.0], Al2O3=[10.0], MgO=[20.0], CaO=[25.0]),)
r2 = pd.DataFrame(dict(SiO2=[40.0], Al2O3=[12.0], MgO=[23.0], CaO=[25.0]),)
r3 = pd.DataFrame(dict(SiO2=[10.0], Ti=[10.0], Ge=[20.0], Tm=[25.0]),)
r4 = pd.DataFrame(dict(SiO2=[10.0], Al2O3=[3.0], Ti=[10.0], Ge=[20.0], Tm=[25.0]),)
aggdf = r1
for r in [r2, r3, r4]:
    aggdf = aggdf.append(r, sort=False)
    aggdf = aggdf.reset_index(drop=True)
aggdf
mean_imputed = aggdf.copy()
mean_imputed.loc[:, :] = simple_impute(aggdf.values, method='mean')
mean_imputed
eps_imputed = aggdf.copy()
eps_imputed.loc[:, :] = simple_impute(aggdf.values, method='eps')
eps_imputed
limit_imputed = aggdf.copy()
limit_imputed.loc[:, :] = simple_impute(aggdf.values, method='threshold', limits=0.01*np.nanmin(aggdf.values, axis=0))
limit_imputed<jupyter_output><empty_output><jupyter_text>Parametric Imputation using Multivariate Regression 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:* 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* 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). 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.<jupyter_code>def replace_below_detection(arr, limits=[0.1]):
    
    if len(limits) == 1:
        arr[arr < lim] = np.nan
    else:
        for ix, lim in enumerate(limits):
            arr[arr[:, ix] < lim, ix] = np.nan
        
    return arr


for pctl in [0.1, 1, 10., 20.]:
    print(f'\nLimits at {pctl} percentile level.')
    arrshape = (1000, 5)
    arr = np.random.rand(*arrshape) * [1.0, 0.2, 0.01, 1.0, 0.5]
    arr = np.divide(arr, np.nansum(arr, axis=1)[:, np.newaxis])

    limits = np.percentile(arr, pctl, axis=0)

    arr2 = arr.copy()
    arr2[:, 2:] = replace_below_detection(arr2[:, 2:], limits=limits[2:])
    imputed_arr, s, n = EMCOMP(arr2, limits)

    print(f'Proportion BDL: \t\t{arr2[~np.isfinite(arr2)].size / arr2.size:2.3f}')
    print(f'Maximum Abs. Difference: \t{(imputed_arr-arr).max()*100:2.2f}%')<jupyter_output>Limits at 0.1 percentile level.
Proportion BDL: 		0.001
Maximum Abs. Difference: 	0.25%

Limits at 1 percentile level.
Proportion BDL: 		0.006
Maximum Abs. Difference: 	0.72%

Limits at 10.0 percentile level.
Proportion BDL: 		0.060
Maximum Abs. Difference: 	9.28%

Limits at 20.0 percentile level.
Proportion BDL: 		0.120
Maximum Abs. Difference: 	15.05%
{'fraction_non_alphanumeric': 0.09849700860936816, 'fraction_numerical': 0.026119947468262076, 'mean_word_length': 4.47444089456869, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '501120', 'license': 'MIT', 'n_tokens_mistral': 2167, 'n_tokens_neox': 1983, 'n_words': 781, 'path': '03_Missing_Data.ipynb', 'repo_name': 'morganjwilliams/exploratory-geochemistry', 'size': '6879'}