Two-sample Kolmogorov-Smirnov test to compare pre and post imputation

[Zethson]

Description of feature

We could implement the two-sample Kolmogorov-Smirnov test which we could use to test pre and post imputation distributions

We tested the null hypothesis considering that the two sets of observations were drawn from the same unknown probability distribution, using a two-sample Kolmogorov-Smirnov test, as implemented by the ks.test function in the stats R package. The two sets of observations considered were either the variable before and after introducing the missing values, or the variable before introducing the missing values compared with the variable after imputation, or the variable after introducing missing values and the variable after imputation. Note that in a typical imputation situation, true values are not available, and thus only comparison between the set of non-missing values and the set of values after imputation would be possible.

[eroell]

Btw do you have the reference for this quote above?

[Zethson]

https://link.springer.com/article/10.1186/s12874-023-02125-x

sorry I should have included it.

[eroell]

So a rough outline here would be

• add this using most likely the scipy implementation: from the handling, it probably should feel similar to how e.g. Little's Test is implemented in ehrapy
• Write a few tests checking that the computation for a tiny synthetic sample works, and check that the arguments are handled properly
• Write a helpful docstring
• Consider adding to bias notebook

[Zethson]

And maybe add this to the bias notebook?

[nicolassidoux]

Hi guys,
First contributions here!

My plan was to create a new function in preprocessing/_quality_control as:

def ks_test(adata: AnnData,
            adata_imputed: AnnData,
            layer: str | None = None,
            method: Literal["auto", "exact", "asymp"] = "auto",
            p_threshold: float = 0.95
            ) -> list[str]

with:

Args:
adata: Encoded annotated data matrix, before imputation.
adata_imputed: Encoded annotated data matrix, after imputation.
layer: Layer to apply the test to. Uses X matrix if set to None.
method: Defines the method used for calculating the p-value.
The following options are available (default is 'auto'):
* 'auto' : use 'exact' for small size arrays, 'asymp' for large
* 'exact' : use exact distribution of test statistic
* 'asymp' : use asymptotic distribution of test statistic
p_threshold: P-value threshold to consider a feature as significantly impacted

Returns:
The list of features which calculated p-value is less than specified p_threshold, potientially indicating a
non-negligible impact of the imputation on the data.

In regards to test I wanted to create two functions:

• Test "same datasets":
  • Create an AnnData dataset A with a random amount of vars and obs
  • Fill it with random values
  • Copy the dataset A to dataset B
  • For each cell of the dataset B, replace the contained value with NaN, based on a probability
  • Perform a KS test for A and B
  • Return from ks_test() should be empty
• Test "different datasets":
  • Create an AnnData dataset A with a random amount of vars and obs
  • Fill it with random values
  • DON'T copy dataset A to dataset B, but instead create a new one of same shape but filled with different random values
  • For each cell of the dataset B, replace the contained value with NaN, based on a probability
  • Perform a KS test for A and B
  • Return from ks_test() should NOT be empty

So, do you guys think i'm having the right approach here? If yes, I was thinking to generate the NaN probability as a value in the range of [0.01, 0.2]. Do that feel realistic? Also, the test results depend on the chosen p-value threshold. After some tests, it seems that 0.95 gives good results. Again, what do you think?

[nicolassidoux]

Also, out of curiosity I wanted to check the effect of imputation on KS test. It looks like the calculated p-values are correlated with the size of the dataset, the type of imputation and the proportion of missing data.
Code used to generate the data:

import numpy as np
import scipy as sp
import pandas as pd
import anndata as ad
import ehrapy as ep
import statistics
from lamin_utils import logger


rng = np.random.default_rng()
n_obs = 58000
n_vars = 80
total_size = n_obs * n_vars

logger.disabled = True

full_adata = ad.AnnData(pd.DataFrame(
    np.random.uniform(0.0, 100.0, size=(n_obs, n_vars)),
    columns=[f'Var{i + 1}' for i in range(n_vars)]))

proba_missing = 0.001
while True:
    hollowed_adata = full_adata.copy()

    for i in range(hollowed_adata.shape[0]):
        for j in range(hollowed_adata.shape[1]):
            if rng.uniform() < proba_missing:
                hollowed_adata.X[i, j] = np.nan

    missing_count = ep.pp.qc_metrics(hollowed_adata)[0]["missing_values_abs"].sum()
    #ep.pp.simple_impute(hollowed_adata)
    ep.pp.knn_impute(hollowed_adata)
    meanpval = statistics.mean(sp.stats.ks_2samp(full_adata.to_df(), hollowed_adata.to_df(), nan_policy="omit").pvalue)
    print(missing_count / total_size * 100.0, ",", meanpval)
    if meanpval < 0.05 or proba_missing >= 1.0:
        exit()
    else:
        proba_missing += 0.001

FYI, this is what I got:
Warning: X scales not the same


Maybe it will be of use?

[eroell]

def ks_test(adata: AnnData, # call adata_1 instead?
            adata_imputed: AnnData, # call adata_2 instead?
            layer: str | None = None, # indicate layer for both, so adata_1_layer, adata_2_layer
            method: Literal["auto", "exact", "asymp"] = "auto",
            p_threshold: float = 0.95
            ) -> list[str]

with:

Args:
adata: Encoded annotated data matrix, before imputation. # see naming above
adata_imputed: Encoded annotated data matrix, after imputation. # see naming above
layer: Layer to apply the test to. Uses X matrix if set to None. # see comment on specifying for both adatas above
method: Defines the method used for calculating the p-value.
The following options are available (default is 'auto'):
* 'auto' : use 'exact' for small size arrays, 'asymp' for large
* 'exact' : use exact distribution of test statistic
* 'asymp' : use asymptotic distribution of test statistic
p_threshold: P-value threshold to consider a feature as significantly impacted

Returns:
The list of features which calculated p-value is less than specified p_threshold, potientially indicating a
non-negligible impact of the imputation on the data.

Write test statistic and p-value into .var["ks_test_statistic"] and .var["ks_p_value"], and a boolean number about the p-value being reached into .var["ks_significant"] instead of returning a list of features. The list of features could then be quite easily be grabbed by the user with this being done

[Zethson]

Some work is on this branch: https://github.com/theislab/ehrapy/commits/enhancement/issue-652/