Add first version of permutation interence method in the model_selection module.

[LegrandNico]

This is my first trial for this method (contributing to #124 ). The function is not working as expected yet, I think this is due to the way the correlation coefficients are extracted (?). Also, it is not feature-complete (it is not possible to provide permutation idx yet).

[jameschapman19]

import numpy as np
from typing import Optional, Tuple, Dict, Union

from sklearn.utils import check_random_state
from tqdm import tqdm
from cca_zoo.models._cca_base import _CCA_Base


def permutation_test_score(
        estimator: _CCA_Base, X: np.ndarray, Y: np.ndarray, latent_dims: int = 1,
        n_perms: int = 1000, Z: Optional[np.ndarray] = None, W: Optional[np.ndarray] = None,
        sel: Optional[np.ndarray] = None, partial: bool = True,
        parameters: Optional[Dict] = None,
        random_state: Union[int, np.random.RandomState] = None,
) -> Tuple[float, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Permutation inference for canonical correlation analysis (CCA) _[1].

    This code is adapted from the Matlab function accompagning the paper:
    https://github.com/andersonwinkler/PermCCA/blob/6098d35da79618588b8763c5b4a519438703dba4/permcca.m#L131-L164

    Parameters
    ----------
    estimator : _CCA_Base
        The object to use to fit the data. This must be one of the CCA models from
        py:class:`cca_zoo-models` and implementing a `fit` method.
    Y : np.ndarray
        Left set of variables, size N by P.
    X : np.ndarray
        Right set of variables, size N by Q.
    latent_dims : int
        The number of latent dimensions infered during the model fitting. Defaults to
        `1`.
    n_perms : int
        An integer representing the number of permutations. Default is 1000 permutations.
    Z : np.ndarray
        (Optional) Nuisance variables for both (partial CCA) or left
        side (part CCA) only.
    W : np.ndarray
        (Optional) Nuisance variables for the right side only (bipartial CCA).
    sel : np.ndarray
        (Optional) Selection matrix or a selection vector, to use Theil's residuals
        instead of Huh-Jhun's projection. If specified as a vector, it can be made
        of integer indices or logicals. The R unselected rows of Z (S of W) must be full
        rank. Use -1 to randomly select N-R (or N-S) rows.
    partial : bool
        (Optional) Boolean indicating whether this is partial (true) or part (false) CCA.
        Default is true, i.e., partial CCA.
    parameters : dict | None
        (Optional) Any additional keyword arguments required by the given estimator.

    Returns
    -------
    p : float
        p-values, FWER corrected via closure.
    r : np.ndarray
        Canonical correlations.
    A : np.ndarray
        Canonical coefficients (X).
    B : np.ndarray
        Canonical coefficients (Y).
    U : np.ndarray
        Canonical variables (X).
    V : np.ndarray
        Canonical variables (Y).

    References
    ----------
    .. [1] Winkler AM, Renaud O, Smith SM, Nichols TE. Permutation Inference for
        Canonical Correlation Analysis. NeuroImage. 2020; 117065.
    """

    random_state = check_random_state(random_state)
    lW, cnt = np.zeros(latent_dims), np.zeros(latent_dims)
    n_obs = X.shape[0]
    if parameters is None:
        parameters = {}

    # Initial fit of the CCA model (without any permutation)
    init_model = estimator(latent_dims=latent_dims, random_state=random_state, **parameters).fit((X, Y))
    A, B = init_model.weights
    U, V = X @ np.hstack((A, null(A.T))), Y @ np.hstack((B, null(B.T)))
    x_idx = np.arange(n_obs)
    y_idx = np.arange(n_obs)
    for i in tqdm(range(n_perms)):

        # If user didn't supply a set of permutations, permute randomly both Y and X.
        if i > 0:
            random_state.shuffle(x_idx)
            random_state.shuffle(y_idx)

        # For each canonical variable
        for k in range(latent_dims):
            # Fit the CCA model using the permuted datasets
            perm_model = estimator(latent_dims=(latent_dims - k), random_state=random_state, **parameters)
            perm_model.fit((U[x_idx, k:], V[y_idx, k:]))

            # Estimate correlation coefficient for this CCA fit
            r_perm = perm_model.score((U[x_idx, k:], V[y_idx, k:]))

            lWtmp = -1 * np.cumsum(np.log(1 - r_perm ** 2)[::-1])[::-1]
            lW[k] = lWtmp[0]

        if i == 0:
            #copy otherwise lw1 and lW share memory
            lw1 = lW.copy()
        cnt = cnt + (lW >= lw1)

    # compute p-values
    p = np.maximum.accumulate(cnt / n_perms)

    return p, A, B, U, V


def null(a, rtol=1e-5):
    # https://stackoverflow.com/questions/19820921/a-simple-matlab-like-way-of-finding-the-null-space-of-a-small-matrix-in-numpy
    u, s, v = np.linalg.svd(a)
    rank = (s > rtol * s[0]).sum()
    return v[rank:].T.copy()

Really nice work translating the matlab. I think I've found the sources of the bugs you mention. Take a look and see what you think.

[jameschapman19]

from cca_zoo.model_selection import permutation_test_score
from sklearn.utils.validation import check_random_state
from cca_zoo.models import CCA
from cca_zoo.data import generate_covariance_data

n = 50
rng = check_random_state(0)

# this generates data with ground truth number of latent dimensions
(X, Y), (_, _) = generate_covariance_data(1000, [5, 5], latent_dims=4, correlation=[1, 1, 0, 0], random_state=rng)


def test_perm_test():
    p, A, B, U, V = permutation_test_score(X=X, Y=Y, estimator=CCA, latent_dims=4, random_state=rng)
    assert p[0] < 0.05
    assert p[1] < 0.05
    assert p[2] > 0.05
    assert p[3] > 0.05

Wonder if something like this might make a good test.

[jameschapman19]

This is really nice work!

Sorry for my slowness getting back to you. Had to properly read the MATLAB code to see how it works. Let me know your thoughts on the above. Would be good to try to side-by-side some data with the MATLAB code as a sense check.

[jameschapman19]

A particularly fiddly bit is that the permutations are on the U and V matrices which are subspaces of X and Y.