Circular-linear correlation for regressing onto circular phase variables based on the method established by: Kempter R, Leibold C, Buzsáki G, Diba K, Schmidt R. Quantifying circular-linear associations: hippocampal phase precession. J Neurosci Methods. 2012;207(1):113–24. doi:10.1016/j.jneumeth.2012.03.007.

kempter.py

"""
Circular-linear correlation for regressing onto circular phase variables 
based on the method established by:

    Kempter R, Leibold C, Buzsáki G, Diba K, Schmidt R. Quantifying 
        circular-linear associations: hippocampal phase precession. 
        J Neurosci Methods. 2012;207(1):113–24. 
        doi:10.1016/j.jneumeth.2012.03.007.

Implemented By: @jdmonaco (github) / @j_d_monaco (twitter)
Last Updated: 2022-02-19

This software is provided AS IS under the terms of the Open Source MIT License. 
See http://www.opensource.org/licenses/mit-license.ph
"""

import numpy as np
import scipy.optimize as opt
import scipy.special as sp


def phaseregress(x, theta):
    """
    Compute circular-linear regression of phase against e.g. distance.
    """
    return KempterPhaseCorrelation(x, theta).regress()


def circmean(theta, weights=None, binsize=None, axis=-1):
    """
    Mean resultant vector angle for a set of phase angles (or 
    binned angular data).

    Arguments:
    theta -- array of radian angle values
    weights -- optional weights (i.e., for binned angle counts)
    binsize -- optional bin size for bias correction of binned data
    axis -- axis for averaging the input angles

    Returns:
    array of mean resultant vector angles
    """
    if weights is None:
        y_bar = np.sin(theta).sum(axis=axis) / theta.shape[axis]
        x_bar = np.cos(theta).sum(axis=axis) / theta.shape[axis]
    else:
        totw = weights.sum(axis=axis)
        y_bar = np.sum(weights * np.sin(theta), axis=axis) / totw
        x_bar = np.sum(weights * np.cos(theta), axis=axis) / totw
        if binsize is not None:
            correction = binsize / (2 * np.sin(binsize / 2))
            y_bar *= correction
            x_bar *= correction

    return np.arctan2(y_bar, x_bar)


class KempterPhaseCorrelation(object):

    """
    Implementation of Kempter et al. (2012) circular-linear correlation.
    """

    def __init__(self, x, phase, maxslope=1.0):
        x, phase = np.atleast_1d(x, phase)
        valid = np.logical_and(np.isfinite(x), np.isfinite(phase))
        self.phi = phase[valid]
        self.x = x[valid]
        self.xnorm = (self.x - self.x.min()) / self.x.ptp()
        self.alim = (-maxslope, maxslope)

    def regress(self):
        """Compute the circular-linear correlation.

        Returns:
        (slope, intercept, correlation r, p-value) tuple
        """
        self.maximize_residuals()
        self.estimate_intercept()
        self.correlation_and_pvalue()
        return 2 * np.pi * self.a_hat, self.phi0, self.rho_hat, self.p_value

    def _R(self, a):
        res = self.phi - 2 * np.pi * self.xnorm * a
        return -1.0 * np.sqrt(np.cos(res).mean()**2 + np.sin(res).mean()**2)

    def maximize_residuals(self):
        """Equation (1): Maximize residuals to estimate slope."""
        res = opt.minimize_scalar(self._R, bounds=self.alim, method='bounded')
        if not res.success:
            print('Error: {}', res.message, error=True)
            self.a_hat = 0.0
            return
        self.a_hat = res.x / self.x.ptp()  # rescale to data range

    def estimate_intercept(self):
        """Equation (2): Estimate the intercept."""
        res = self.phi - 2 * np.pi * self.x * self.a_hat
        self.phi0 = circmean(res)

    def correlation_and_pvalue(self):
        """Equations (3-4): Calculate correlation strength and p-value."""
        sdphi = np.sin(self.phi - circmean(self.phi))
        sdphi2 = np.power(sdphi, 2)

        theta = 2 * np.pi * np.abs(self.a_hat) * self.x
        sdtheta = np.sin(theta - circmean(theta))
        sdtheta2 = np.power(sdtheta, 2)

        lambda20 = np.sum(sdphi2)
        lambda02 = np.sum(sdtheta2)
        lambda22 = np.sum(sdphi2 * sdtheta2)

        self.rho_hat = (sdphi * sdtheta).sum() / np.sqrt(lambda20 * lambda02)

        z = self.rho_hat * np.sqrt(lambda20 * lambda02 / lambda22)
        self.p_value = sp.erfc(np.abs(z) / np.sqrt(2))