Circular Regression#

[1]:
import arviz as az
import bambi as bmb
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
[2]:
az.style.use("arviz-white")

Directional statistics, also known as circular statistics or spherical statistics, refers to a branch of statistics dealing with data which domain is the unit circle, as opposed to “linear” data which support is the real line. Circular data is convenient when dealing with directions or rotations. Some examples include temporal periods like hours or days, compass directions, dihedral angles in biomolecules, etc.

The fact that a Sunday can be both the day before or after a Monday, or that 0 is a “better average” for 2 and 358 degrees than 180 are illustrations that circular data and circular statistical methods are better equipped to deal with this kind of problem than the more familiar methods 1.

There are a few circular distributions, one of them is the VonMises distribution, that we can think as the cousin of the Gaussian that lives in circular space. The domain of this distribution is any interval of length \(2\pi\). We are going to adopt the convention that the interval goes from \(-\pi\) to \(\pi\), so for example 0 radians is the same as \(2\pi\). The VonMises is defined using two parameters, the mean \(\mu\) (the circular mean) and the concentration \(\kappa\), with \(\frac{1}{\kappa}\) being analogue of the variance. Let see a few example of the VonMises family:

[3]:
x = np.linspace(-np.pi, np.pi, 200)
mus = [0., 0., 0.,  -2.5]
kappas = [.001, 0.5,  3, 0.5]
for mu, kappa in zip(mus, kappas):
    pdf = stats.vonmises.pdf(x, kappa, loc=mu)
    plt.plot(x, pdf, label=r'$\mu$ = {}, $\kappa$ = {}'.format(mu, kappa))
plt.yticks([])
plt.legend(loc=1);
../_images/notebooks_circular_regression_4_0.png

When doing linear regression a commonly used link function is \(2 \arctan(u)\) this ensure that values over the real line are mapped into the interval \([-\pi, \pi]\)

[4]:
u = np.linspace(-12, 12, 200)
plt.plot(u, 2*np.arctan(u))
plt.xlabel("Reals")
plt.ylabel("Radians");
../_images/notebooks_circular_regression_6_0.png

Bambi supports circular regression with the VonMises family, to exemplify this we are going to use a dataset from the following experiment. 31 periwinkles (a kind of sea snail) were removed from it original place and released down shore. Then, our task is to model the direction of motion as function of the distance travelled by them after being release.

[5]:
data = bmb.load_data("periwinkles")
data.head()
[5]:
distance direction
0 107 1.169371
1 46 1.151917
2 33 1.291544
3 67 1.064651
4 122 1.012291

Just to compare results, we are going to use the VonMises family and the normal (default) family.

[6]:
model_vm = bmb.Model("direction ~ distance", data, family="vonmises")
idata_vm = model_vm.fit(include_mean=True)

model_n = bmb.Model("direction ~ distance", data)
idata_n = model_n.fit(include_mean=True)
WARNING (aesara.tensor.basic_opt): Optimization Warning: The Op i1 does not provide a C implementation. As well as being potentially slow, this also disables loop fusion.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
WARNING (aesara.tensor.basic_opt): Optimization Warning: The Op i0 does not provide a C implementation. As well as being potentially slow, this also disables loop fusion.
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Intercept, distance, direction_kappa]
100.00% [8000/8000 00:03<00:00 Sampling 4 chains, 0 divergences]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Intercept, distance, direction_sigma]
100.00% [8000/8000 00:02<00:00 Sampling 4 chains, 0 divergences]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.
The acceptance probability does not match the target. It is 0.8855, but should be close to 0.8. Try to increase the number of tuning steps.
[7]:
az.summary(idata_vm, var_names=["~direction_mean"])
[7]:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 1.692 0.321 1.106 2.299 0.008 0.006 1888.0 1722.0 1.0
distance -0.011 0.004 -0.019 -0.002 0.000 0.000 2212.0 2012.0 1.0
direction_kappa 2.630 0.616 1.496 3.749 0.014 0.010 2066.0 1835.0 1.0
[8]:
_, ax = plt.subplots(1,2, figsize=(8, 4), sharey=True)
posterior_mean = bmb.families.link.tan_2(idata_vm.posterior["direction_mean"])
ax[0].plot(data.distance, posterior_mean.mean(("chain", "draw")))
az.plot_hdi(data.distance, posterior_mean, ax=ax[0])

ax[0].plot(data.distance, data.direction, "k.")
ax[0].set_xlabel("Distance travelled (in m)")
ax[0].set_ylabel("Direction of travel (radians)")
ax[0].set_title("VonMises Family")

posterior_mean = idata_n.posterior["direction_mean"]
ax[1].plot(data.distance, posterior_mean.mean(("chain", "draw")))
az.plot_hdi(data.distance, posterior_mean, ax=ax[1])

ax[1].plot(data.distance, data.direction, "k.")
ax[1].set_xlabel("Distance travelled (in m)")
ax[1].set_title("Normal Family");
../_images/notebooks_circular_regression_12_0.png

We can see that there is a negative relationship between distance and direction. This could be explained as Periwinkles travelling in a direction towards the sea travelled shorter distances than those travelling in directions away from it. From a biological perspective, this could have been due to a propensity of the periwinkles to stop moving once they are close to the sea.

We can also see that if inadvertently we had assumed a normal response we would have obtained a fit with higher uncertainty and more importantly the wrong sign for the relationship.

As a last step for this example we are going to do a posterior predictive check. In the figure below we have to panels showing the same data, with the only difference that the on the right is using a polar projection and the KDE are computing taking into account the circularity of the data.

We can see that our modeling is failing at capturing the bimodality in the data (with mode around 1.6 and \(\pm \pi\)) and hence the predicted distribution is wider and with a mean closer to \(\pm \pi\).

[9]:
fig = plt.figure(figsize=(12, 5))
ax0 = plt.subplot(121)
ax1 = plt.subplot(122, projection='polar')

model_vm.predict(idata_vm, kind="pps")
pp_samples = az.extract_dataset(idata_vm, group="posterior_predictive", num_samples=200)["direction"]
colors = ["C0" , "k", "C1"]

for ax, circ in zip((ax0, ax1), (False, "radians", colors)):
    for s in pp_samples:
        az.plot_kde(s.values,  plot_kwargs={"color":colors[0], "alpha": 0.25}, is_circular=circ, ax=ax)
    az.plot_kde(idata_vm.observed_data["direction"].values,
                plot_kwargs={"color":colors[1], "lw":3}, is_circular=circ, ax=ax)
    az.plot_kde(idata_vm.posterior_predictive["direction"].values,
                plot_kwargs={"color":colors[2], "ls":"--", "lw":3}, is_circular=circ, ax=ax)

custom_lines = [Line2D([0], [0], color=c) for c in colors]

ax0.legend(custom_lines, ["posterior_predictive", "Observed", 'mean posterior predictive'])
ax0.set_yticks([])
fig.suptitle("Directions (radians)", fontsize=18);
../_images/notebooks_circular_regression_14_0.png

We have shown an example of regression where the response variable is circular and the covariates are linear. This is sometimes refereed as linear-circular regression in order to distinguish it from other cases. Namely, when the response is linear and the covariates (or at least one of them) is circular the name circular-linear regression is often used. And when both covariates and the response variables are circular, we have a circular-circular regression. When the covariates are circular they are usually modelled with the help of sin and cosine functions. You can read more about this kind of regression and other circular statistical methods in the following books.

[10]:
%load_ext watermark
%watermark -n -u -v -iv -w
Last updated: Mon Jun 20 2022

Python implementation: CPython
Python version       : 3.9.7
IPython version      : 8.3.0

scipy     : 1.7.3
arviz     : 0.13.0.dev0
bambi     : 0.9.0
matplotlib: 3.5.1
sys       : 3.9.7 (default, Sep 16 2021, 13:09:58)
[GCC 7.5.0]
pandas    : 1.4.2
numpy     : 1.21.5

Watermark: 2.3.0

[ ]: