Statistical-Scripts/anova.py at main · sina-programer/Statistical-Scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
from functools import partial
from scipy import stats
import itertools as it
import pandas as pd
import numpy as np


def sum_of_between_errors(series, delta=0):
    # return series.count() * np.square(series.mean() - delta)
    return np.count_nonzero(~np.isnan(series)) * np.square(np.nanmean(series) - delta)


def anova(df: pd.DataFrame, target, *factors):
    factors = list(factors)
    gdf = df.groupby(factors)[target]
    gdfs = {f: df.groupby(f)[target] for f in factors}

    mu = df[target].mean()
    mdfs = {f: g.mean() for f, g in gdfs.items()}
    se = partial(sum_of_between_errors, delta=mu)

    n = len(df)
    k = len(gdf)
    ks = {f: len(g) for f, g in gdfs.items()}
    output = {}

    total = {}
    total['df'] = n - 1
    total['ss'] = float(np.square(df[target] - mu).sum())
    total['ms'] = total['ss'] / total['df']
    total['f'] = np.nan
    total['p'] = np.nan
    output['total'] = total

    error = {}
    error['df'] = n - k
    error['ss'] = float(gdf.transform(lambda x: np.square(x - x.mean())).sum())
    error['ms'] = error['ss'] / error['df']
    error['f'] = np.nan
    error['p'] = np.nan
    output['error'] = error

    model = {}
    model['df'] = k - 1
    model['ss'] = float(gdf.apply(se).sum())
    model['ms'] = model['ss'] / model['df']
    model['f'] = model['ms'] / error['ms']
    model['p'] = float(stats.f.sf(model['f'], model['df'], error['df']))
    output['model'] = model

    for i, factor in enumerate(factors):
        info = {}
        info['df'] = ks[factor] - 1
        info['ss'] = float(gdfs[factor].apply(se).sum())
        info['ms'] = info['ss'] / info['df']
        info['f'] = info['ms'] / error['ms']
        info['p'] = float(stats.f.sf(info['f'], info['df'], error['df']))
        output[(i+1,)] = info

    rng = list(range(1, len(factors)+1))
    for r in rng[1:]:
        for comb in it.combinations(rng, r):
            fs = np.take(factors, np.subtract(comb, 1)).tolist()
            gf = df.groupby(fs)[target]
            mf = gf.mean()

            def interaction(sdf):
                return sdf.count() * np.square(
                    mf[sdf.name]
                    - sum([mdfs[fs[i]][sdf.name[i]] for i in range(len(fs))])
                    + mu * (len(fs)-1)
                ).sum()

            info = {}
            info['df'] = int(np.prod([output[(i,)]['df'] for i in comb]))
            info['ss'] = float(gf.apply(interaction).sum())
            info['ms'] = info['ss'] / info['df']
            info['f'] = info['ms'] / error['ms']
            info['p'] = float(stats.f.sf(info['f'], info['df'], error['df']))
            output[comb] = info

    return output