|
defoverlay_effect_size(ax):height = ax.get_ylim()[ 1] * 0.5ax.hlines(height, 0, 0.2, 'red', lw= 5) ax.hlines(height, 0.2, 0.8, 'blue', lw= 5) ax.hlines(height, 0.8, ax.get_xlim()[ 1], 'green', lw= 5)ax = pm.plot_posterior(kruschke_trace[ 2000:], varnames=[ 'effect_size'],color= '#87ceeb')[ 0]overlay_effect_size(ax)
Effect size (Cohen's d, none to small,medium, large) could be anywhere from essentially nothing to large (95% HPD [0.0, 0.77])。 IQ改善0-4 该药很可能无关紧要。 没有生物学意义的证据。 例子2：手机消毒问题两种常用的方法相比，我的“特别方法”能更好的消毒我的手机吗？ the experiment design 随机将手机分配到六组之一：4“特别”方法+ 2“对照”方法。 count 形成的细菌菌落数，比较前后的计数。 renamed_treatments = dict()renamed_treatments[ 'FBM_2'] = 'FM1'renamed_treatments[ 'bleachwipe'] = 'CTRL1' renamed_treatments[ 'ethanol'] = 'CTRL2' renamed_treatments[ 'kimwipe'] = 'FM2' renamed_treatments[ 'phonesoap'] = 'FM3' renamed_treatments[ 'quatricide'] = 'FM4' # Reload the data one more time. data = pd.read_csv( 'smartphone_sanitization_manu.csv', na_values=[ '#DIV/0!']) deldata[ 'perc_reduction colonies'] # Exclude cellblaster data data = data[data[ 'treatment'] != 'CB30']data = data[data[ 'treatment'] != 'cellblaster'] # Rename treatments data[ 'treatment'] = data[ 'treatment'].apply( lambdax: renamed_treatments[x]) # Sort the data according to the treatments. treatment_order = [ 'FM1', 'FM2', 'FM3', 'FM4', 'CTRL1', 'CTRL2']data[ 'treatment'] = data[ 'treatment'].astype( 'category')data[ 'treatment'].cat.set_categories(treatment_order, inplace= True)data = data.sort_values([ 'treatment']).reset_index(drop= True) # Encode the treatment index. data[ 'treatment_idx'] = data[ 'treatment'].apply( lambdax: treatment_order.index(x))data[ 'perc_change_colonies'] = (data[ 'colonies_post'] - data[ 'colonies_pre']) / data[ 'colonies_pre'] # # View the first 5 rows. # data.head(5) # # filter the data such that we have only PhoneSoap (PS-300) and Ethanol (ET) # data_filtered = data[(data['treatment'] == 'PS-300') | (data['treatment'] == 'QA')] # data_filtered = data_filtered[data_filtered['site'] == 'phone'] # data_filtered.sample(10)数据defplot_colonies_data():fig = plt.figure(figsize=( 10, 5)) ax1 = fig.add_subplot( 2, 1, 1) sns.swarmplot(x= 'treatment', y= 'colonies_pre', data=data, ax=ax1) ax1.set_title( 'pre-treatment') ax1.set_xlabel( '') ax1.set_ylabel( 'colonies') ax2 = fig.add_subplot( 2, 1, 2) sns.swarmplot(x= 'treatment', y= 'colonies_post', data=data, ax=ax2) ax2.set_title( 'post-treatment') ax2.set_ylabel( 'colonies') ax2.set_ylim(ax1.get_ylim()) plt.tight_layout() returnfigfig = plot_colonies_data()plt.show()
说明 计数是泊松分布。
withpm.Model() aspoisson_estimation: mu_pre = pm.DiscreteUniform( 'pre_mus', lower= 0, upper= 10000,shape=len(treatment_order)) pre_mus = mu_pre[data[ 'treatment_idx'].values] # fancy indexing!!pre_counts = pm.Poisson( 'pre_counts', mu=pre_mus,observed=data[ 'colonies_pre']) mu_post = pm.DiscreteUniform( 'post_mus', lower= 0, upper= 10000,shape=len(treatment_order)) post_mus = mu_post[data[ 'treatment_idx'].values] # fancy indexing!!post_counts = pm.Poisson( 'post_counts', mu=post_mus, observed=data[ 'colonies_post']) perc_change = pm.Deterministic( 'perc_change', 100* (mu_pre - mu_post) / mu_pre) MCMC Inference Button (TM)withpoisson_estimation: poisson_trace = pm.sample( 20000) pm.traceplot(poisson_trace, varnames=[ 'pre_mus', 'post_mus'])plt.show()
|
