Repository: WenjieDu/PyGrinder (language: python)

Issue:
Add Little test for MCAR 
### 1. Feature description

https://www.jstor.org/stable/2290157

### 2. Motivation

To detect a missing pattern whether is MCAR.

### 3. Your contribution

Will create a PR to add this.


Discussion / hints:

Fix patch:
diff --git a/pygrinder/__init__.py b/pygrinder/__init__.py
index e87e54a..26a35db 100644
--- a/pygrinder/__init__.py
+++ b/pygrinder/__init__.py
@@ -24,7 +24,7 @@
 __version__ = "0.1.1"
 
 try:
-    from pygrinder.missing_completely_at_random import mcar
+    from pygrinder.missing_completely_at_random import mcar, mcar_little_test
     from pygrinder.missing_at_random import mar_logistic
     from pygrinder.missing_not_at_random import mnar_x, mnar_t
     from pygrinder.utils import (
@@ -37,6 +37,7 @@
 __all__ = [
     "__version__",
     "mcar",
+    "mcar_little_test",
     "mar_logistic",
     "mnar_x",
     "mnar_t",
diff --git a/pygrinder/missing_completely_at_random/__init__.py b/pygrinder/missing_completely_at_random/__init__.py
index 90a7f57..bf726d1 100644
--- a/pygrinder/missing_completely_at_random/__init__.py
+++ b/pygrinder/missing_completely_at_random/__init__.py
@@ -5,8 +5,10 @@
 # Created by Wenjie Du <wenjay.du@gmail.com>
 # License: GLP-v3
 
+from .little_test import mcar_little_test
 from .mcar import mcar
 
 __all__ = [
     "mcar",
+    "mcar_little_test",
 ]
diff --git a/setup.py b/setup.py
index b2365d2..41fa2cf 100644
--- a/setup.py
+++ b/setup.py
@@ -17,20 +17,24 @@
     url="https://github.com/WenjieDu/PyGrinder",
     download_url="https://github.com/WenjieDu/PyGrinder/archive/main.zip",
     keywords=[
-        "missing data",
-        "missing values",
         "data corruption",
         "incomplete data",
-        "partial observation",
         "data mining",
         "pypots",
         "missingness",
+        "partially observed",
+        "irregular sampled",
+        "partially-observed time series",
+        "incomplete time series",
+        "missing data",
+        "missing values",
+        "pypots",
     ],
     packages=find_packages(exclude=["tests"]),
     include_package_data=True,
     install_requires=[
         "numpy",
-        "scikit_learn",
+        "scipy",
         "pandas",
     ],
     setup_requires=["setuptools>=38.6.0"],


Test patch:
diff --git a/pygrinder/missing_completely_at_random/little_test.py b/pygrinder/missing_completely_at_random/little_test.py
new file mode 100644
index 0000000..0c55621
--- /dev/null
+++ b/pygrinder/missing_completely_at_random/little_test.py
@@ -0,0 +1,64 @@
+# Created by Wenjie Du <wenjay.du@gmail.com>
+# License: GPL-v3
+
+from typing import Union
+
+import numpy as np
+import pandas as pd
+from scipy.stats import chi2
+
+
+def mcar_little_test(X: Union[pd.DataFrame, np.ndarray]) -> float:
+    """Little's MCAR Test. Refer to :cite:`little1988TestMCAR` for more details.
+
+    Notes
+    -----
+    This implementation is inspired by
+    https://github.com/RianneSchouten/pyampute/blob/master/pyampute/exploration/mcar_statistical_tests.py
+    """
+
+    if isinstance(X, np.ndarray):
+        assert len(X.shape) == 2
+        dataset = pd.DataFrame(X)
+    elif isinstance(X, pd.DataFrame):
+        dataset = X.copy()
+    else:
+        raise RuntimeError(f"X should be np.ndarray or pd.DataFrame, but got {type(X)}")
+
+    vars = dataset.dtypes.index.values
+    n_features = dataset.shape[1]
+
+    # mean and covariance estimates
+    # ideally, this is done with a maximum likelihood estimator
+    global_mean = dataset.mean()
+    global_covariance = dataset.cov()
+
+    # set up missing data patterns
+    r = 1 * dataset.isnull()
+    mdp = np.dot(r, list(map(lambda x: pow(2, x), range(n_features))))
+    sorted_mdp = sorted(np.unique(mdp))
+    n_pat = len(sorted_mdp)
+    correct_mdp = list(map(lambda x: sorted_mdp.index(x), mdp))
+    dataset["mdp"] = pd.Series(correct_mdp, index=dataset.index)
+
+    # calculate statistic and df
+    pj = 0
+    d2 = 0
+    for i in range(n_pat):
+        dataset_temp = dataset.loc[dataset["mdp"] == i, vars]
+        select_vars = ~dataset_temp.isnull().any()
+        pj += np.sum(select_vars)
+        select_vars = vars[select_vars]
+        means = dataset_temp[select_vars].mean() - global_mean[select_vars]
+        select_cov = global_covariance.loc[select_vars, select_vars]
+        mj = len(dataset_temp)
+        parta = np.dot(
+            means.T, np.linalg.solve(select_cov, np.identity(select_cov.shape[1]))
+        )
+        d2 += mj * (np.dot(parta, means))
+
+    df = pj - n_features
+
+    # perform test and output p value
+    p_value = 1 - chi2.cdf(d2, df)
+    return p_value
diff --git a/tests/test_pygrinder.py b/tests/test_pygrinder.py
index bd69795..4ed79c1 100644
--- a/tests/test_pygrinder.py
+++ b/tests/test_pygrinder.py
@@ -11,9 +11,12 @@
 
 from pygrinder import (
     mcar,
+    mcar_little_test,
     mar_logistic,
     mnar_x,
     mnar_t,
+    masked_fill,
+    cal_missing_rate,
 )
 
 DEFAULT_MISSING_RATE = 0.1
@@ -26,22 +29,24 @@ def test_0_mcar(self):
         X_intact, X_with_missing, missing_mask, indicating_mask = mcar(
             X, p=DEFAULT_MISSING_RATE, nan=NAN
         )
-        shape_product = np.product(X_intact.shape)
-        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
-        assert (
-            round(actual_missing_rate, 1) == DEFAULT_MISSING_RATE
-        ), f"Actual missing rate is {actual_missing_rate}, not given {DEFAULT_MISSING_RATE}"
         assert np.sum(X_with_missing[(1 - missing_mask).astype(bool)]) == NAN * np.sum(
             1 - missing_mask
         )
+        X_with_missing = masked_fill(X_with_missing, 1 - missing_mask, np.nan)
+        actual_missing_rate = cal_missing_rate(X_with_missing)
+        assert (
+            round(actual_missing_rate, 1) == DEFAULT_MISSING_RATE
+        ), f"Actual missing rate is {actual_missing_rate}, not given {DEFAULT_MISSING_RATE}"
+        test_pvalue = mcar_little_test(X_with_missing.reshape(128, -1))
+        print(f"MCAR Little test p_value: {test_pvalue}")
 
     def test_1_mar(self):
         X = np.random.randn(128, 36)
         X_intact, X_with_missing, missing_mask, indicating_mask = mar_logistic(
             X, obs_rate=0.1, missing_rate=0.2
         )
-        shape_product = np.product(X_intact.shape)
-        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        X_with_missing = masked_fill(X_with_missing, 1 - missing_mask, np.nan)
+        actual_missing_rate = cal_missing_rate(X_with_missing)
         assert (
             round(actual_missing_rate, 1) > 0
         ), f"Actual missing rate is {actual_missing_rate}"
@@ -53,8 +58,8 @@ def test_2_mnar(self):
         X_intact, X_with_missing, missing_mask, indicating_mask = mnar_x(
             X, offset=0, nan=NAN
         )
-        shape_product = np.product(X_intact.shape)
-        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        X_with_missing = masked_fill(X_with_missing, 1 - missing_mask, np.nan)
+        actual_missing_rate = cal_missing_rate(X_with_missing)
         assert (
             round(actual_missing_rate, 1) > 0
         ), f"Actual missing rate is {actual_missing_rate}"
@@ -63,8 +68,8 @@ def test_2_mnar(self):
         X_intact, X_with_missing, missing_mask, indicating_mask = mnar_t(
             X, cycle=20, pos=10, scale=3, nan=NAN
         )
-        shape_product = np.product(X_intact.shape)
-        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        X_with_missing = masked_fill(X_with_missing, 1 - missing_mask, np.nan)
+        actual_missing_rate = cal_missing_rate(X_with_missing)
         assert (
             round(actual_missing_rate, 1) > 0
         ), f"Actual missing rate is {actual_missing_rate}"