# 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="<https://1939159422-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FG1KwOmfbAgmKed8MHinY%2Fuploads%2Fcovj94eGSgTMJHhMDhz4%2FMatrix_PDAC_16082023.xlsx?alt=media&token=0c0701da-3afc-45cb-b72b-1d4663d81469>" %}

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.