Simulates correlation between effect sizes of original studies and replication studies

t1_t2_correlation_sim.py

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()