# 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);


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")


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")

[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...
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: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 3 seconds.
Auto-assigning NUTS sampler...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Intercept, distance, direction_sigma]

100.00% [8000/8000 00:01<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.

[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.686 0.325 1.118 2.304 0.008 0.006 1966.0 1640.0 1.0
distance -0.010 0.004 -0.019 -0.003 0.000 0.000 2205.0 2290.0 1.0
direction_kappa 2.622 0.623 1.538 3.837 0.012 0.008 2871.0 2219.0 1.0
[8]:

_, ax = plt.subplots(1,2, figsize=(8, 4), sharey=True)
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_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");


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 = idata_vm.posterior_predictive["direction"].stack(samples=("chain", "draw")).T[::50]
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([])