Statistics with Pingouin

Pingouin is a statistics package in Python that wraps pandas and scipy-stats, creating more-complete statistics reports and presenting the results in a readable way.

Setup

Import Libraries

import numpy as np
import pandas as pd
import pingouin as pg
import seaborn as sns
import matplotlib.pyplot as plt
from pathlib import Path
import xarray as xr
import owncloud

Download the dataset and load it into a Pandas DataFrame

Path('data').mkdir(exist_ok=True, parents=True)

owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/3Uf2gScrvuTPQhB').get_file('/', f'data/steinmetz_2017-01-08_Muller.nc')

True

Section 1: T-Tests

T-tests compare the means of two samples of data generated from a normally-distributed population and compute the probability that they have the same mean.

Let’s do some analysis on some fake data to get a feel for these two statistics tools.

Generate the Data: Run the code below to create the dataset df.

Exercises

randomizer = np.random.RandomState(17)  # Makes sure the pseudorandom number generators reproduce the same data for us all.
df = pd.DataFrame()
df['A'] = randomizer.normal(0, 1, size=20)
df['B'] = randomizer.normal(0.2, 1, size=20)
df['C'] = randomizer.normal(0.7, 1, size=20)
df['D'] = df.A * 0.3 + randomizer.normal(0, 0.3, size=20)
df.plot.box()

Analyze the Data with T-Tests in Pingouin

Example: A vs 0, One-Sampled T-Test: Is the mean of the normally-distributed population that the dataset A is generated from unlikely to be zero?

pg.ttest(df['A'], 0)

Exercise: B vs 1, One-Sampled T-Test: Is the mean of the normally-distributed population that the dataset B is generated from unlikely to be one?

pg.ttest(df['B'], 1)

A vs B, Independent Samples T-Test: Is the mean of the normally-distributed population that the dataset A is generated from unlikely to be the same as the mean of the normally-distributed population that the dataset B is generated from?

pg.ttest(df['A'], df['B'])

Exercise: A vs C, Independent Samples T-Test: Is the mean of the normally-distributed population that the dataset A is generated from unlikely to be the same as the mean of the normally-distributed population that the dataset C is generated from?

pg.ttest(df['A'], df['C'])

Exercise: A vs C, Paired Samples T-Test: Is the mean of the differences between each pair of samples generated from the two normally-distributed populations A and C unlikely to be 0?

pg.ttest(df['A'], df['C'], paired=True)

Exercise: A vs D, Paired Samples T-Test: Is the mean of the differences between each pair of samples generated from the two normally-distributed populations A and D unlikely to be 0?

pg.ttest(df['A'], df['D'], paired=True)

Section 2: T-Tests on Real-World Data: Pupil Position Analysis

In the experiments reported by Steinmetz et al, 2019 in Nature, mice performed a discrimination task where they moved the position of a stimulus using a steering wheel. During the experiment, a camera recorded the pupil position of the subject in the x and y directions.

Load and Filter the Data

Run the code below, which loads the file, extracts the variables, filters the data, and transforms it into a dataframe.

Exercises

dset = xr.load_dataset('data/steinmetz_2017-01-08_Muller.nc')
dd = dset[['pupil_x', 'pupil_y', 'contrast_left', 'contrast_right', 'response_type']]
dd = dd.where(
    (dset.active_trials==True) 
    & (dset.response_type != 0)
    & (dset.time >= 0.5)
    & (dset.time <= 1.3)
    , drop=True)
dd['contrast_diff'] = dd['contrast_right'] - dd['contrast_left']
trials = dd.median(dim='time')
rdf = trials.to_dataframe()
rdf.head()

Extract data from the dataframe rdf, plot its distribution(s), and use a t-test to answer the questions below:

Example: For the trials where contrast_left was greater or equal to 50 and contrast_right was 0, was the average pupil_x position for the trial zero?

left = rdf[(rdf.contrast_left >= 50) & (rdf.contrast_right == 0)]
sns.kdeplot(left.pupil_x, label='left')
plt.axvline(0, color='black', linestyle='dotted')
plt.legend()
pg.ttest(left.pupil_x, 0)

Exercise: Where the mice looking at the stimulus that was presented? In trials where there was only one stimulus present (i.e. contrast_left or contrast_right was zero), and the stimulus had at least a contrast level of 50 (quite visible), was there a difference in the average pupil_x position?

Solution

mask_left = np.logical_and(rdf['contrast_left'] >= 50, rdf['contrast_right'] == 0)
df_left = rdf[mask_left]
sns.kdeplot(df_left['pupil_x'], label = 'left')

mask_right = np.logical_and(rdf['contrast_right'] >= 50, rdf['contrast_left'] == 0)
df_right = rdf[mask_right]
sns.kdeplot(df_right['pupil_x'], label = 'right')

plt.legend()

pg.ttest(df_left['pupil_x'], df_right['pupil_x'], paired = True)

	T	dof	alternative	p-val	CI95%	cohen-d	BF10	power
T-test	2.900766	41.883443	two-sided	0.005908	[0.19, 1.07]	0.855661	7.42	0.784252

Exercise: As a sanity check, was there a difference in vertical eye position between the left-stimulus trials and the right-stimulus trials? In trials where there was only one stimulus present (i.e. contrast_left or contrast_right was zero), and the stimulus had at least a contrast level of 50 (quite visible), was there a difference in the average pupil_y position?

Solution

mask_left = np.logical_and(rdf['contrast_left'] >= 50, rdf['contrast_right'] == 0)
df_left = rdf[mask_left]
sns.kdeplot(df_left['pupil_y'], label = 'left')

mask_right = np.logical_and(rdf['contrast_right'] >= 50, rdf['contrast_left'] == 0)
df_right = rdf[mask_right]
sns.kdeplot(df_right['pupil_y'], label = 'right')
plt.legend()

pg.ttest(df_left['pupil_y'], df_right['pupil_y'], paired = False)

	T	dof	alternative	p-val	CI95%	cohen-d	BF10	power
T-test	0.81158	38.482526	two-sided	0.422024	[-0.2, 0.47]	0.247714	0.39	0.125031

Exercise: Did mice look in the direction they turned the wheel? Was there a difference in average pupil_x position between left-turning trials (i.e. response_type is -1) and right-turning trials (i.e. response_type is 1)?

Solution

mask_left = rdf['response_type'] == -1
df_left = rdf[mask_left]
sns.kdeplot(df_left['pupil_x'], label = 'left')

mask_right = rdf['response_type'] == 1
df_right = rdf[mask_right]
sns.kdeplot(df_right['pupil_x'], label = 'right')

plt.legend()

pg.ttest(df_left['pupil_x'], df_right['pupil_x'], paired = False)

	T	dof	alternative	p-val	CI95%	cohen-d	BF10	power
T-test	-2.472363	78.49085	two-sided	0.015585	[-0.59, -0.06]	0.50892	3.159	0.610308

Exercise: As a sanity check, did vertical eye position matter? Was there a difference in average pupil_y position between left-turning trials (i.e. response_type is -1) and right-turning trials (i.e. response_type is 1)?

Solution

mask_left = rdf['response_type'] == -1
df_left = rdf[mask_left]
sns.kdeplot(df_left['pupil_y'], label = 'left')
mask_right = rdf['response_type'] == 1
df_right = rdf[mask_right]
sns.kdeplot(df_right['pupil_y'], label = 'right')
plt.legend()
pg.ttest(df_left['pupil_y'], df_right['pupil_y'], paired = False)

	T	dof	alternative	p-val	CI95%	cohen-d	BF10	power
T-test	1.218473	61.22414	two-sided	0.227719	[-0.1, 0.42]	0.275322	0.443	0.22797

Section 3: Working with Long DataFrames

If data is organized as a long dataframe, it is also possible to have pingouin simply do all the comparisons between the groups. Not only is this very convenient; it also makes it possible for pingouin to do richer analyses.

Run the cell below to melt the DataFrame df into a long DataFrame dfl to be used in the following exercises.

Exercises

dfl = df.melt(ignore_index=False).reset_index()
dfl.sample(5)

Example: Before doing pairwise tests, it’s often helpful to run an ANOVA analysis, which can tell you quickly whether a pairwise difference exists within a collection of potential pairs. Below, run an ANOVA test comparing the dv ‘value’ between values of ‘variable’ to see if a significant difference is expected.

pg.anova(dfl, dv='value', between=['variable'], detailed=True, )

Exercise: Do a t-tests for every possible pairwise comparison of “variable” (e.g. “A” vs “B”, “A” vs “C”, etc). (Tip: it works similar to the pg.anova function.)

pg.pairwise_tests(dfl, dv = 'value', between = ['variable'])

Exercise: Repeat the tests, but this time, adjust the p-value to compensate for multiple comparisons using the ‘fdr’ method. Are the significant differences still significant? How do the corrected p-values compare to the uncorrected p-values?

pg.pairwise_tests(dfl, dv = 'value', between = ['variable'], padjust = 'fdr_bh')

Exercise: Do a pairwise test again on the long dataframe, comparing the value of ‘value’ between each value ‘variable’. This time, use a non-parametric comparison (not assuming a normal distribution of each variable), and still adjust for multiple comparisons. Are your results the same?

pg.pairwise_tests(dfl, dv = 'value', between = ['variable'], padjust = 'fdr_bh', parametric = False)

Use the rdf dataframe from the pupil measurement data. Is there a significant different in average horizontal position, when comparing all the contrast_diff levels in the study pairwise?

pg.pairwise_tests(rdf, dv = 'pupil_x', between = ['contrast_diff'], padjust = 'fdr_bh', parametric = False)

Section 4: A Gallery of Plots : 2D Histograms with sns.jointplot()

Let’s try out a few different plots, for comparing two different continuous variables’ distributions.

Exercises

Example: Make a joint plot comparing the pupil_x and pupil_y variables.

sns.jointplot(rdf, x="pupil_x", y="pupil_y");

Exercise: Make a joint plot comparing the pupil_x and pupil_y variables, using kind='kde'.

Solution

sns.jointplot(rdf, x='pupil_x', y='pupil_y', kind = 'kde')

Exercise: Make a joint plot comparing the pupil_x and pupil_y variables, using kind='kde' and setting hue to 'response_type'.

Solution

sns.jointplot(rdf, x='pupil_x', y='pupil_y', kind = 'kde', hue='response_type')

Exercise: Make a joint plot comparing the pupil_x and pupil_y variables, using kind='hex'.

Solution

sns.jointplot(rdf, x='pupil_x', y='pupil_y', kind='hex')

Section 5: Further Reading

Nice article on Pingouin here: https://towardsdatascience.com/the-new-kid-on-the-statistics-in-python-block-pingouin-6b353a1db57c

Nice summary of the different effect size metrics and when to pick which: https://www.socscistatistics.com/effectsize/default3.aspx

Section 6: Extra: Correlations and Simpson’s Paradox

Run the cell below to download the iris dataset for this section:

Exercises

from bokeh.sampledata import iris
flowers = iris.flowers
flowers.head()

Example: Is there a correlation between flowers’ sepal widths and their petal lengths?

pg.corr(flowers['sepal_width'], flowers['petal_length'])

Exercise: Use seaborn’s lmplot() function to plot a regression line showing the sepal widths and petal lengths. Does the direction of this line match the results? Looking closely at the data’s distribution, what is missing from this story?

Solution

sns.lmplot(data=flowers, x='sepal_width', y = 'petal_length')

Exercise: This time factor in the species variable into the lmplot() visualization. What kind of relationship do you see between sepal_width and petal_length now?

Solution

sns.lmplot(data=flowers, x='sepal_width', y='petal_length', hue='species')

Demo

Example: What’s happening here is that a third variable, "species", is influencing both the petal and the sepal shape; this effect is called “Simpson’s Paradox”. Unfortunately, at this time (as far as I know) pengouin doesn’t have a way to split a correlation into multiple groups; one way to do this, though, is to use a df.groupby() analysis, concatenating the results together.

analyses = []
for name, group in flowers.groupby('species'):
    analysis = pg.corr(group['sepal_width'], group['petal_length'])
    analysis['name'] = name
    analyses.append(analysis)
results = pd.concat(analyses)
results

To summarize, accounting for the right variables is a big part of science; we should always be aware that the patterns we see in our data can be turned completely upside down when new aspects of a problem are understood and factored in!