Repository: jakobrunge/tigramite (language: python)

Issue:
Developer
Contains the code, the corresponding tutorial and unit tests for the UAI2023 paper on pairwise independence testing


Discussion / hints:

Fix patch:
diff --git a/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb b/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb
new file mode 100644
index 00000000..8466d094
--- /dev/null
+++ b/tutorials/causal_discovery/tigramite_tutorial_pairwise_mult_ci.ipynb
@@ -0,0 +1,279 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multivariate conditional independence testingwith `TIGRAMITE`\n",
+    "\n",
+    "TIGRAMITE is a time series analysis python module. It allows to reconstruct graphical models (conditional independence graphs) from discrete or continuously-valued time series based on the PCMCI framework and create high-quality plots of the results.\n",
+    "\n",
+    "This tutorial explains how the multivariate conditional independence test PairwiseMultCI works. For the theoretical details, we refer to\n",
+    " - Tom Hochsprung, Jonas Wahl, Andreas Gerhardus, Urmi Ninad, and Jakob Runge.\n",
+    "      Increasing Effect Sizes of Pairwise Conditional Independence Tests between Random Vectors. UAI2023, 2023.\n",
+    "\n",
+    "We imagine that such a multivariate test becomes relevant for vector-valued causal discovery, that is, causal discovery where certain variables are grouped together and one is only interested in learning causal relationsships between groups of variables.\n",
+    "That research area is still active, for recent work, see\n",
+    "- Wahl, J., Ninad, U., & Runge, J. (2023). Foundations of Causal Discovery on Groups of Variables. arXiv preprint arXiv:2306.07047.\n",
+    "- Wahl, J., Ninad, U., & Runge, J. (2022). Vector causal inference between two groups of variables. arXiv preprint arXiv:2209.14283."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Imports\n",
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "import tigramite\n",
+    "\n",
+    "from tigramite.independence_tests.parcorr import ParCorr\n",
+    "from tigramite.independence_tests.pairwise_CI import PairwiseMultCI"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Multivariate CI-Testing\n",
+    "\n",
+    "### Overview\n",
+    "\n",
+    "Suppose one has random vectors $X, Y, Z$ (each potentially multivariate) and suppose that one  wants to decide whether $X$ is independent of $Y$ given $Z$ (that is, the conditional densities (assuming they exist) factorize: $\\forall x,y,z:\\;p_{X,Y|Z=z}(x,y) = p_{X|Z=z}(x)\\cdot p_{Y|Z=z}(y)$.)\n",
+    "\n",
+    "One approach to test independence between $X$ and $Y$ given $Z$ is to test whether all components $X_i$ and $Y_j$ are independent given $Z$.\n",
+    "The conditional independence test PairwiseMultCI builds on this pairwise idea. PairwiseMultCI tests independence between each pair of components $X_i$ and $Y_j$ given $(Z,S_{ij}).$ Here, $S_{ij}$ is a set that consists of components of either $X$ or $Y$ that are independent of $Y_j$ or $X_i$. One can show that including these components increases the effect size of the corresponding test.\n",
+    "\n",
+    "In practice, these sets $S_{ij}$ are usually not known a priori and hence need to be learned.\n",
+    "PairwiseMultCI first learns these conditional independencies on one part of the sample and then uses these conditional independencies to do the testing on the second part of the sample. The sample splitting ratio (``pre_step_sample_fraction``) is a hyperparameter, as is the significance level ($\\alpha_{pre}$) for learning independencies in the first step. (We think that only a small part of the sample should be used for learning independencies in the first step, moreover, $\\alpha_{pre}$ should be rather large to only mark something as independent if there is strong evidence).\n",
+    "\n",
+    "We remark that PairwiseMultCI is a \"meta-algorithm\", that is, it is a framework that works with several univariate test statistics. Because of that, PairwiseMultCI also has the parameter cond_ind_test, which can be any other conditional independence test from tigramite/independence_tests"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Simulations - Part 1\n",
+    "\n",
+    "We start by generating a dataset where $X$ and $Z$ are univariate, and $Y=(Y_1,Y_2)$ is bivariate."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T = 100\n",
+    "z = np.random.normal(0, 1, T). reshape(T, 1)\n",
+    "x = np.random.normal(0, 1, T). reshape(T, 1) + 0.4 * z\n",
+    "y1 = np.random.normal(0, 1, T). reshape(T, 1) + 0.2 * z\n",
+    "y2 = y1 + 0.3 * np.random.normal(0, 1, T). reshape(T, 1) + 0.5 * x\n",
+    "y = np.hstack((y1, y2))\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We then print out the maximum of the univariate test statistics and the aggregated p-value for PairwiseMultCI as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "test_stat 0.8330221744959563\n",
+      "pval 6.30177446825862e-21\n"
+     ]
+    }
+   ],
+   "source": [
+    "# alpha_pre: significance level for the first step\n",
+    "# pre_step_sample_fraction: relative size of the sample for the first step\n",
+    "# cond_ind_test: respective univariate conditional independence test\n",
+    "ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = ParCorr())\n",
+    "test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
+    "print(\"test_stat\", test_stat)\n",
+    "print(\"pval\", pval)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Simulations - Part 2\n",
+    "\n",
+    "We now illustrate that PairwiseMultCI performs especially well relative to other approaches when the within-$Y$ (or within-$X$) dependence is strong. For that we consider similar models as above and calculate the rejection rate over $100$ replications. We also include one competitor  approach that just tests pairwise independencies (as explained in the second paragraph of the Overview-Section)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Large within-$Y$ dependence:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Rejection Rate PairwiseMultCI: 0.53\n",
+      "Rejection Rate Competitor: 0.11\n"
+     ]
+    }
+   ],
+   "source": [
+    "np.random.seed(2023)\n",
+    "T = 100\n",
+    "# number replications:\n",
+    "repls = 100\n",
+    "alpha = 0.05\n",
+    "rejection_bool = np.zeros(repls)\n",
+    "rejection_bool_competitor = np.zeros(repls)\n",
+    "par_corr = ParCorr()\n",
+    "for i in np.arange(0, repls):\n",
+    "    # data generation:\n",
+    "    z = np.random.normal(0, 1, T). reshape(T, 1)\n",
+    "    x = np.random.normal(0, 1, T). reshape(T, 1) + 0.4 * z\n",
+    "    y1 = np.random.normal(0, 1, T). reshape(T, 1) + 0.2 * z\n",
+    "    y2 = y1 + 0.3 * np.random.normal(0, 1, T). reshape(T, 1) + 0.1 * x\n",
+    "    y = np.hstack((y1, y2))\n",
+    "    # PairwiseMultCI\n",
+    "    ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
+    "    test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval <= alpha):\n",
+    "        rejection_bool[i] = 1  \n",
+    "           \n",
+    "    ## competitor:\n",
+    "    pval1 = par_corr.run_test_raw(x, y1, z)[1]\n",
+    "    pval2 = par_corr.run_test_raw(x, y2, z)[1]  \n",
+    "    if (np.min((pval1, pval2)) * 2 <= alpha):\n",
+    "        rejection_bool_competitor[i] = 1 \n",
+    "print(\"Rejection Rate PairwiseMultCI:\", np.mean(rejection_bool))\n",
+    "print(\"Rejection Rate Competitor:\",np.mean(rejection_bool_competitor))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We see that PairwiseMultCI strongly outperforms the competitor."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### No within-$Y$ dependence:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Rejection Rate PairwiseMultCI: 0.77\n",
+      "Rejection Rate Competitor: 0.85\n"
+     ]
+    }
+   ],
+   "source": [
+    "np.random.seed(2023)\n",
+    "T = 100\n",
+    "# number replications:\n",
+    "repls = 100\n",
+    "alpha = 0.05\n",
+    "rejection_bool = np.zeros(repls)\n",
+    "rejection_bool_competitor = np.zeros(repls)\n",
+    "par_corr = ParCorr()\n",
+    "for i in np.arange(0, repls):\n",
+    "    # data generation:\n",
+    "    z = np.random.normal(0, 1, T). reshape(T, 1)\n",
+    "    x = np.random.normal(0, 1, T). reshape(T, 1) + 0.4 * z\n",
+    "    y1 = np.random.normal(0, 1, T). reshape(T, 1) + 0.2 * z\n",
+    "    y2 = 0.3 * np.random.normal(0, 1, T). reshape(T, 1) + 0.1 * x\n",
+    "    y = np.hstack((y1, y2))\n",
+    "    # PairwiseMultCI\n",
+    "    ci_test = PairwiseMultCI(alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = par_corr)\n",
+    "    test_stat, pval = ci_test.run_test_raw(x = x, y = y, z = z)\n",
+    "    if (pval <= alpha):\n",
+    "        rejection_bool[i] = 1  \n",
+    "           \n",
+    "    ## competitor:\n",
+    "    pval1 = par_corr.run_test_raw(x, y1, z)[1]\n",
+    "    pval2 = par_corr.run_test_raw(x, y2, z)[1]  \n",
+    "    if (np.min((pval1, pval2)) * 2 <= alpha):\n",
+    "        rejection_bool_competitor[i] = 1 \n",
+    "print(\"Rejection Rate PairwiseMultCI:\", np.mean(rejection_bool))\n",
+    "print(\"Rejection Rate Competitor:\",np.mean(rejection_bool_competitor))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We see that in case of no dependence between $Y_1$ and $Y_2$, PairwiseMultCI performs slightly worse than the competitor."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Remark:\n",
+    "\n",
+    "It is also possible to instantiate PairwiseMultCI with the option significance = \"fixed_thres\". Strictly speaking, setting the significance of the cond_ind_test_object is something else, however, at the moment, we will work with fixed_thres if at least at one occasion significance = fixed_thres. When working with fixed_thres, instead of working with the significance level $\\alpha_{pre},$ one needs to specify a threshold for the pre-step, called fixed_thres_pres. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Conclusion:\n",
+    "\n",
+    "- PairwiseMultCI is only relevant when $X$ and/or $Y$ are multivariate.\n",
+    "- PairwiseMultCI relies on the underlying assumption that independence is only violated for a few components $X_i$ and $Y_j$.\n",
+    "- PairwiseMultCI works especially well (relative to other approaches) when the within-$X$ or within-$Y$ dependence is large."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "e1",
+   "language": "python",
+   "name": "e1"
+  },
+  "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.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}


Test patch:
diff --git a/tests/test_PairwiseMultCI.py b/tests/test_PairwiseMultCI.py
new file mode 100644
index 00000000..501fb0d9
--- /dev/null
+++ b/tests/test_PairwiseMultCI.py
@@ -0,0 +1,98 @@
+"""
+Tests for independence_tests.py.
+"""
+from __future__ import print_function
+import numpy as np
+import pytest
+from scipy import stats
+
+from tigramite.independence_tests.parcorr import ParCorr
+from tigramite.independence_tests.robust_parcorr import RobustParCorr
+from tigramite.independence_tests.pairwise_CI import PairwiseMultCI
+from tigramite.independence_tests.cmiknn import CMIknn
+from tigramite.independence_tests.regressionCI import RegressionCI
+# PairwiseMultCI TESTING ##################################################
+@pytest.fixture(params=[
+    # Generate PairwiseMultCI test instances
+    # test PairwiseMultCI for different combinations of:
+    # significance (of the univariate tests),  alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre
+    ('fixed_thres', 0.5, 0.2, ParCorr, 1, 0.5),
+    ('fixed_thres', 0.5, 0.5, ParCorr, 1, 0.5),
+    ('fixed_thres', 0.8, 0.2, ParCorr, 1, 0.5),
+    ('fixed_thres', 0.8, 0.5, ParCorr, 1, 0.5),
+    ('analytic', 0.5, 0.2, ParCorr, None, None),
+    ('shuffle_test', 0.5, 0.2, ParCorr, None, None),
+    ('analytic', 0.5, 0.5, ParCorr, None, None),
+    ('shuffle_test', 0.5, 0.5, ParCorr, None, None),
+    ('analytic', 0.8, 0.2, ParCorr, None, None),
+    ('shuffle_test', 0.8, 0.2, ParCorr, None, None),
+    ('analytic', 0.8, 0.5, ParCorr, None, None),
+    ('shuffle_test', 0.8, 0.5, ParCorr, None, None),
+    ('analytic', 0.5, 0.2, RobustParCorr, None, None),
+    ('shuffle_test', 0.5, 0.2, RobustParCorr, None, None),
+    ('analytic', 0.5, 0.5, RobustParCorr, None, None),
+    ('shuffle_test', 0.5, 0.5, RobustParCorr, None, None),
+    ('analytic', 0.8, 0.2, RobustParCorr, None, None),
+    ('shuffle_test', 0.8, 0.2, RobustParCorr, None, None),
+    ('analytic', 0.8, 0.5, RobustParCorr, None, None),
+    ('shuffle_test', 0.8, 0.5, RobustParCorr, None, None),
+    ('shuffle_test', 0.5, 0.2, CMIknn, None, None),
+    ('shuffle_test', 0.5, 0.5, CMIknn, None, None),
+    ('shuffle_test', 0.8, 0.2, CMIknn, None, None),
+    ('shuffle_test', 0.8, 0.5, CMIknn, None, None),
+])
+
+def pairwise_mult_ci(request):
+    # Unpack the parameters
+    sig, alpha_pre, pre_step_sample_fraction, cond_ind_test, cond_ind_test_thres, cond_ind_test_thres_pre = request.param
+    # Generate the par_corr_wls independence test
+    if sig != "fixed_thres":
+        return PairwiseMultCI(cond_ind_test = cond_ind_test(significance = sig),
+                          alpha_pre = alpha_pre,
+                          pre_step_sample_fraction = pre_step_sample_fraction)
+    else:
+        return PairwiseMultCI(cond_ind_test = cond_ind_test(significance = sig),
+                          alpha_pre = None,
+                          pre_step_sample_fraction=pre_step_sample_fraction,
+                          significance= sig,
+                          fixed_thres_pre = 2
+                          )
+
+@pytest.fixture(params=[
+    # Generate PairwiseMultCI test instances
+    # test PairwiseMultCI for different combinations of:
+    # seed,  true_dep, T (=sample size)
+    (123, 0, 100),
+    (123, 0, 1000),
+    (123, 0.2, 100),
+    (123, 0.2, 1000),
+    (123, 0.5, 100),
+    (123, 0, 1000),
+    (46, 0, 100),
+    (46, 0, 1000),
+    (46, 0.2, 100),
+    (46, 0.2, 1000),
+    (46, 0.5, 100),
+    (46, 0, 1000),
+])
+
+def data_sample_c_cc_c(request):
+    # Set the random seed
+    seed, true_dep, T = request.param
+    np.random.seed(seed)
+    x = np.random.normal(0, 1, T).reshape(T, 1)
+    y1 = np.random.normal(0, 1, T).reshape(T, 1)
+    y2 = true_dep * x + y1 + 0.3 * np.random.normal(0, 1, T).reshape(T, 1)
+    y = np.hstack((y1,y2))
+    z = np.random.normal(0, 1, T).reshape(T, 1)
+    # Return data xyz
+    return x, y, z
+
+def test_pairwise_mult_ci(pairwise_mult_ci, data_sample_c_cc_c):
+    # Get the data sample values
+    x, y, z = data_sample_c_cc_c
+    # Get the analytic significance
+    test_result = pairwise_mult_ci.run_test_raw(x = x, y = y, z = z, alpha_or_thres=1)
+    val = test_result[0]
+    pval = test_result[1]
+    np.testing.assert_allclose(pval, 0.5, atol=0.5)
diff --git a/tigramite/independence_tests/gpdc.py b/tigramite/independence_tests/gpdc.py
index b91e1c20..0df386bc 100644
--- a/tigramite/independence_tests/gpdc.py
+++ b/tigramite/independence_tests/gpdc.py
@@ -27,6 +27,9 @@ class GaussProcReg():
     a \*.npz file containing the null distribution for different sample sizes.
     This file can then be supplied as null_dist_filename.
 
+    Assumes one-dimensional X, Y. But can be combined with PairwiseMultCI to
+    obtain a test for multivariate X, Y.
+
     Parameters
     ----------
     null_samples : int
diff --git a/tigramite/independence_tests/gpdc_torch.py b/tigramite/independence_tests/gpdc_torch.py
index ada3944f..be6d272b 100644
--- a/tigramite/independence_tests/gpdc_torch.py
+++ b/tigramite/independence_tests/gpdc_torch.py
@@ -25,6 +25,9 @@ class GaussProcRegTorch():
     a \*.npz file containing the null distribution for different sample sizes.
     This file can then be supplied as null_dist_filename.
 
+    Assumes one-dimensional X, Y. But can be combined with PairwiseMultCI to
+    obtain a test for multivariate X, Y.
+
     Parameters
     ----------
     null_samples : int
diff --git a/tigramite/independence_tests/gsquared.py b/tigramite/independence_tests/gsquared.py
index 48163558..b31c01ae 100644
--- a/tigramite/independence_tests/gsquared.py
+++ b/tigramite/independence_tests/gsquared.py
@@ -23,7 +23,8 @@ class Gsquared(CondIndTest):
     above 1000-2000 samples (depends on data). For smaller sample sizes use the
     CMIsymb class which includes a local permutation test.
 
-    Assumes one-dimensional X, Y.
+    Assumes one-dimensional X, Y. But can be combined with PairwiseMultCI to
+    obtain a test for multivariate X, Y.
 
     This method requires the scipy.stats package.
 
diff --git a/tigramite/independence_tests/pairwise_CI.py b/tigramite/independence_tests/pairwise_CI.py
new file mode 100644
index 00000000..222d6b9c
--- /dev/null
+++ b/tigramite/independence_tests/pairwise_CI.py
@@ -0,0 +1,337 @@
+"""Tigramite causal discovery for time series."""
+
+# Author: Jakob Runge <jakob@jakob-runge.com>
+#
+# License: GNU General Public License v3.0
+
+from __future__ import print_function
+from scipy import stats
+import numpy as np
+import sys
+import warnings
+
+from tigramite.independence_tests.independence_tests_base import CondIndTest
+from tigramite.independence_tests.parcorr import ParCorr
+
+class PairwiseMultCI(CondIndTest):
+    r""" Multivariate CI-test that aggregates univariate tests
+
+    Basically, the algorithm consists of two main steps:
+
+    - Step 1: On the first part of the sample, conditional independencies of the form :math:`X_i \perp Y_j|Z`
+      are learned
+
+    - Step 2: On the second part of the sample, conditional independencies
+      :math:`X_i \perp Y_j|(Z, S_{ij})` are tested, where the set :math:`S_{ij}` consists of components
+      of :math:`X` that are independent of :math:`Y_j` given :math:`Z` (or components of :math:`Y` that
+      are independent of :math:`X_i` given :math:`Z`). Using the Bonferroni method, the univariate tests
+      are then aggregated.
+
+    The main reasoning behind this two-step procedure is that the conditional independence tests in the second step
+    have larger effect sizes. One can show, that in cases where the within-:math:`X` or within-:math:`Y` dependence
+    is large, this method outperforms naive pairwise independence testing (i.e., without having the sets :math:`S_{ij}`)
+
+    For the details, we refer to:
+
+    - Tom Hochsprung, Jonas Wahl*, Andreas Gerhardus*, Urmi Ninad*, and Jakob Runge.
+      Increasing Effect Sizes of Pairwise Conditional Independence Tests between Random Vectors. UAI2023, 2023.
+
+    Parameters
+    ----------
+    alpha_pre: float
+        Significance level for the first step of the algorithm.
+        If cond_ind_test is instantiated with significance = "fixed_thres",
+        then the value of alpha_pre is not used, and can, for example, be set to None
+
+    pre_step_sample_fraction: float
+        Relative size for the first part of the sample
+        (abbreviation for "size block one")
+
+    cond_ind_test: conditional independence test object
+        This can be ParCorr or other classes from
+        ``tigramite.independence_tests`` or an external test passed as a
+        callable. This test can be based on the class
+        tigramite.independence_tests.CondIndTest.
+
+    fixed_thres_pre: In case cond_ind_test or self is instantiated with significance
+        = "fixed_thres" (at least one such setting suffices), the pre_step works with a threshold instead of significance
+        level
+
+
+    **kwargs :
+        Arguments passed on to Parent class CondIndTest.
+    """
+
+    # documentation
+    @property
+    def measure(self):
+        """
+        Concrete property to return the measure of the independence test
+        """
+        return self._measure
+
+    def __init__(self, cond_ind_test = ParCorr(), alpha_pre = 0.5, pre_step_sample_fraction = 0.2, fixed_thres_pre = None, **kwargs):
+        self._measure = 'pairwise_CI'
+        self.cond_ind_test = cond_ind_test
+        self.alpha_pre = alpha_pre
+        self.pre_step_sample_fraction = pre_step_sample_fraction
+        self.fixed_thres_pre = fixed_thres_pre
+        self.two_sided = False
+        CondIndTest.__init__(self, **kwargs)
+
+
+    def calculate_dep_measure_and_significance(self, array, xyz, data_type = None, ci_test_thres = None):
+        """Return aggregated p-value of all univariate independence tests.
+
+        First, learns conditional independencies, uses these conditional independencies to have new conditional
+        independence tests with larger effect sizes.
+
+        Parameters
+        ----------
+        array : array-like
+            data array with X, Y, Z in rows and observations in columns
+
+        xyz : array of ints
+            XYZ identifier array of shape (dim,).
+
+        data_type: array of 0's and 1's, same dimensions as array
+            relevant, for example for the CI-test RegressionCI,
+            to indicate whether variables are discrete (=1) or
+            continuous (=0).
+
+
+        Returns
+        -------
+        pval : float
+            an aggregated p-value of all univariate tests.
+        """
+
+        fixed_thres_bool = False
+        if self.significance == "fixed_thres":
+            fixed_thres_bool = True
+            self.cond_ind_test.significance = "fixed_thres"
+        if self.cond_ind_test.significance == "fixed_thres":
+            fixed_thres_bool = True
+            self.significance = "fixed_thres"
+
+        if (fixed_thres_bool) and (self.fixed_thres_pre == None):
+            raise ValueError("If significance == 'fixed_thres', fixed_thres_pre for the"
+                             " pre-step needs to be defined in initializing PairwiseMultCI.")
+
+        x_indices = np.where(xyz == 0)[0]
+        y_indices = np.where(xyz == 1)[0]
+        z_indices = np.where(xyz == 2)[0]
+
+        dim_x = x_indices.shape[0]
+        dim_y = y_indices.shape[0]
+        # what if unconditional test ...
+        dim_z = z_indices.shape[0]
+        T = array.shape[1]
+        size_first_block = np.floor(self.pre_step_sample_fraction * T)
+        size_first_block = size_first_block.astype(int)
+
+        # split the sample
+        array_s1 = array[:, 0:size_first_block]
+        x_s1 = array_s1[x_indices]
+        y_s1 = array_s1[y_indices]
+        z_s1 = array_s1[z_indices]
+
+        array_s2 = array[:, size_first_block : T]
+        x_s2 = array_s2[x_indices]
+        y_s2 = array_s2[y_indices]
+        z_s2 = array_s2[z_indices]
+
+        if data_type is not None:
+            x_type = data_type[x_indices, :].T
+            y_type = data_type[y_indices, :].T
+            z_type = data_type[z_indices, :].T
+        else:
+            x_type = None
+            y_type = None
+            z_type = None
+        if x_type is not None:
+            x_type_s1 = x_type[0:size_first_block, :]
+            x_type_s2 = x_type[size_first_block : T, :]
+        else:
+            x_type_s1 = None
+            x_type_s2 = None
+        if y_type is not None:
+            y_type_s1 = y_type[0:size_first_block, :]
+            y_type_s2 = y_type[size_first_block : T, :]
+        else:
+            y_type_s1 = None
+            y_type_s2 = None
+        if z_type is not None:
+            z_type_s1 = z_type[0:size_first_block, :]
+            z_type_s2 = z_type[size_first_block: T, :]
+        else:
+            z_type_s1 = None
+            z_type_s2 = None
+
+        ## Step 1: estimate conditional independencies
+        if fixed_thres_bool == False:
+            p_vals_pre = np.zeros((dim_x, dim_y))
+        else:
+            vals_pre = np.zeros((dim_x, dim_y))
+        for j in np.arange(0, dim_x):
+            for jj in np.arange(0, dim_y):
+                if fixed_thres_bool == False:
+                    p_vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
+                                                                   y_s1[jj].reshape(size_first_block, 1),
+                                                                   z_s1.reshape(size_first_block, dim_z),
+                                                                   x_type = x_type_s1,
+                                                                   y_type = y_type_s1,
+                                                                   z_type = z_type_s1)[1]
+                else:
+                    vals_pre[j, jj] = self.cond_ind_test.run_test_raw(x_s1[j].reshape(size_first_block, 1),
+                                                                   y_s1[jj].reshape(size_first_block, 1),
+                                                                   z_s1.reshape(size_first_block, dim_z),
+                                                                   x_type = x_type_s1,
+                                                                   y_type = y_type_s1,
+                                                                   z_type = z_type_s1,
+                                                                   alpha_or_thres=999.)[0] # just a dummy
+        if fixed_thres_bool == False:
+            indep_set = np.where(p_vals_pre > self.alpha_pre)
+        else:
+            indep_set = np.where(np.abs(vals_pre) >= self.fixed_thres_pre)
+
+        # Step 2: test conditional independencies with increased effect sizes
+        if self.cond_ind_test.significance != "fixed_thres":
+            p_vals_main = np.zeros((dim_x, dim_y))
+        else:
+            dependent_main = np.zeros((dim_x, dim_y))
+        test_stats_main = np.zeros((dim_x, dim_y))
+
+        for j in np.arange(0, dim_x):
+            for jj in np.arange(0, dim_y):
+                indicesY = np.zeros(0)
+                indicesX = np.zeros(0)
+
+                if (j in indep_set[0]):
+                    indicesX_locs = np.where(indep_set[0] == j)
+                    indicesY = np.setdiff1d(indep_set[1][indicesX_locs], jj)
+                if (jj in indep_set[1]):
+                    indicesY_locs = np.where(indep_set[1] == jj)
+                    indicesX = np.setdiff1d(indep_set[0][indicesY_locs], j)
+                lix = indicesX.shape[0]
+                liy = indicesY.shape[0]
+                if lix + liy > 0:
+                    if (lix > liy):
+                        if x_type_s2 is not None:
+                            z_type = np.hstack((z_type_s2, x_type_s2[:, indicesX]))
+                        if fixed_thres_bool == False:
+                            test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                        y_s2[jj].reshape(T - size_first_block, 1),
+                                                        np.hstack((z_s2.reshape(T - size_first_block, dim_z),
+                                                                    x_s2[indicesX].reshape(T - size_first_block, lix))),
+                                                        x_type = x_type_s2,
+                                                        y_type = y_type_s2,
+                                                        z_type = z_type_s2)
+                        else:
+                            test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                        y_s2[jj].reshape(T - size_first_block, 1),
+                                                        np.hstack((z_s2.reshape(T - size_first_block, dim_z),
+                                                                    x_s2[indicesX].reshape(T - size_first_block, lix))),
+                                                        x_type = x_type_s2,
+                                                        y_type = y_type_s2,
+                                                        z_type = z_type_s2,
+                                                        alpha_or_thres=999.) # just a dummy
+                    elif (lix <= liy):
+                        if y_type_s2 is not None:
+                            z_type_s2 = np.hstack((z_type_s2, y_type_s2[:, indicesY]))
+                        if fixed_thres_bool == False:
+                            test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                        y_s2[jj].reshape(T - size_first_block, 1),
+                                                        np.hstack((z_s2.reshape(T - size_first_block, dim_z),
+                                                                    y_s2[indicesY].reshape(T - size_first_block, liy))),
+                                                        x_type = x_type_s2,
+                                                        y_type = y_type_s2,
+                                                        z_type = z_type_s2)
+                        else:
+                            test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                        y_s2[jj].reshape(T - size_first_block, 1),
+                                                        np.hstack((z_s2.reshape(T - size_first_block, dim_z),
+                                                                    y_s2[indicesY].reshape(T - size_first_block, liy))),
+                                                        x_type = x_type_s2,
+                                                        y_type = y_type_s2,
+                                                        z_type = z_type_s2,
+                                                        alpha_or_thres=999.) # just a dummy
+
+                else:
+                    if fixed_thres_bool == False:
+                        test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                                      y_s2[jj].reshape(T - size_first_block, 1),
+                                                                      z_s2.reshape(T - size_first_block, dim_z),
+                                                                      x_type = x_type_s2,
+                                                                      y_type = y_type_s2,
+                                                                      z_type = z_type_s2)
+                    else:
+                        test_result = self.cond_ind_test.run_test_raw(x_s2[j].reshape(T - size_first_block, 1),
+                                                                      y_s2[jj].reshape(T - size_first_block, 1),
+                                                                      z_s2.reshape(T - size_first_block, dim_z),
+                                                                      x_type = x_type_s2,
+                                                                      y_type = y_type_s2,
+                                                                      z_type = z_type_s2,
+                                                                      alpha_or_thres=999.) # just a dummy
+                
+                test_stats_main[j, jj] = test_result[0]
+                if fixed_thres_bool == False:
+                    p_vals_main[j, jj] = test_result[1]
+
+
+        # Aggregate p-values
+        test_stats_aggregated = np.max(np.abs(test_stats_main))
+        if self.cond_ind_test.significance != "fixed_thres":
+            p_aggregated = np.min(np.array([np.min(p_vals_main) * dim_x * dim_y, 1]))
+        else:
+            p_aggregated = None
+
+        return test_stats_aggregated, p_aggregated
+
+
+    def get_dependence_measure(self, array, xyz, data_type=None, ci_test_thres = None):
+
+        self.dep_measure, self.signif = self.calculate_dep_measure_and_significance(array = array, xyz = xyz, data_type = data_type, ci_test_thres = ci_test_thres)
+
+        return self.dep_measure
+
+    def get_analytic_significance(self, value = None, T = None, dim = None, xyz = None, data_type=None):
+
+        return self.signif
+
+if __name__ == '__main__':
+
+    import tigramite
+    from tigramite.data_processing import DataFrame
+    import tigramite.data_processing as pp
+    import numpy as np
+    from tigramite.independence_tests.robust_parcorr import RobustParCorr
+    from tigramite.independence_tests.cmiknn import CMIknn
+    alpha = 0.05
+    ci = PairwiseMultCI(fixed_thres_pre = 0.5, cond_ind_test = ParCorr(significance = "fixed_thres"))
+    # ci = PairwiseMultCI(cond_ind_test=ParCorr(significance="fixed_thres"))
+    # ci = PairwiseMultCI(cond_ind_test=ParCorr(significance="analytic"))
+    T = 100
+    reals = 1
+    rate = np.zeros(reals)
+    np.random.seed(123)
+    for t in range(reals):
+        # continuous example
+        x = np.random.normal(0, 1, T).reshape(T, 1)
+        y1 = np.random.normal(0, 1, T).reshape(T, 1)
+        y2 = 5 * x + y1 + 0.3 * np.random.normal(0, 1, T).reshape(T, 1)
+        z = np.random.normal(0, 1, T).reshape(T, 1)
+        # discrete example
+        #x = np.random.binomial(n=10, p=0.5, size=T).reshape(T, 1)
+        #y1 = np.random.binomial(n=10, p=0.5, size=T).reshape(T, 1)
+        #y2 = np.random.binomial(n=10, p=0.5, size=T).reshape(T, 1) +  x + y1
+        # z = np.random.binomial(n=10, p=0.5, size=T).reshape(T, 1)
+
+        test_stat, pval, dependent = ci.run_test_raw(x = x, y = np.hstack((y1,y2)), z = z, alpha_or_thres=0.91)#, z = z, alpha_pre = 0.5, pre_step_sample_fraction = 0.2, cond_ind_test = base_ci)# , x_type = np.ones((T, 1)), y_type = np.ones((T, 2)) , z_type = np.zeros((T, 1)))
+        print(dependent)
+        if (pval <= alpha):
+            rate[t] = 1
+        print("test_stat ", test_stat)
+        print("pval ", pval)
+    print(np.mean(rate))
+
diff --git a/tigramite/independence_tests/parcorr.py b/tigramite/independence_tests/parcorr.py
index 95c80e8b..ba4e65a5 100644
--- a/tigramite/independence_tests/parcorr.py
+++ b/tigramite/independence_tests/parcorr.py
@@ -35,6 +35,9 @@ class ParCorr(CondIndTest):
     For the ``significance='analytic'`` Student's-*t* distribution with
     :math:`T-D_Z-2` degrees of freedom is implemented.
 
+    Assumes one-dimensional X, Y.Use ParCorrMult for multivariate X, Y.
+
+
     Parameters
     ----------
     **kwargs :
diff --git a/tigramite/independence_tests/regressionCI.py b/tigramite/independence_tests/regressionCI.py
index 3ead1011..074f2e5c 100644
--- a/tigramite/independence_tests/regressionCI.py
+++ b/tigramite/independence_tests/regressionCI.py
@@ -17,7 +17,19 @@
 class RegressionCI(CondIndTest):
     r"""Flexible parametric conditional independence tests for continuous, categorical, or mixed data.
 
-    Assumes one-dimensional X, Y.
+    Asymptotically equivalent to the tests for mixed data suggested in
+
+    Tsagris, Michail, et al. "Constraint-based causal discovery with mixed
+    data." International journal of data science and analytics 6
+    (2018): 19-30.
+
+    For linear regression RegressionCI implements a likelihood ratio test,
+    while the above employs a F-statistic. Furthermore, while our
+    implementation utilizes the Chi^2 null distribution, theirs uses the
+    F-distribution. 
+
+    Assumes one-dimensional X, Y. But can be combined with PairwiseMultCI to
+    obtain a test for multivariate X, Y.
 
     Notes
     -----
@@ -105,10 +117,11 @@ def get_dependence_measure(self, array, xyz, data_type):
 
         def convert_to_one_hot(data, nb_classes):
             """Convert an iterable of indices to one-hot encoded labels."""
-            
-            targets = np.array(data).reshape(-1)
+            targets = np.array(data).reshape(-1).tolist()
             # categories need to be mapped to 0, 1, ... in this function
-            targets = targets - np.min(targets)
+            s = sorted(set(targets))
+            targets = [s.index(i) for i in targets]
+            targets = np.array(targets)
             return np.eye(nb_classes)[targets]
 
         def do_componentwise_one_hot_encoding(X, var_type):
diff --git a/tigramite/independence_tests/robust_parcorr.py b/tigramite/independence_tests/robust_parcorr.py
index ecb6c771..c21b3cd6 100644
--- a/tigramite/independence_tests/robust_parcorr.py
+++ b/tigramite/independence_tests/robust_parcorr.py
@@ -19,6 +19,9 @@ class RobustParCorr(CondIndTest):
         normal marginals, ordinary least squares (OLS) regression, and a test for
         non-zero linear Pearson correlation on the residuals.
 
+        Assumes one-dimensional X, Y. But can be combined with PairwiseMultCI to
+        obtain a test for multivariate X, Y.
+
         Notes
         -----
         To test :math:`X \perp Y | Z`, firstly, each marginal is transformed to be