Multi-sample comparisons in Python

Required packages

The required packages for this section are pandas, scipy and statsmodels. These can be installed with the following command in the command window (Windows) / terminal (Mac).

pip install pandas statsmodels scipy scikit-posthocs

Loading the data

We will again use the demo lipidomics dataset:

348KB

Lipidomics_dataset.xlsx

Load the dataset into a Pandas DataFrame named df as described in the basic plotting section:

import pandas as pd
df = pd.read_excel("Lipidomics_dataset.xlsx")
df.set_index("Sample Name", inplace=True)

Test for normality

Statsmodels contains a Lilliefors’ test for normality testing, which is a Kolmogorov-Smirnov test with estimated parameters. The function can be used as follows, for example on species PC 34:1 for N and T labelled samples:

from statsmodels.stats.diagnostic import kstest_normal

lipid = "PC 34:1"
group1 = df[df.Label == "T"][lipid]
group2 = df[df.Label == "PAN"][lipid]
group3 = df[df.Label == "N"][lipid]

print(kstest_normal(group1, "norm"))
print(kstest_normal(group2, "norm"))
print(kstest_normal(group3, "norm"))

The function returns 2 values: Kolmogorov-Smirnov test statistic and a pvalue; if the pvalue is lower than some threshold, e.g. 0.05, then we can reject the Null hypothesis that the sample comes from a normal distribution.

group T result: (0.08189149328478457, 0.08119506735189941)
group PAN result: (0.11146474201970102, 0.702986809787123)
group N result: (0.08643554971567619, 0.078397863558115)

Since the p-values is >0.05 for all groups, the distributions do not deviate significantly from the normal distribution and we assume normality.

ANOVA test

We'll perform an ANOVA test for the T, PAN and N labelled samples for species PC 34:1. We'll define our model with the ols function (ordinary least squares) from statsmodels, and then we'll pass this model to anova_lm function from statsmodels. Since the ols function doesn't accept spaces or special characters in our species names, we'll replace these with underscores first:

from statsmodels.formula.api import ols
from statsmodels.api import stats

df2 = df.copy()
df2.columns = df2.columns.str.rstrip()
df2.columns = df2.columns.str.replace(' ', '_')
df2.columns = df2.columns.str.replace(':', '_')

model = ols("PC_34_1 ~ Label", data=df2).fit()
ow_anova_table = stats.anova_lm(model, typ=2)
print(ow_anova_table) 

The function returns the following result:

                sum_sq     df          F    PR(>F)
Label     3.197941e+05    2.0  10.250407  0.000055
Residual  3.494196e+06  224.0        NaN       NaN

The One Way ANOVA Test returns a p < .05 which indicates that H0 can be rejected. Accordingly there is a significant difference between the mean of at least one group compared to other groups.

Levene's Homogenity Test

Next, we'll use Levene's Homogenity Test to determine wether the different groups can be assumed to have similar variances:

from scipy.stats import levene

levene(df.loc[df.Label == "T", "PC 34:1"],
       df.loc[df.Label == "PAN", "PC 34:1"],
       df.loc[df.Label == "N", "PC 34:1"])

This function returns the following result:

LeveneResult(statistic=6.072591927982112, pvalue=0.0027022136179025353)

Since Levene's Homogenity Test resulted in a p < .05, this indicates that H0 can be rejected. Accordingly there is a significant difference between the variance of at least one group compared to other groups. This means we can proceed with a Tamhane T2 test. Had the levene's test had a p > .05, H0 could not have been rejected and we could not have assumed a difference between the variances of the groups. In that case we would have used a post hoc Tukey test. We illustrate both cases bellow (keeping in mind the Tamhane T2 test would be the more appropriate test in this case)

Unequal Variance: Tamhane T2 test

We can use the Tamhane T2 test from scikit-posthocs:

import scikit_posthocs as sp

tamhane_table = sp.posthoc_tamhane(df, val_col="PC 34:1", group_col="Label")
print(tamhane_table) 

Which gives us the result:

            N       PAN         T
N    1.000000  0.001194  0.121347
PAN  0.001194  1.000000  0.019051
T    0.121347  0.019051  1.000000

We see that PAN vs N and PAN vs T have a p > 0.05, while T vs N does not.

Equal Variance: Post Hoc Tukey Test

We can use the Tamhane T2 test from scikit-posthocs:

import scikit_posthocs as sp

tukey_table = sp.posthoc_tukey(df, val_col="PC 34:1", group_col="Label")
print(tukey_table) 

Which gives us the result:

            N       PAN         T
N    1.000000  0.001000  0.121107
PAN  0.001000  1.000000  0.002523
T    0.121107  0.002523  1.000000

We see that PAN vs N and PAN vs T have a p > 0.05, while T vs N does not.

Kruskal-Wallis

In case the data is not normally distributed, the non parametric Kruskal-Wallis test is used instead of the ANOVA test. Let's perform another normality test, this time for species PC 36:5:

from statsmodels.stats.diagnostic import kstest_normal

lipid = "PC 36:5"
group1 = df[df.Label == "T"][lipid]
group2 = df[df.Label == "PAN"][lipid]
group3 = df[df.Label == "N"][lipid]

print(kstest_normal(group1, "norm"))
print(kstest_normal(group2, "norm"))
print(kstest_normal(group3, "norm"))

Which returns the results:

(0.12467614039830022, 0.0009999999999998899)
(0.09264557901905313, 0.9091272492376122)
(0.09253781472904482, 0.0448171307073585)

Since the p-values is <0.05 for 2 groups, these 2 distributions deviate significantly from the normal distribution and we can not assume normality. We will proceed with the non parametric Kruskal-Wallis test (we reuse the group1-3 variables calculated above for the normality test):

from scipy.stats import kruskal

kruskal(group1, group2, group3)

Which return the results:

KruskalResult(statistic=54.7072751290757, pvalue=1.319672436175926e-12)

The Kruskal-Wallis test resulted in a p < .05 which indicates that H0 can be rejected. Accordingly there is a significant difference between the distribution of at least one group compared to other groups. This means we can perform a Dunn posthoc test.

import scikit_posthocs as sp
dunn_table = sp.posthoc_dunn(df, val_col="PC 34:1", group_col="Label")
print(dunn_table)

Which gives us the result:

            N       PAN         T
N    1.000000  0.000104  0.113663
PAN  0.000104  1.000000  0.002757
T    0.113663  0.002757  1.000000

We see that PAN vs N and PAN vs T have a p > 0.05, while T vs N does not.