Skip to content

Instantly share code, notes, and snippets.

@tyarkoni
Last active February 23, 2016 02:15
Show Gist options
  • Save tyarkoni/e20ab0105a1ddb4a1648 to your computer and use it in GitHub Desktop.
Save tyarkoni/e20ab0105a1ddb4a1648 to your computer and use it in GitHub Desktop.
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