Non Parameteric Tests - Permutation Test

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
from scipy.stats import permutation_test

treat = np.array([114, 203, 217, 254, 256, 284, 296])
control  = np.array([4, 7, 24, 25, 48, 71, 294])

def statistic(x, y, axis):
    return np.median(x, axis) - np.median(y, axis)

res = permutation_test((treat, control), statistic, vectorized=True,
                       n_resamples=np.inf, alternative='greater') ##n_samples=np.infはexact test

res_median, res_pvalue, res_null = res.statistic, res.pvalue, res.null_distribution
print(res_pvalue)
>>> 0.017482517482517484

1
2
np.count_nonzero(res_null >= res_median - max(1e-14, res_median * 1e-14)) / len(res_null)
>>> 0.017482517482517484

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def two_sample_bootstrap(x, y, statistic, size, percentile):
    boot_res = []
    for i in range(size):
        tmp_x = np.random.choice(x, len(x), replace=True) ## x[np.random.randint(0, len(x), len(x))]の方が早い
        tmp_y = np.random.choice(y, len(y), replace=True)
        tmp_stats = statistic(tmp_x, tmp_y)
        boot_res.append(tmp_stats)
    
    return np.quantile(boot_res, [1 - percentile, percentile])

def statistic_median(x, y):
    return np.median(x) - np.median(y)

two_sample_bootstrap(treat, control, statistic_median, 100000, 0.975)
>>> (132, 260) ### Rのsimpleboot, bootと結果は一致

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
from scipy.stats import permutation_test

def statistic(x, y):
    return np.median(x) - np.median(y)

x = [3.67, 1.72, 3.46, 2.60, 2.03, 2.10, 3.01]
y = [2.11, 1.79, 2.71, 1.89, 1.69, 1.71, 2.01]

res_sci = permutation_test((x, y), statistic, vectorized=False,
                       n_resamples=np.inf,
                       permutation_type='samples',
                       alternative='two-sided')
print(res_sci.pvalue)
>>> 0.0625

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
x = [3.67, 1.72, 3.46, 2.60, 2.03, 2.10, 3.01]
y = [2.11, 1.79, 2.71, 1.89, 1.69, 1.71, 2.01]
z = np.array([x, y])

def two_sample_permutation(z, statistic, size):
    boot_res = []
    z_len = z.shape[1]
    array_idx = np.arange(0, z_len)
    for i in range(size):
        flip_coin = np.random.randint(0, 2, z_len)
        tmp_x = z[flip_coin, array_idx]
        tmp_y = z[1 - flip_coin, array_idx]
        tmp_stats = statistic(tmp_x, tmp_y)
        boot_res.append(tmp_stats)
    
    return boot_res

res = two_sample_permutation(z=z, statistic=statistic, size=1000)
obs = np.median(y) - np.median(x)
np.count_nonzero(np.abs(res) >= np.abs(obs)) / len(res)
>>> 0.062

1
2
3
4
5
6
7
8
df = pd.read_stata('https://github.com/QuantEcon/lecture-python/blob/master/source/_static/lecture_specific/ols/maketable1.dta?raw=true')
df = df.dropna()
y, x = df['logpgp95'].values, df['avexpr'].values

plt.scatter(x, y)
plt.title('size = {}, correlation = {:.2}'.format(len(x), pearsonr(x, y)[0]))
plt.xlabel('avexpr')
plt.ylabel('logpgp95');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
### Permutation Test
def statistic_r(x, y):
    return pearsonr(x, y)[0]

res = permutation_test((x, y), statistic_r, vectorized=False,
                       permutation_type='pairings',
                       alternative='greater')

### 可視化
count, bins, _ = plt.hist(res.null_distribution, bins=100)
plt.vlines(obs_corr, 0, max(count), 
           colors='red', linestyles='dashed'
          )
plt.title('Permuation Test: Corr coef Distirbution: p-value = {}'.format(res.pvalue))

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
import statsmodels.api as sm
import itertools
from multiprocessing import Pool
import collections

N = [15, 30]
const = [1]
beta = np.linspace(0, 1, 11)
sig_level = 0.05
nrep = 2000
seed = 42

def flatten(l):
    for el in l:
        if isinstance(el, collections.abc.Iterable) and not isinstance(el, (str, bytes)):
            yield from flatten(el)
        else:
            yield el

def func_dgp(size, const, beta):
    eps = np.random.normal(0, 1, size)
    x = np.random.normal(0, 1, size)
    y = const + x * beta + eps

    return y, x

def compute_ttest(y, X):
    _X = sm.add_constant(X)
    model = sm.OLS(y,_X)
    result = (model.fit()).t_test([0,1])
    return result.effect, result.pvalue

def compute_permtest(y, X):
    def perm_statistics(_y, _X):
        return np.cov(_y, _X, ddof=-1)[0,1]/np.var(_X, ddof=-1)

    res = permutation_test((y, X), perm_statistics, vectorized=False,
                            n_resamples=2000,
                            permutation_type='pairings',
                            alternative='two-sided')
    
    return res.statistic, res.pvalue


def compute_power_analysis(params):
    record = []
    t_power, p_power = [], []
    for i in range(nrep):
        y, X = func_dgp(*params)
        t_coef, t_pval = compute_ttest(y, X)
        t_power.append(t_pval)
        p_coef, p_pval = compute_permtest(y, X)
        p_power.append(p_pval)

    power_list = [np.mean(np.array(t_power) < sig_level), np.mean(np.array(p_power) < sig_level)]

    params = list(params)
    params.extend(*zip([power_list]))

    record.extend(list(flatten(params)))
    
    return record

## table 
p = Pool(6)
res = p.map(compute_power_analysis, [param for param in list(itertools.product(*[N, const, beta]))])
table_1 = pd.DataFrame(res, columns=['obs_size', 'const', 'beta', 't_test_power', 'perm_test_power'])

## p-value distribution
p = Pool(6)
res_pval_dist =  p.map(compute_power_analysis_pval, [(15, 1, 0, 2000)])
plt.hist(res_pval_dist, bins=np.linspace(0, 1, 11), density=True);

from scipy import stats
stats.kstest(res_pval_dist,stats.norm.cdf)
>>> KstestResult(statistic=0.8413447460685429, pvalue=0.31731050786291415)