<jupyter_start><jupyter_text>Comparing Permutation and Bootstrap In this notebook we compare the permutation and bootstrap JER control methods in the two sample case.<jupyter_code>import numpy as np
import numpy.matlib as npm
import sys
import matplotlib.pyplot as plt 

import sanssouci as sa

import pyrft as pr
import random<jupyter_output><empty_output><jupyter_text>Initialize the example (using the equi-corr example here)<jupyter_code>random.seed(a=2021)

m = 2000                    # number of variables     (eg : number of voxels in an image)
n = 80                      # number of observations  (eg : number of images)
pi0 = 0.99                  # proportion of noise (true null hypotheses)
p0 = int(np.round(pi0 * m)) # number of true null hypotheses (non active voxels)
rho = 0.3                   # equi-correlation
s = 1                       # signal to noise ratio (because sigma is 1 CHECK!)

# generate noise
Z = np.random.randn(n, m)
w = np.random.randn(n)
W = Z.copy()
for i in range(m):
  W[:, i] = w # W has the same w for each subject to ensure the dependence
Y = np.sqrt(1 - rho) * Z + np.sqrt(rho) * W

# add signal (for false null hypotheses)
categ = np.random.binomial(1, 0.5, size=n) # Equal probability of each - but ususally not balanced in fact!
data = Y.copy()
w1 = np.where(categ == 1)
# Add signal to the data when the category is 1! :)
data[w1, (p0 + 1):m] = Y[w1, (p0 + 1): m] + s

# Record where the signal is as a boolean vector
bool_signal = np.ones((m, 1)) == 0
bool_signal[(p0 + 1): m] = True
sum(bool_signal)<jupyter_output><empty_output><jupyter_text>JER control using permutation<jupyter_code># Initialize the number of bootstraps and the alpha level
B = 100
alpha = 0.1
pval0 = sa.get_perm_p(data, categ, B=B , row_test_fun=sa.row_welch_tests)
piv_stat = sa.get_pivotal_stats(pval0, K=m)
lambda_perm = np.quantile(piv_stat, alpha)

print('Lambda Perm:', lambda_perm)
thr_perm = sa.t_linear(lambda_perm, np.arange(1, m + 1), m)<jupyter_output>Lambda Perm: 0.19195942582772574<jupyter_text>JER control using bootstrap<jupyter_code># Generate the design matrix from the category vector
X = pr.group_design(categ); 

# Set the contrast matrix
C = np.array([1, -1]); 

# Run the bootstrap JER algorithm
minPperm, orig_pvalues, pivotal_stats, bs = pr.boot_contrasts(data.transpose(), X, C, B)

# Obtain the lambda calibration
lambda_boot = np.quantile(pivotal_stats, alpha)
print('Lambda Boot:', lambda_boot)

thr_boot = sa.t_linear(lambda_boot, np.arange(1, m + 1), m)
import matplotlib.pyplot as plt 
plt.hist(minPperm)<jupyter_output><empty_output><jupyter_text>Calculate the p-values (note that the bootstrap and permutation methods as coded use different two sample test so the p-values are slightly different)<jupyter_code>rwt = sa.row_welch_tests(data, categ)
pvals_perm = rwt['p_value']
print(pvals_perm.shape)
pval_sort_idx_perm = np.argsort(pvals_perm)

pval_sort_idx_boot = np.argsort(np.ravel(orig_pvalues.field))
pvals_boot = np.ravel(orig_pvalues.field)[pval_sort_idx_boot]<jupyter_output>(2000,)<jupyter_text>Compare permutation and bootstrap p-values The p values are very close between the two methods (with use different 2 sample tests) but are slightly different, especially in the tails<jupyter_code>print(pvals_perm[0: 10])
print(pvals_boot[0: 10])
print(pvals_perm[1000: 1010])
print(pvals_boot[1000: 1010])

plt.plot(pvals_perm[pval_sort_idx_boot], pvals_boot)
plt.xlabel('Permutation p-values')
plt.ylabel('Bootstrap p-values')
plt.title('Comparing the p-values from the two methods')<jupyter_output>[0.03247533 0.97933297 0.59796329 0.94204189 0.2767466  0.38967197
 0.84673441 0.01325814 0.05675138 0.28444079]
[5.61548585e-11 2.86722202e-10 3.66596531e-10 5.21172883e-10
 1.71791337e-09 1.04356601e-08 1.51298616e-08 2.25915850e-08
 3.08375017e-08 3.89408796e-08]
[0.61238534 0.00709827 0.05340741 0.0684976  0.54325549 0.28449149
 0.80367945 0.03700122 0.01571187 0.10648225]
[0.25240828 0.25315707 0.25386105 0.25489101 0.25569662 0.25624043
 0.25695789 0.25760985 0.25772658 0.25785108]<jupyter_text>Generate a histogram of the p-values Note that for low numbers of voxels this is not uniform because the data is highly correlated<jupyter_code>plt.hist(pvals_perm)<jupyter_output><empty_output><jupyter_text>Plot the FPR and TP curve bounds<jupyter_code>max_FP_perm = sa.curve_max_fp(pvals_perm, thr_perm)
max_FP_boot = sa.curve_max_fp(pvals_boot, thr_boot)

# Generate the vector [0,...,m]
one2npvals = np.arange(1, m + 1)

# Choose the number of p-values (always the smallest ones first) to plot
lowestnumber = 100

# Ensure that selected number is not greater than the total number of p-values
lowestnumber = np.min([lowestnumber, m])

# Dividing the envelope by the number of elements in the set gives a bound on the false discovery proportion
max_FDP_boot = max_FP_boot[0: lowestnumber] / one2npvals[0: lowestnumber] 
min_TP_boot  = one2npvals[0: lowestnumber] - max_FP_boot[0: lowestnumber]
max_FDP_perm = max_FP_perm[0: lowestnumber] / one2npvals[0: lowestnumber] 
min_TP_perm  = one2npvals[0: lowestnumber] - max_FP_perm[0: lowestnumber]

# Calculate the truth (to determine if it is correctly bounded!)
sorted_signal = bool_signal[pval_sort_idx_boot]
TP = np.zeros(lowestnumber)
for I in np.arange(lowestnumber):
    TP[I] = np.sum(sorted_signal[0: I + 1])
    
# Calculate the true FDP for each subset
FP = np.zeros(lowestnumber)
for I in np.arange(lowestnumber):
    FP[I] = np.sum(abs(sorted_signal[0: I + 1] - 1))
true_FDP = FP / one2npvals[0: lowestnumber] 
    
# Initialize the figure
figure = plt.figure(figsize=(10, 4))

# Plot the false discovery proportion and its bound
plt.subplot(121)
plt.plot(max_FDP_boot, label = 'Bootstrap FDP bound')
plt.plot(max_FDP_perm, label = 'Permutation FDP bound')
plt.plot(true_FDP, label = 'True FDP')
plt.title('Upper bound on FDP amongst smallest p-values')
plt.xlim(1, lowestnumber)
plt.xlabel('k')
plt.ylabel('FDP(p_{(1)}, \dots, p_{(k)}')
plt.legend(loc="lower right")

# Plot the true postives and their bound
plt.subplot(122)
plt.plot(min_TP_perm, label = 'Bootstrap TP bound')
plt.plot(min_TP_perm, label = 'Permutation TP bound')
plt.plot(TP, label = 'True Positives')
plt.title('Lower bound on TP')
plt.legend(loc="lower right")
plt.xlim(1, lowestnumber)
plt.xlabel('k')
plt.ylabel('TP(p_{(1)}, \dots, p_{(k)}')
#figure, axes = plt.subplots(nrows=1, ncols=2) 
figure.tight_layout(pad=1.0)<jupyter_output>The input p-values were not sorted -> this is done now<jupyter_text>Analyzing the impact of lambda changes on the FDP curve Suprisingly the value of lambda doesn't appear to make much of a difference!<jupyter_code># Initialize the number of bootstraps and the alpha level
B = 100
alpha = 0.1
pval0 = sa.get_perm_p(data, categ, B=B , row_test_fun=sa.row_welch_tests)
piv_stat = sa.get_pivotal_stats(pval0, K=m)
lambda_perm = np.quantile(piv_stat, alpha)

print('Lambda Perm:', lambda_perm)
thr_perm = sa.t_linear(lambda_perm, np.arange(1, m + 1), m)

n2use = 100
lambda_set = np.arange(10) / 10
print(lambda_set)
lambda_set = lambda_set[1:]
for I in np.arange(len(lambda_set)):
    thresh = sa.t_linear(lambda_set[I], np.arange(1, m + 1), m)
    max_FP_perm = sa.curve_max_fp(pvals_perm, thresh)
    max_FDP_perm = max_FP_perm[0:n2use]/one2npvals[0: n2use] 
    plt.plot(max_FDP_perm)
plt.plot(true_FDP[0: n2use], label = 'True FDP')
plt.legend(loc="lower right")
lambda_set = np.arange(10) / 10
lambda_set = lambda_set[1:]
print(lambda_set)
n2use = 100
for I in np.arange(len(lambda_set)):
    thresh = sa.t_linear(lambda_set[I], np.arange(1, m + 1), m)
    max_FP_perm = sa.curve_max_fp(pvals_perm, thresh)
    min_TP_perm = np.arange(1, m + 1) - max_FP_perm
    plt.plot(min_TP_perm[0: n2use])
plt.plot(TP[0: n2use], label = 'True FDP')
plt.legend(loc="lower right")<jupyter_output>[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now<jupyter_text>Analyzing for very low lambda<jupyter_code>lambda_set = np.arange(10) / 10000
for I in np.arange(len(lambda_set)):
    thresh = sa.t_linear(lambda_set[I], np.arange(1, m + 1), m)
    max_FP_perm = sa.curve_max_fp(pvals_perm, thresh)
    plt.plot(max_FDP_perm)
plt.plot(true_FDP, label = 'True FDP')<jupyter_output>The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now
The input p-values were not sorted -> this is done now<jupyter_text>One Sample Example<jupyter_code>random.seed(a=2021)

m = 2000                    # number of variables     (eg : number of voxels in an image)
n = 80                      # number of observations  (eg : number of images)
pi0 = 0.99                  # proportion of noise (true null hypotheses)
p0 = int(np.round(pi0 * m)) # number of true null hypotheses (non active voxels)
rho = 0.3                   # equi-correlation

# generate noise
Z = np.random.randn(n, m)
w = np.random.randn(n)
W = Z.copy()
for i in range(m):
  W[:, i] = w # W has the same w for each subject to ensure the dependence
Y = np.sqrt(1 - rho) * Z + np.sqrt(rho) * W

# add signal (for false null hypotheses)
categ = np.zeros(n) # Equal probability of each - but ususally not balanced in fact!
data = Y.copy()<jupyter_output><empty_output><jupyter_text>Running bootstrap procedure<jupyter_code># Generate the design matrix from the category vector
X = pr.group_design(categ)

# Set the contrast matrix
C = np.array(1)

B = 100; alpha = 0.1

# Run the bootstrap JER algorithm
minPperm, orig_pvalues, pivotal_stats, _ = pr.boot_contrasts(data.transpose(), X, C, n_bootstraps = B)

# Obtain the lambda calibration
lambda_boot = np.quantile(pivotal_stats, alpha)
print('Lambda Boot:', lambda_boot)

thr_boot = sa.t_linear(lambda_boot, np.arange(1, m + 1), m)
print(len(pivotal_stats))
# Running permutation procedure
# Initialize the number of bootstraps and the alpha level
B = 100
alpha = 0.1
pval0 = sa.get_permuted_p_values_one_sample(data, B=1000)
piv_stat = sa.get_pivotal_stats(pval0)
lambda_perm = np.quantile(piv_stat, alpha)

print('Lambda Perm:', lambda_perm)
thr_perm = sa.t_linear(lambda_perm, np.arange(1, m + 1), m)
thr_simes = sa.t_linear(alpha, np.arange(1, m + 1), m)
# sa.min_tp(cluster_p_values, thr)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11068943706514864, 'fraction_numerical': 0.05683563748079877, 'mean_word_length': 4.587077233720343, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '378373', 'license': 'MIT', 'n_tokens_mistral': 4510, 'n_tokens_neox': 4026, 'n_words': 1193, 'path': 'examples/.ipynb_checkpoints/compare_perm_and_boot-checkpoint.ipynb', 'repo_name': 'alexblnn/pyrft', 'size': '11067'}