Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Simulates correlation between effect sizes of original studies and replication studies
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
g1_d_mu = 0.4
g1_d_sd = 0.4
prop_null = 0.3
n_subs = 20
n_studies = 400
grp1n = np.round(n_studies * prop_null)
grp2n = n_studies - grp1n
def select_sig(x, n, p=0.05):
t, p = ss.ttest_1samp(x, np.zeros(len(x)), axis=1)
sig_inds = np.where((p < 0.05) & ( t > 0))[0]
xt = x[sig_inds[:n], :]
d = xt.mean(1) / xt.std(1)
return (d, sig_inds)
# Null-effect study group
x = np.random.normal(size=(grp1n*50, n_subs))
d0_t1 = select_sig(x, grp1n)[0]
x0_t2 = np.random.normal(size=(grp1n, n_subs))
d0_t2 = x0_t2.mean(1) / x0_t2.std(1)
print('Group 1 r: %.2f' % np.corrcoef(d0_t1, d0_t2)[0, 1])
plt.scatter(d0_t1, d0_t2, c='b')
# Real effect study group
_d = np.random.normal(g1_d_mu, g1_d_sd, size=grp2n*10)
cov = np.eye(grp2n*10)
x1_t1 = np.random.multivariate_normal(_d, cov, size=n_subs).T
d1_t1, val_inds = select_sig(x1_t1, grp2n)
x1_t2 = np.random.multivariate_normal(_d, cov, size=n_subs).T[val_inds][:grp2n]
d1_t2 = x1_t2.mean(1) / x1_t2.std(1)
plt.scatter(d1_t1, d1_t2, c='g')
print('Group 2 r: %.2f' % np.corrcoef(d1_t1, d1_t2)[0, 1])
print('Collapsed r: %.2f' % np.corrcoef(np.r_[d0_t1, d1_t1], np.r_[d0_t2, d1_t2])[0, 1])
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment