# 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:

{% file src="/files/O8GcvqGYOzpXcfU7NayN" %}

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

```python
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:

<pre class="language-python"><code class="lang-python"><strong>from statsmodels.stats.diagnostic import kstest_normal
</strong>
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"))
</code></pre>

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:

```python
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:

```python
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:

```python
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:

```python
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:

<pre class="language-python"><code class="lang-python"><strong>from statsmodels.stats.diagnostic import kstest_normal
</strong>
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"))
</code></pre>

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

```python
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.

```python
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.


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://laboratory-of-lipid-metabolism-a.gitbook.io/omics-data-visualization-in-r-and-python/metabolites-and-lipids-univariate-statistical-analysis-in-python/multi-sample-comparisons-in-python.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.