<jupyter_start><jupyter_text>Load both well-level data and synthesis summary data<jupyter_code>df = pd.read_csv('synthetic_data.csv')
df.head()
synthesis_df = pd.read_csv('synthetic_synthesis_data.csv')
synthesis_df.head()
# Get baseline synthesis mean and standard deviaiton for pure yield
'Pure Yield: {} +/- {} nmol'.format(synthesis_df.pure_yield.mean(), synthesis_df.pure_yield.std())
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(synthesis_df.pure_yield)
ax.set_xlabel('Pure Yield', fontsize=14)
ax.set_ylabel('Density', fontsize=14);<jupyter_output>/Users/aaron.wiegel/anaconda3/envs/py37/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval<jupyter_text>Simple power simulationWe want to test whether a process change will improve our pure yield. Luckily, we know our baseline "population" pure yield, assuming our pure yields do not drift over time. We then declare the null and alternative hypotheses for this test:$H_0: \mu = 29.95$ nmol$H_a: \mu > 29.95$ nmolWe can set the significance level for whatever resulting hypothesis test we would like to perform at the customary level of 5%. This determines our false positive rate, i.e., how often we would reject the null hypothesis when nothing has actually changed. How do we avoid false negatives, though? (When we don't reject the null hypothesis when we should.) In other words, how many times should we run experiments before we can be reasonably sure that we haven't actually improved pure yield?For that we can use power simulations. (Sometimes it's possible to calculate this directly, but I think simulations better illustrate the concepts.) We need to figure out an acceptable false negative rate (customarily set to 20%). Importantly, we also need to assume something about the _expected_ effect size. To start, we will assume a 10% effect size but relax this assumption later.<jupyter_code>effect_size = 1.10
false_negative_rate = 0.2
population_pure_yield = synthesis_df['pure_yield'].mean()
n_experiments = 3
n_experiments_max = 31<jupyter_output><empty_output><jupyter_text>Let's start by considering an experiment with 3 replicates (and recall that a synthesis, not a well, is a replicate.) We want to find pure yield under the null hypothesis, and the pure yield under the alternative hypothesis (assuming an effect size of 10%)<jupyter_code># Set random number seed to ensure reproducibility
np.random.seed(42)
results = []
for i in range(n_experiments, n_experiments_max):
    for _ in range(1000):
        sample_null = synthesis_df.sample(n=i)
        sample_alt = synthesis_df.sample(n=i)

        sample_alt['pure_yield'] = sample_alt['pure_yield'] * effect_size

        t_null = scipy.stats.ttest_1samp(sample_null['pure_yield'], population_pure_yield)[0]
        t_alt = scipy.stats.ttest_1samp(sample_alt['pure_yield'], population_pure_yield)[0]

        data = {
            'n_experiments': i,
            'pure_yield_null': sample_null['pure_yield'].mean(),
            'pure_yield_null_std': sample_null['pure_yield'].std(),
            't_stat_null': t_null,
            # Convert two-sided results to one-sided
            'p_value_null': 1 - scipy.stats.t.cdf(t_null, i - 1),
            'pure_yield_alt': sample_alt['pure_yield'].mean(),
            'pure_yield_alt_std': sample_alt['pure_yield'].std(),
            't_stat_alt': t_alt,
            # Convert two-sided results to one-sided
            'p_value_alt': 1 - scipy.stats.t.cdf(t_alt, i - 1),
        }
        results.append(data)
    
sample_df = pd.DataFrame(results)
sample_df.head()
mask = sample_df.n_experiments == 3
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(sample_df[mask].pure_yield_null, hist=False, ax=ax, label='Null')
sns.distplot(sample_df[mask].pure_yield_alt, hist=False, ax=ax, label='Alternative')
ax.set_xlabel('Pure Yield', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
ax.set_xlim([20, 45])
fig.savefig('null_alt_comp_n_3.png', dpi=150, bbox_inches='tight');
sample_df.describe()
sample_df['reject_null_alt'] = sample_df['p_value_alt'] < 0.05
sample_df['reject_null_null'] = sample_df['p_value_null'] < 0.05
sample_df['reject_null_alt'].mean(), sample_df['reject_null_null'].mean()
sample_df.loc[sample_df.n_experiments == 3, 'reject_null_alt'].mean()
results = []
for n_experiments, group in sample_df.groupby('n_experiments'):
    data = {
        'n_experiments': n_experiments,
        'reject_null_alt': group['reject_null_alt'].mean(),
        'reject_null_null': group['reject_null_null'].mean()
    }
    results.append(data)
power_df = pd.DataFrame(results)
power_df['false_negative_rate'] = 1 - power_df['reject_null_alt']
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(power_df['n_experiments'], power_df['false_negative_rate'], label='False negative')
ax.plot(power_df['n_experiments'], power_df['reject_null_null'], label='False positive')
ax.plot([3, 30], [0.2, 0.2], '--')
ax.legend()
ax.set_xlabel('Number of experiments', fontsize=14)
ax.set_ylabel('Frequency', fontsize=14)
fig.savefig('type1_type2_comparison.png', dpi=150, bbox_inches='tight');
experiment_mask = sample_df.n_experiments == 22
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(sample_df[experiment_mask].pure_yield_null, hist=False, ax=ax, label='Null')
sns.distplot(sample_df[experiment_mask].pure_yield_alt, hist=False, ax=ax, label='Alternative')
ax.set_xlabel('Pure Yield', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
ax.set_xlim([20, 45])
fig.savefig('null_alt_comp_n_22.png', dpi=150, bbox_inches='tight');<jupyter_output><empty_output><jupyter_text>Issues with a particular side<jupyter_code>df['other_impurity_pass'] = df['other_impurity'] < 0.03
bottom_mask = df.well_position_label.str.contains('H')
df.loc[bottom_mask, 'other_impurity_pass'].mean()
df.loc[~bottom_mask, 'other_impurity_pass'].mean()
results = []
for synthesis_id, group in df.groupby('synthesis_id'):
    data = {
        'synthesis_id': synthesis_id,
        'bottom_pure_yield': group.loc[bottom_mask, 'pure_yield'].mean(),
        'top_pure_yield': group.loc[~bottom_mask, 'pure_yield'].mean(),
        'overall_pass_rate': group['other_impurity_pass'].mean(),
        'bottom_pass_rate': group.loc[bottom_mask, 'other_impurity_pass'].mean(),
        'top_pass_rate': group.loc[~bottom_mask, 'other_impurity_pass'].mean(),
    }
    results.append(data)
position_df = pd.DataFrame(results)
position_df['bottom_pass_rate'].hist()
position_df['top_pass_rate'].hist()
position_df['pass_rate_diff'] = position_df['bottom_pass_rate'] - position_df['top_pass_rate']
position_df['pass_rate_diff'].hist()
position_df['pass_rate_diff'].mean()
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(position_df['overall_pass_rate'], ax=ax)
ax.set_xlabel('Other impurity pass rate', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
fig.savefig('other_impurity_pass_rate.png', dpi=150, bbox_inches='tight')
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(position_df['top_pass_rate'], ax=ax, hist=False, label='Top 7 rows')
sns.distplot(position_df['bottom_pass_rate'], ax=ax, hist=False, label='Bottom row')
ax.set_xlabel('Other impurity pass rate', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
fig.savefig('other_impurity_pass_rate_comparison.png', dpi=150, bbox_inches='tight')
 fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(position_df['pass_rate_diff'], ax=ax, label='Bottom - Top', hist=False)
ax.set_xlabel('Other impurity pass rate difference', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
fig.savefig('other_impurity_pass_rate_difference.png', dpi=150, bbox_inches='tight')
results = []
for i in range(2, 5):
    for _ in range(1000):
        sample = position_df.sample(n=i)
        data = {
            'sample_size': i,
            'pass_rate_difference_mean': sample['pass_rate_diff'].mean()
        }
        results.append(data)
sample_df = pd.DataFrame(results)
fig, ax = plt.subplots(figsize=(5, 5))
sns.distplot(position_df['pass_rate_diff'], ax=ax, label='Synthesis Pass Rate', hist=False)
sns.distplot(sample_df.loc[sample_df['sample_size'] == 2,
                       'pass_rate_difference_mean'], ax=ax, label='2 Syntheses', hist=False)
sns.distplot(sample_df.loc[sample_df['sample_size'] == 4,
                           'pass_rate_difference_mean'], ax=ax, label='4 Syntheses', hist=False)
ax.set_xlabel('Pass Rate Difference', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
fig.savefig('sampling_pass_rate_differences.png', dpi=150, bbox_inches='tight')
fig, ax = plt.subplots(figsize=(6, 6))
fig.suptitle('Bottom-to-Top correlation on 96-well plates', fontsize=14)
ax.set_xlabel('Pure yield on bottom row', fontsize=12)
ax.set_ylabel('Pure yield on top 7 rows', fontsize=12)
ax.set_xlim([10, 50])
ax.set_ylim([10, 50])
plt.scatter(position_df['bottom_pure_yield'],
            position_df['top_pure_yield'], alpha=0.1)
fig.savefig('bottom_top_pure_yield.png', bbox_inches='tight', dpi=150);
fig, ax = plt.subplots(figsize=(6, 6))
fig.suptitle('Bottom-to-Top correlation on 96-well plates', fontsize=14)
ax.set_xlabel('Other impurity pass rate on bottom row', fontsize=12)
ax.set_ylabel('Other impurity pass rate on top 7 rows', fontsize=12)
plt.scatter(position_df['bottom_pass_rate'],
            position_df['top_pass_rate'], alpha=0.1)
fig.savefig('bottom_top_pass_rate.png', bbox_inches='tight', dpi=150);
results = []
population_pass_rate_diff = -0.234
synthesis_ids = df.synthesis_id.unique()
for i in range(2, 6):
    for _ in range(1000):
        random_ids = np.random.choice(synthesis_ids, replace=False, size=i)
        
        bottom_pass = []
        bottom_pass_fixed = []
        top_pass = []
        
        for id_ in random_ids:
            synthesis_mask = df.synthesis_id == id_
            fixed_sample = df[synthesis_mask & ~bottom_mask].sample(n=12, replace=True)
            
            bottom_pass.append(df.loc[synthesis_mask & bottom_mask, 'other_impurity_pass'].mean())
            bottom_pass_fixed.append(fixed_sample['other_impurity_pass'].mean())
            top_pass.append(df.loc[synthesis_mask & ~bottom_mask, 'other_impurity_pass'].mean())
        
        pass_rate_diff = np.array(bottom_pass) - np.array(top_pass)
        pass_rate_diff_fixed = np.array(bottom_pass_fixed) - np.array(top_pass)
        
        t_null = scipy.stats.ttest_1samp(pass_rate_diff, population_pass_rate_diff)[0]
        t_alt = scipy.stats.ttest_1samp(pass_rate_diff_fixed, population_pass_rate_diff)[0]
        
        data = {
            'n_experiments': i,
            'pass_rate_diff_null': np.mean(pass_rate_diff),
            't_stat_null': t_null,
            # Convert two-sided results to one-sided
            'p_value_null': 1 - scipy.stats.t.cdf(t_null, i - 1),
            'pass_rate_diff_alt': np.mean(pass_rate_diff_fixed),
            't_stat_alt': t_alt,
            # Convert two-sided results to one-sided
            'p_value_alt': 1 - scipy.stats.t.cdf(t_alt, i - 1),
        }
        
        results.append(data)

pass_rate_df = pd.DataFrame(results)
pass_rate_df['reject_null_alt'] = pass_rate_df['p_value_alt'] < 0.05
pass_rate_df['reject_null_null'] = pass_rate_df['p_value_null'] < 0.05
results = []
for n_experiments, group in pass_rate_df.groupby('n_experiments'):
    data = {
        'n_experiments': n_experiments,
        'false_positive_rate': group['reject_null_null'].mean(),
        'false_negative_rate': 1 - group['reject_null_alt'].mean()
    }
    results.append(data)
power_pass_df = pd.DataFrame(results)
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(power_pass_df['n_experiments'], power_pass_df['false_negative_rate'], label='False negative')
ax.plot(power_pass_df['n_experiments'], power_pass_df['false_positive_rate'], label='False positive')
ax.plot([2, 5], [0.2, 0.2], '--')
ax.legend()
ax.set_xlabel('Number of experiments', fontsize=14)
ax.set_ylabel('Frequency', fontsize=14)
fig.savefig('type1_type2_pass_rate_comparison.png', dpi=150, bbox_inches='tight');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13558909444985395, 'fraction_numerical': 0.019961051606621226, 'mean_word_length': 4.908437200383509, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 17, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '601989', 'license': 'MIT', 'n_tokens_mistral': 4523, 'n_tokens_neox': 4203, 'n_words': 855, 'path': 'PowerSimulations.ipynb', 'repo_name': 'aawiegel/PowerSimulationDemo', 'size': '12324'}