import math
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

import pingouin as pg

# delete boring Future warnings from pingouin
import warnings
warnings.simplefilter(action='ignore')

import sys
print("Python:", sys.version)
print("numpy:", np.__version__)
print("pandas:", pd.__version__)
print("seaborn:", sns.__version__)

Python: 3.12.2 (tags/v3.12.2:6abddd9, Feb  6 2024, 21:26:36) [MSC v.1937 64 bit (AMD64)]
numpy: 1.26.4
pandas: 2.2.1
seaborn: 0.13.2

!pip show pingouin

Name: pingouin
Version: 0.5.3
Summary: Pingouin: statistical package for Python
Home-page: https://pingouin-stats.org/index.html
Author: Raphael Vallat
Author-email: raphaelvallat9@gmail.com
License: GPL-3.0
Location: C:\Users\Sébastien\Documents\second_scight\website\analysis_animal_exp\.env\Lib\site-packages
Requires: matplotlib, numpy, outdated, pandas, pandas-flavor, scikit-learn, scipy, seaborn, statsmodels, tabulate
Required-by:

# Set the maximum number of rows to display
pd.set_option('display.max_rows', 15)

# import the raw data
data = pd.read_csv('./ifng_data.csv', sep=";", na_values='na')

# mapping the dose to the group label
map_groups = {
    'a': 0,
    'b': 1.6,
    'c': 8,
    'd': 40,
    'e': 200,
}

data['dose_µg'] = data['group'].map(map_groups)

# creating a new variable based on the group id and animal number
data['animal_id'] = data['group'] + '|' + data['group_#'].apply(str)
data = data.drop(columns=['group', 'group_#'])

# pseudonymisation
stim_map = {
    'aaaa': 'CMP1',
    'bbbb': 'CMP2',
    'cccc': 'PC',
    'medium': 'medium',
}

data['stim'] = data['stim'].map(stim_map)

# we also encode the stimulation condition for the ANOVA
data['stim_id'] = data['stim'] + '|' + data['conc_µg_ml-1'].apply(str)
data

ix = data['animal_id'] == 'a|5'
data[ix]

# drop animal A5 from the whole analysis
data.drop(data[ix].index, inplace=True)

data.info()

<class 'pandas.core.frame.DataFrame'>
Index: 640 entries, 0 to 655
Data columns (total 7 columns):
 #   Column        Non-Null Count  Dtype  
---  ------        --------------  -----  
 0   sample        640 non-null    object 
 1   stim          640 non-null    object 
 2   conc_µg_ml-1  560 non-null    float64
 3   ifng_pg_ml-1  640 non-null    float64
 4   dose_µg       640 non-null    float64
 5   animal_id     640 non-null    object 
 6   stim_id       640 non-null    object 
dtypes: float64(3), object(4)
memory usage: 40.0+ KB

data[data['conc_µg_ml-1'].isna()]

data_PC = data[data['stim'] == 'PC']

def draw_boxnstripplot(
        x='dose_µg',
        y='ifng_pg_ml-1',
        data=data,
        hue='sample',
        palette='Accent',
    ):
    """
    Draw a stripplot and a customized boxplot, with hue on sample (SP vs. LN).
    """

    ax = sns.boxplot(
        x=x,
        y=y,
        data=data,
        hue=hue,
        showcaps=False,
        showfliers=False,
        boxprops=dict(linewidth=2, edgecolor='k', facecolor='white',),
        whiskerprops=dict(linewidth=2, color='k',),
        capprops=dict(linewidth=2, color='k',),
        medianprops=dict(linewidth=2, color='darkgrey',),
    )

    sns.stripplot(
        x=x,
        y=y,
        data=data,
        hue=hue,
        dodge=True,
        palette=palette,
        size=5,
        alpha=0.85,
        ax=ax,
    )

    sns.despine()
    # generate a nicer legend
    handles, labels = ax.get_legend_handles_labels()
    n_items = len(labels)//2
    plt.legend(
        handles[n_items:], labels[n_items:],
        fontsize=8,
        title_fontsize=9,
        fancybox=False,
    )

plt.figure(figsize=(6,3.5))

draw_boxnstripplot(data=data_PC, palette='Dark2')

plt.ylim((5, 2e5))
plt.xlabel("dose of drug (µg)")
plt.ylabel("[IFN-$\gamma$](pg/mL)")
plt.yscale("log")
plt.title("restimulation with Positive Control")
plt.gca().get_legend().set_title("sample")
plt.tight_layout();

data_PC.groupby(['sample','dose_µg']).describe()['ifng_pg_ml-1']

# 2-way ANOVA
data_PC.anova(
    dv="ifng_pg_ml-1",
    between=['sample', 'dose_µg'],
).round(3)

pg.pairwise_tests(
    data=data_PC,
    dv="ifng_pg_ml-1",
    between=['sample', 'dose_µg'],
).round(3)

animals = pd.read_csv("./animals.csv", sep=';', comment='!')

day_exp = 44887 # day of the study
animals['age_days'] = day_exp - animals['day_of_birth'].astype(float)
animals['animal_id'] = animals['group'] + '|' + animals['group_#'].apply(str)

animals.drop(columns=['group_#', 'id', 'day_of_birth'], inplace=True)

animals.head()

animals.groupby(['group', 'sex'])['age_days'].describe()

plt.figure(figsize=(3,2))
sns.boxplot(
    x='group',
    y='age_days',
    data=animals,
    hue='sex',
    showfliers=True,
);
# plt.tight_layout();

pg.anova(data=animals, dv='age_days', between=['sex', 'group'])

# merge this data with the full ELISA data set
data_animals = data.merge(animals, on='animal_id')
data_animals.head(3)

sns.relplot(
    x='age_days',
    y='ifng_pg_ml-1',
    style='sex',
    hue='conc_µg_ml-1',
    size='dose_µg',
    sizes=[1, 8, 40, 90, 250],
    data=data_animals.fillna(10), # concentration for PC
    col='stim',
    row='sample',
    height=2.5,
    aspect=1,
    legend='auto',
    alpha=.75,
    palette='nipy_spectral',
)
plt.yscale('log')
# plt.tight_layout();

cytometry = pd.read_csv("cytometry.csv", sep=';', comment='!')
cytometry['animal_id'] = cytometry['group'] + '|' + cytometry['group_#'].apply(str)
data_cytometry = data.merge(cytometry, on='animal_id')
data_cytometry.head()

plt.figure(figsize=(3,2))
sns.boxplot(
    x='group',
    y='viability_pct',
    data=cytometry,
)

plt.ylim((80, 100));

# for celltype in ('hladr_pct', 'cd19_pct', 'cd8_pct', 'cd4_pct'):
for celltype in ('cd8_pct', 'cd4_pct'):
    sns.relplot(
        x=celltype,
        y='ifng_pg_ml-1',
        hue='conc_µg_ml-1',
        style='sample',
        size='dose_µg',
        sizes=[1, 6, 30, 80, 200],
        data=data_cytometry.fillna(10), # fake concentration for PC
        col='stim',
        height=3,
        aspect=.75,
        legend='auto',
        alpha=.70,
        palette='Set1',
    )
    plt.yscale('log');

# filtering data either CMP1 or medium control
data_CMP1 = data[(data['stim'] == 'CMP1') | (data['conc_µg_ml-1'] == 0)]

data_CMP1.groupby(['dose_µg', 'sample', 'stim', 'conc_µg_ml-1']).describe(percentiles=[.5])['ifng_pg_ml-1']

plt.figure(figsize=(10,4))

for i, sampl_ in enumerate(('sp', 'ln')):

    plt.subplot(1,2,i+1)
    draw_boxnstripplot(
        x='dose_µg',
        data=data_CMP1[data_CMP1['sample'] == sampl_],
        hue='conc_µg_ml-1',
        palette='tab20b_r',
    )

    plt.yscale('linear')
    # plt.ylim((0,25000))
    plt.xlabel('dose of drug (µg)')
    plt.ylabel('[IFN$\gamma$](pg/mL)')
    plt.title(f"influence of drug dose on {sampl_} samples restim with CMP1", size=10, fontweight='bold')
    plt.gca().get_legend().set_title("concentration")
plt.tight_layout();

data_CMP2 = data[(data['stim'] == 'CMP2') | (data['conc_µg_ml-1'] == 0)]

data_CMP2.groupby(['dose_µg', 'sample', 'stim', 'conc_µg_ml-1']).describe(percentiles=[.5])['ifng_pg_ml-1']

plt.figure(figsize=(10,4))

for i, sampl_ in enumerate(('sp', 'ln')):

    plt.subplot(1,2,i+1)
    draw_boxnstripplot(
        x='dose_µg',
        data=data_CMP2[data_CMP2['sample'] == sampl_],
        hue='conc_µg_ml-1',
        palette='tab20b_r',
    )

    plt.yscale('linear')
    # plt.ylim((0,25000))
    plt.xlabel('dose of drug (µg)')
    plt.ylabel('[IFN$\gamma$](pg/mL)')
    plt.title(f"influence drug dose on {sampl_} samples restim with CMP2", size=10, fontweight='bold')
    plt.gca().get_legend().set_title("concentration")
plt.tight_layout();

plt.figure(figsize=(10,7))
i=1

for stim_ in ('CMP1', 'CMP2'):
    for sampl_ in ('sp', 'ln'):
        plt.subplot(2,2,i)
        draw_boxnstripplot(
            x='conc_µg_ml-1', data=data[
                (
                (data['stim'] == stim_)
                |
                (data['conc_µg_ml-1'] == 0)
                )
                &
                (data['sample'] == sampl_)
            ],
            hue='dose_µg',
            palette=sns.color_palette('muted', 5)
        )

        plt.yscale('log')
        plt.ylim(5, 2e5)
        plt.xlabel(f"[{stim_}](µg/mL)")
        plt.ylabel('[IFN$\gamma$](pg/mL), $\log_{10}$')
        plt.title(f"{sampl_} samples restimulated with {stim_}", size=10, fontweight='bold')
        plt.gca().get_legend().set_title("dose (µg)")
        i+=1

plt.tight_layout();

data[(data['sample'] == 'sp') & (data['conc_µg_ml-1'] == 1)]['ifng_pg_ml-1'].max()

79350.82434

pd.set_option("display.max_rows", None, "display.max_columns", None)
res_norm = data.groupby(['sample', 'stim', 'conc_µg_ml-1', 'dose_µg'])['ifng_pg_ml-1'].apply(pg.normality)
res_norm.sample(6)

# we just count the proportion of `True` normality
pct_normal_distributed_lin = 100*res_norm['normal'].mean()
print(f"normal distribution founded in {pct_normal_distributed_lin:.1f}% of the box plots (original/untransformed)")

normal distribution founded in 65.7% of the box plots (original/untransformed)

data['log_ifng'] = data['ifng_pg_ml-1'].apply(np.log10)

res_norm_log = data.groupby(['sample', 'stim', 'conc_µg_ml-1', 'dose_µg'])['log_ifng'].apply(pg.normality)

print(f"normal distribution founded in {100*res_norm_log['normal'].mean():.1f}% of the box plots (log-transformed)")

normal distribution founded in 92.9% of the box plots (log-transformed)

plt.figure(figsize=(6,3))
plt.subplot(121)
pg.qqplot(data['ifng_pg_ml-1'])
plt.title("original")

plt.subplot(122)
pg.qqplot(data['log_ifng'])
plt.title("log-transformed")
plt.tight_layout();

data_filtered=(
data
.dropna()                   # remove PC
[data['conc_µg_ml-1'] != 0] # remove medium control
)

data_filtered.head()

# let's save now the data in a csv file for later use in R
data_filtered.to_csv("./data_filtered.csv")

plt.figure(figsize=(9,6))
i=1

for stim_ in ('CMP1', 'CMP2'):
    for sampl_ in ('sp', 'ln'):
        for conc_ in (0.1, 1.0, 10):
            for dose_ in data['dose_µg'].unique():
                data_ = data[
                        (data['stim'] == stim_)
                        &
                        (data['sample'] == sampl_)
                        &
                        (data['conc_µg_ml-1'] == conc_)
                        &
                        (data['dose_µg'] == dose_)
                    ]

                ax_ = plt.subplot(6,10,i)
                pg.qqplot(data_['ifng_pg_ml-1'], confidence=False, ax=ax_)
                plt.xlabel('')
                plt.ylabel('')
                for text_ in  ax_.texts: text_.set_fontsize(5) # decrease size of R² text
                plt.title(f"{stim_}|{sampl_}|{conc_}|{dose_}", fontsize=6)
                plt.xticks([],[]) #remove x-axis ticks
                plt.yticks([],[]) #remove y-axis ticks
                i+=1
plt.suptitle("individual QQ-plots for untransformed values")
plt.tight_layout();

plt.figure(figsize=(9,6))
i=1

for stim_ in ('CMP1', 'CMP2'):
    for sampl_ in ('sp', 'ln'):
        for conc_ in (0.1, 1.0, 10):
            for dose_ in data['dose_µg'].unique():
                data_ = data[
                        (data['stim'] == stim_)
                        &
                        (data['sample'] == sampl_)
                        &
                        (data['conc_µg_ml-1'] == conc_)
                        &
                        (data['dose_µg'] == dose_)
                    ]

                ax_ = plt.subplot(6,10,i)
                pg.qqplot(data_['log_ifng'], confidence=False, ax=ax_)
                plt.xlabel('')
                plt.ylabel('')
                for text_ in  ax_.texts: text_.set_fontsize(5) # decrease size of R² text
                plt.title(f"{stim_}|{sampl_}|{conc_}|{dose_}", fontsize=6)
                plt.xticks([],[]) #remove x-axis ticks
                plt.yticks([],[]) #remove y-axis ticks
                i+=1
plt.suptitle("individual QQ-plots for log-transformed values")
plt.tight_layout();

# untransformed data
data_filtered.groupby(['sample', 'stim_id']).apply(pg.homoscedasticity, dv='ifng_pg_ml-1', group='dose_µg')

# log-transformed data
data_filtered.groupby(['sample', 'stim_id']).apply(pg.homoscedasticity, dv='log_ifng', group='dose_µg')

pg.box_m(data_filtered, dvs=['ifng_pg_ml-1'], group='dose_µg')

pg.box_m(data_filtered, dvs=['log_ifng'], group='dose_µg')

library(tidyverse)
library(tidyr)
library(ggpubr)

df<-read_csv("data_filtered.csv")

bxp <- ggboxplot(
    df, x = "dose_µg", y = "log_ifng",
    color = "stim_id", palette = "jco",
    facet.by = "sample", short.panel.labs = FALSE
)
bxp

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
New names:
• `` -> `...1`
Rows: 480 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): sample, stim, animal_id, stim_id
dbl (5): ...1, conc_µg_ml-1, ifng_pg_ml-1, dose_µg, log_ifng

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

library(rstatix)

res.aov <- anova_test(
      data = df, dv = log_ifng, wid = animal_id,
      between = dose_µg, within = c(sample, stim_id)
)
get_anova_table(res.aov)

Attachement du package : 'rstatix'


L'objet suivant est masqué depuis 'package:stats':

    filter

# Two-way ANOVA at each dose group level
two.way <- df %>%
    group_by(sample) %>%
    anova_test(dv = log_ifng, between = dose_µg, wid = animal_id, within = stim_id)

get_anova_table(two.way)

# same in Python
data_filtered.groupby(['sample']).apply(
    pg.mixed_anova,
    dv='log_ifng',
    between='dose_µg',
    within='stim_id',
    subject='animal_id',
    effsize='ng2', # generalized eta-squared
    correction=False
)

df %>%
    group_by(sample, stim_id) %>%
    anova_test(dv=log_ifng, wid=animal_id, between=dose_µg) %>%
    get_anova_table()

# same in python
data_filtered.groupby(['sample', 'stim_id']).apply(
    pg.anova,
    dv='log_ifng',
    between='dose_µg',
)

df %>%
    group_by(sample, stim_id) %>%
    pairwise_t_test(
    log_ifng ~ dose_µg, paired = FALSE,
    p.adjust.method = "bonferroni"
    )

pwc = data_filtered.groupby(['sample', 'stim_id']).apply(
    pg.pairwise_tests,
    dv="log_ifng",
    between='dose_µg',
    padjust='fdr_bh',
    parametric=True,
    # effsize='none',
)

pwc.round(4).head(10) # note that we use the values from this table in the final figure

def asterix(ax=None, sample=None, stim=None, ticks=False, data=data):
    """
    Function
    --------
    draw_significance_line(data, pairwise_comparisons, ticks=False)

    Purpose
    -------
    This function automatically generates horizontal lines on a graph 
    to indicate statistically significant differences between groups 
    in a box or strip plot.

    Parameters
    ----------
    data
        A data structure containing the data for the box/strip plots
        (format depends on your plotting library).
    pairwise_comparisons
        A table containing the results of pairwise comparisons between 
        groups. This table should specify which comparisons are 
        statistically significant.
    ticks (optional)
        A boolean value indicating whether to add small ticks at the 
        ends of the horizontal lines (default: False).
    
    Functionality
    -------------
    The function iterates through all combinations of sample and stim 
    in the data. For each combination, it identifies statistically 
    significant differences based on the provided pairwise_comparisons 
    table. If a significant difference is found, the function draws a 
    horizontal line between the two groups on the graph. The height of 
    the line is positioned slightly above the highest data point in the 
    compared groups. Additionally, the function displays the level of 
    significance (number of asterisks) next to the line.
    
    Overall, this function simplifies the process of visually highlighting 
    significant differences in box/strip plots.
    """

    def n_stars(pval):
        """
        Encodes the number of asterices corresponding to P value.
        """
        if   pval <= 1e-4 : return '*'*4
        elif pval <= 1e-3 : return '*'*3
        elif pval <= 1e-2 : return '*'*2
        elif pval <= .05  : return '*'
        else              : return 'ns'
    
    def plot_significant(x1, x2, y, pval=1):
        """
        Actually plot the line and asterices.
        """

        ax.plot(
            [x1, x1, x2, x2],
            [y-(.04*ticks), y, y, y-(.04*ticks)],
            lw=1,
            color='k'
        )
        ax.text(
            (x1 + x2)*.5,
            y-.03,
            n_stars(pval),
            ha='center',
            va='bottom',
            color='red',
            fontdict={'size': 9},
        )

    # map the treatment/dose group to the x offset for each boxplot group
    x_offset = {
        0  : -.33,
        1.6: -.16,
        8  : 0,
        40 : .16,
        200: .33
    }

    # main loop of the function that goes thought each combination of 
    # stimulus and concentration, to extract and plot all the signficant differences
    for conc in (.1, 1, 10):
        y = round(data[data['conc_µg_ml-1'] == conc]['log_ifng'].max(), 1) + .15 # initial height of the line

        # extraction of all significant differences from the previous pairwise comparisons
        stim_id = data['stim_id'].unique()[0]
        pwc_ = pwc.loc[(sample, stim_id)]
        sign_ = pwc_[pwc_['p-corr'] < .05]
        
        # loop over the list and plot
        for idx in sign_.index:
            plot_significant(
                x1 = math.log10(conc) + 1 + x_offset[sign_.loc[idx, 'A']], # 0.1 on the axis is x=0
                x2 = math.log10(conc) + 1 + x_offset[sign_.loc[idx, 'B']], # 1.0 on the axis is x=1, etc.
                y=y,
                pval=sign_.loc[idx, 'p-corr']
            )
            y+=.2 # next lines will be higher

plt.figure(figsize=(10,8))
i=1

for stim_ in ('CMP1', 'CMP2'):
    for sampl_ in ('sp', 'ln'):
        data_ = data[(data['stim'] == stim_) & (data['sample'] == sampl_)]
        ax_ = plt.subplot(2,2,i)
        
        draw_boxnstripplot(
            x='conc_µg_ml-1',
            y='log_ifng',
            data=data_,
            hue='dose_µg',
            palette=sns.color_palette('Set2', 5)
        )

        plt.xlabel(f"[{stim_.upper()}](µg/mL)")
        plt.ylabel('[IFN$\gamma$](pg/mL), $\log_{10}$')
        plt.title(f"{sampl_.upper()} restimulated with {stim_.upper()}", size=10, fontweight='bold')
        ax_.get_legend().set_title("dose (µg)")

        # plot the significant differences for the active sample|stim combination and all concentrations
        asterix(ax=ax_, sample=sampl_, stim=stim_, data=data_, ticks=True)

        i+=1

plt.suptitle("3-way mixed design ANOVA, $F(13, 113.76)=2.609,P=0.003,\eta_g^2=0.014$", fontweight='bold', color='darkblue')
plt.annotate(
    text="pwc: unpaired T tests; p.adjust: Benjamini/Hochberg FDR",
    xy=(1,-.2),
    xycoords='axes fraction',
    color='darkred',
    horizontalalignment='right',
    verticalalignment='top')
plt.tight_layout();

	sample	stim	conc_µg_ml-1	ifng_pg_ml-1	dose_µg	animal_id	stim_id
7	sp	PC	NaN	2357.286815	0.0	a\|1	PC\|nan
15	sp	PC	NaN	4831.622047	0.0	a\|2	PC\|nan
23	sp	PC	NaN	7404.615045	0.0	a\|3	PC\|nan
31	sp	PC	NaN	6741.175581	0.0	a\|4	PC\|nan
47	sp	PC	NaN	7632.145795	0.0	a\|6	PC\|nan
...	...	...	...	...	...	...	...
623	ln	PC	NaN	14867.005840	200.0	e\|5	PC\|nan
631	ln	PC	NaN	890.046823	200.0	e\|6	PC\|nan
639	ln	PC	NaN	9169.993506	200.0	e\|7	PC\|nan
647	ln	PC	NaN	10161.485230	200.0	e\|8	PC\|nan
655	ln	PC	NaN	16477.947980	200.0	e\|9	PC\|nan

		count	mean	std	min	25%	50%	75%	max
sample	dose_µg
ln	0.0	7.0	15701.805526	6739.609791	7496.349090	11174.950206	15170.899100	18821.843260	27251.803560
	1.6	8.0	9536.057609	4977.493491	2711.783721	6669.294884	9774.753105	11924.078515	17960.493060
	8.0	8.0	9316.513183	6739.574502	1758.995678	5108.393259	8511.964890	10960.736607	23807.761460
	40.0	8.0	10142.228410	7894.074925	3422.978052	4126.810321	7569.684480	13610.679860	26195.707890
	200.0	9.0	14583.158003	7392.203196	890.046823	10161.485230	14867.005840	17700.191070	25559.543020
sp	0.0	7.0	6521.857562	2285.459666	2357.286815	5786.398814	7237.272611	7518.380420	9448.885043
	1.6	8.0	2477.304483	1093.774661	390.289530	2118.675463	2810.950628	3172.089082	3726.386554
	8.0	8.0	2615.217522	4216.415500	463.763273	653.130086	1158.616933	1949.494963	12950.361210
	40.0	8.0	3554.253922	4883.213734	398.796385	774.024033	1927.244999	3383.356395	15098.226550
	200.0	9.0	3724.867523	2389.856297	641.664264	2640.331036	2940.849040	4762.757852	8277.768827

	Source	SS	DF	MS	F	p-unc	np2
0	sample	1.318432e+09	1.0	1.318432e+09	45.823	0.000	0.396
1	dose_µg	3.062497e+08	4.0	7.656243e+07	2.661	0.040	0.132
2	sample * dose_µg	5.961853e+07	4.0	1.490463e+07	0.518	0.723	0.029
3	Residual	2.014061e+09	70.0	2.877230e+07	NaN	NaN	NaN

	Contrast	sample	A	B	Paired	Parametric	T	dof	alternative	p-unc	BF10	hedges
0	sample	-	ln	sp	False	True	6.573	78.000	two-sided	0.000	1.818e+06	1.456
1	dose_µg	-	0.0	1.6	False	True	2.311	23.776	two-sided	0.030	2.369	0.840
2	dose_µg	-	0.0	8.0	False	True	2.122	27.018	two-sided	0.043	1.774	0.758
3	dose_µg	-	0.0	40.0	False	True	1.669	27.822	two-sided	0.106	0.967	0.592
4	dose_µg	-	0.0	200.0	False	True	0.762	29.492	two-sided	0.452	0.422	0.260
...	...	...	...	...	...	...	...	...	...	...	...	...
26	sample * dose_µg	sp	1.6	40.0	False	True	-0.609	14.000	two-sided	0.552	0.485	-0.288
27	sample * dose_µg	sp	1.6	200.0	False	True	-1.409	11.486	two-sided	0.185	0.808	-0.624
28	sample * dose_µg	sp	8.0	40.0	False	True	-0.412	14.000	two-sided	0.687	0.453	-0.195
29	sample * dose_µg	sp	8.0	200.0	False	True	-0.657	10.798	two-sided	0.525	0.487	-0.313
30	sample * dose_µg	sp	40.0	200.0	False	True	-0.090	9.905	two-sided	0.930	0.421	-0.043

	group	sex	cage	termination	age_days	animal_id
0	a	m	12	ES	59.0	a\|1
1	a	m	12	ES	59.0	a\|2
2	a	m	16	SA	94.0	a\|3
3	a	m	16	SA	94.0	a\|4
4	a	m	17	SA	63.0	a\|5

Case study | Analysis of an in vivo experiment¶

Introduction¶

Experimental details¶

Preparation of the raw data¶

Importing libraries and modules¶

Load data¶

Data cleaning¶

Sanity check¶

QC analysis¶

Visualization function¶

Summary chart¶

Summary statistics¶

Combining data sources for comprehensive analysis¶

Animal data¶

Flow cytometry¶

Systematic overview and preanalysis¶

Cells restimulated with CMP1¶

CMP2-stimulated cells¶

Assumptions for statistical analysis¶

Data exploration on a logarithmic scale¶

Normality tests¶

Applying Log transformation¶

QQ plots¶

Homoscedasticity of variances¶

Homogeneity of covariances¶

Sphericity¶

Mixed ANOVA¶

Combining R analysis with the Notebook¶

Post-hoc tests¶

Differences between R and Python packages¶

Simple Main Effect analysis¶

False Discovery Rate adjustment¶

Finalized figure¶

	sample	stim	conc_µg_ml-1	ifng_pg_ml-1	dose_µg	animal_id	stim_id
0	sp	medium	0.0	59.405353	0.0	a\|1	medium\|0.0
1	sp	CMP1	10.0	5671.825154	0.0	a\|1	CMP1\|10.0
2	sp	CMP1	1.0	2680.813757	0.0	a\|1	CMP1\|1.0
3	sp	CMP1	0.1	795.581544	0.0	a\|1	CMP1\|0.1
4	sp	CMP2	10.0	27495.080920	0.0	a\|1	CMP2\|10.0
...	...	...	...	...	...	...	...
651	ln	CMP1	0.1	288.879208	200.0	e\|9	CMP1\|0.1
652	ln	CMP2	10.0	3699.793797	200.0	e\|9	CMP2\|10.0
653	ln	CMP2	1.0	1954.964786	200.0	e\|9	CMP2\|1.0
654	ln	CMP2	0.1	282.778393	200.0	e\|9	CMP2\|0.1
655	ln	PC	NaN	16477.947980	200.0	e\|9	PC\|nan

	sample	stim	conc_µg_ml-1	ifng_pg_ml-1	dose_µg	animal_id	stim_id
32	sp	medium	0.0	21.736263	0.0	a\|5	medium\|0.0
33	sp	CMP1	10.0	148.652715	0.0	a\|5	CMP1\|10.0
34	sp	CMP1	1.0	126.871397	0.0	a\|5	CMP1\|1.0
35	sp	CMP1	0.1	113.452095	0.0	a\|5	CMP1\|0.1
36	sp	CMP2	10.0	582.270884	0.0	a\|5	CMP2\|10.0
...	...	...	...	...	...	...	...
99	ln	CMP1	0.1	NaN	0.0	a\|5	CMP1\|0.1
100	ln	CMP2	10.0	NaN	0.0	a\|5	CMP2\|10.0
101	ln	CMP2	1.0	NaN	0.0	a\|5	CMP2\|1.0
102	ln	CMP2	0.1	NaN	0.0	a\|5	CMP2\|0.1
103	ln	PC	NaN	NaN	0.0	a\|5	PC\|nan

		count	mean	std	min	25%	50%	75%	max
group	sex
a	f	1.0	85.000000	NaN	85.0	85.00	85.0	85.00	85.0
a	m	7.0	71.428571	15.714719	59.0	61.00	63.0	81.00	94.0
b	f	2.0	85.000000	0.000000	85.0	85.00	85.0	85.00	85.0
b	m	6.0	72.666667	16.524729	62.0	62.00	62.0	86.00	94.0
c	f	2.0	85.000000	0.000000	85.0	85.00	85.0	85.00	85.0
c	m	6.0	71.500000	12.243366	62.0	63.00	68.0	73.75	94.0
d	f	2.0	60.500000	2.121320	59.0	59.75	60.5	61.25	62.0
d	m	6.0	72.833333	2.401388	68.0	73.25	74.0	74.00	74.0
e	f	2.0	65.000000	4.242641	62.0	63.50	65.0	66.50	68.0
e	m	7.0	69.571429	5.593363	62.0	65.00	73.0	74.00	74.0

	Source	SS	DF	MS	F	p-unc	np2
0	sex	90.280589	1.0	90.280589	0.729665	0.399545	0.022996
1	group	332.166006	4.0	83.041501	0.671157	0.616959	0.079699
2	sex * group	833.096395	4.0	208.274099	1.683310	0.178957	0.178443
3	Residual	3835.595238	31.0	123.728879	NaN	NaN	NaN

	Effect	DFn	DFd	F	p	p<.05	ges
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<dbl>
1	dose_µg	4.00	35.00	7.345	2.10e-04	*	0.374
2	sample	1.00	35.00	13.879	6.85e-04	*	0.054
3	stim_id	2.48	86.74	399.037	1.00e-47	*	0.523
4	dose_µg:sample	4.00	35.00	0.160	9.57e-01		0.003
5	dose_µg:stim_id	9.91	86.74	6.417	2.73e-07	*	0.066
6	sample:stim_id	3.25	113.76	3.634	1.30e-02	*	0.005
7	dose_µg:sample:stim_id	13.00	113.76	2.609	3.00e-03	*	0.014

		Source	SS	DF1	DF2	MS	F	p-unc	ng2	eps
sample
ln	0	dose_µg	24.135926	4	35	6.033981	4.951279	2.861356e-03	0.324774	NaN
	1	stim_id	42.580516	5	175	8.516103	198.006717	4.519549e-70	0.459037	0.398419
	2	Interaction	4.533131	20	175	0.226657	5.269959	2.972874e-10	0.082853	NaN
sp	0	dose_µg	24.234661	4	35	6.058665	8.042147	1.042409e-04	0.443506	NaN
	1	stim_id	45.646615	5	175	9.129323	395.360919	2.452373e-93	0.600177	0.486179
	2	Interaction	2.294174	20	175	0.114709	4.967656	1.448528e-09	0.070152	NaN

group	dose (µg)
A	untreated
B	1.6
C	8
D	40
E	200

				count	mean	std	min	50%	max
dose_µg	sample	stim	conc_µg_ml-1
0.0	ln	CMP1	0.1	7.0	817.140381	796.610850	75.714679	577.098512	2400.309755
			1.0	7.0	2210.753717	2540.936318	115.467230	1872.883475	7678.569123
			10.0	7.0	3771.122427	3054.144372	304.119215	3519.635227	9707.649227
		medium	0.0	7.0	34.623566	18.816274	11.949131	28.705061	72.093906
	sp	CMP1	0.1	7.0	960.868304	223.934137	795.581544	917.903459	1456.012620
...	...	...	...	...	...	...	...	...	...
200.0	ln	medium	0.0	9.0	72.109845	55.625232	12.998240	57.308153	173.686927
	sp	CMP1	0.1	9.0	362.844449	304.948962	50.678500	294.487313	910.011775
			1.0	9.0	1073.192492	591.718820	170.867993	898.653612	2132.258257
			10.0	9.0	1948.155756	1257.182240	554.271580	1753.861420	3902.806026
		medium	0.0	9.0	49.876183	23.371521	22.478801	45.144443	103.555885

				count	mean	std	min	50%	max
dose_µg	sample	stim	conc_µg_ml-1
0.0	ln	CMP2	0.1	7.0	4105.579009	4129.743499	218.355170	3108.493300	12614.761430
			1.0	7.0	17519.735045	25335.639807	1111.065441	8779.501359	72953.894280
			10.0	7.0	37922.774760	45267.067137	1666.928359	21093.249890	135413.978600
		medium	0.0	7.0	34.623566	18.816274	11.949131	28.705061	72.093906
	sp	CMP2	0.1	7.0	5095.828155	3617.860716	2086.342701	3041.426623	10341.507490
...	...	...	...	...	...	...	...	...	...
200.0	ln	medium	0.0	9.0	72.109845	55.625232	12.998240	57.308153	173.686927
	sp	CMP2	0.1	9.0	718.298342	554.521045	128.962498	483.514389	1814.468536
			1.0	9.0	2699.733319	1800.224535	901.143809	2218.807523	5999.015402
			10.0	9.0	8944.243014	8225.876811	1701.793320	6359.198704	23246.583770
		medium	0.0	9.0	49.876183	23.371521	22.478801	45.144443	103.555885

			W	pval	equal_var
sample	stim_id
ln	CMP1\|0.1	levene	0.642492	0.635804	True
	CMP1\|1.0	levene	0.697553	0.598806	True
	CMP1\|10.0	levene	0.679567	0.610761	True
	CMP2\|0.1	levene	3.054415	0.029305	False
	CMP2\|1.0	levene	1.995815	0.116646	True
	CMP2\|10.0	levene	1.577794	0.201878	True
sp	CMP1\|0.1	levene	3.933091	0.009698	False
	CMP1\|1.0	levene	6.944668	0.000317	False
	CMP1\|10.0	levene	3.387082	0.019177	False
	CMP2\|0.1	levene	13.245935	0.000001	False
	CMP2\|1.0	levene	6.974131	0.000307	False
	CMP2\|10.0	levene	3.756314	0.012069	False