mixed_level_model.py

# Following https://m-clark.github.io/docs/mixedModels/anovamixed.html

import numpy as np
from scipy import stats
from pandas import DataFrame
import statsmodels.api as sm
from statsmodels.stats import anova
import statsmodels.formula.api as smf
import pingouin

treat = ['treat'] * 5 + ['control'] * 5
pre = np.array([20, 10, 60, 20, 10, 50, 10, 40, 20, 10])
post = np.array([70, 50, 90, 60, 50, 20, 10, 30, 50, 10])
change = post - pre

df = DataFrame(dict(id=np.arange(1, len(pre) + 1), treat=treat, pre=pre, post=post))
print(df.head())

dflong = \
    df.melt(
        ['id', 'treat'], var_name='time', value_name='score'
        ).sort_values(
            ['id', 'time'], ascending=[True, False]
            ).reset_index(drop=True)
print(dflong.head())

###############################################################################
# Post only
print('\nPost only')

#
# ttest
#
post_treat = df.query('treat == "treat"')['post']
post_control = df.query('treat == "control"')['post']
t, p = stats.ttest_ind(post_treat, post_control)
print(f'ttest: t²={t ** 2}\n       p={p}')

#
# ANOVA
#
f, p = stats.f_oneway(post_treat, post_control)
print(f'ANOVA: f={f}\n       p={p}')
lm = smf.ols('post ~ treat', data=df).fit()
table = sm.stats.anova_lm(lm, typ=2)
print(f'ANOVA: F={table.F[0]}\n       p={table["PR(>F)"][0]}')
table = pingouin.anova(df, 'post', between='treat')
print(f'ANOVA: F={table.F[0]}\n       p={table["p-unc"][0]}')

###############################################################################
# Pre-post

print('\nPre-post')
is_treat = np.array(treat) == 'treat'
change_treat = post[is_treat] - pre[is_treat]
change_control = post[~is_treat] - pre[~is_treat]
t, p = stats.ttest_ind(change_treat, change_control)
print(f'ttest: t²={t ** 2}\n       p={p}')

lm = smf.ols('post ~ pre + treat', data=df).fit()
table = sm.stats.anova_lm(lm, typ=2)
print(f'ANOVA: F={table.F[0]}\n       p={table["PR(>F)"][0]}')
table = pingouin.ancova(df, 'post', between='treat', covar='pre')
print(f'ANOVA: F={table.F[0]}\n       p={table["p-unc"][0]}')

# "a t-test on the change score is the interaction term from the repeated
# measures ANOVA on the same data"
table = pingouin.mixed_anova(dflong, 'score', within='time', subject='id',
                             between='treat', correction=False)
print(f'ANOVA: F={table.F[2]}\n       p={table["p-unc"][2]}')
lm = smf.ols('score ~ treat*time + (1|id)', dflong).fit()
table = sm.stats.anova_lm(lm, typ=2)
print(f'ANOVA: F={table.F[2]}\n       p={table["PR(>F)"][2]}')

md = smf.mixedlm('score ~ treat*time', dflong, groups=dflong["id"]).fit()
print(f'MIXED: F={table.F[2]}\n       p={table["PR(>F)"][2]}')