cwas_simulations_part.py

import sys
#sys.path.append("/Users/zarrar/Dropbox/Code/C-PAC/CPAC/cwas")
sys.path.insert(0, "/home2/data/Projects/CPAC_Regression_Test/2013-05-30_cwas/C-PAC/CPAC/cwas")
sys.path.append("/home/data/PublicProgram/epd-7.2-2-rh5-x86_64/lib/python2.7/site-packages")

import multiprocessing
import numpy as np
from mdmr import mdmr
from sklearn.metrics.pairwise import euclidean_distances


# Simulate
def dtest(n=50, d=0.0, r=0.0, model_covariate=True, niters=100, nperms=4999):
    import mkl
    mkl.set_num_threads(2)
    
    d   = float(d)
    r   = float(r)
        
    # Data/Distances
    pvals = np.zeros(niters)
    Fvals = np.zeros(niters)
    for i in xrange(niters):
        # Design
        
        ## Categorical
        gp  = np.repeat([0, 1], n/2)
        np.random.shuffle(gp)
        x   = gp*d + np.random.standard_normal(n)
        
        ## Continuous
        # see http://stackoverflow.com/questions/16024677/generate-correlated-data-in-python-3-3
        # and http://stats.stackexchange.com/questions/19367/creating-two-random-sequences-with-50-correlation?lq=1
        uncorrelated    = np.random.standard_normal((2,n))
        motion          = uncorrelated[0]
        y               = r*motion + np.sqrt(1-r**2)*uncorrelated[1]
        
        ## Design Matrix
        if model_covariate:
            design = np.vstack((np.ones(n), gp, motion)).T
        else:
            design = np.vstack((np.ones(n), gp)).T
                
        # Date
        points = np.vstack((x,y)).T
        
        # Distances
        dmat  = euclidean_distances(points)
        dmats = dmat[np.newaxis,:,:]
        
        # Only the group effect is the variable of interest
        cols = [1]
        
        # Call MDMR
        pval, Fval, _, _ = mdmr(dmats, design, cols, nperms)
        
        pvals[i] = pval
        Fvals[i] = Fval
    
    return pvals, Fvals


# We just need to run the above
# for a given n and nperms
# and varying d (effect of interest)
# and varying r but only 3 times (covariate)
# and we model or don't model r (the covariate)

n  = 50
ds = np.linspace(0,1.5,16)
rs = np.linspace(0,1,3)
niters = 1000
nperms = 999

def dtest_wrap(args):
    return dtest(*args, niters=niters, nperms=nperms)

def dtest_wrap_nocov(args):
    return dtest(*args, model_covariate=False, niters=niters, nperms=nperms)

# Get with covariate
p       = multiprocessing.Pool(processes=12)
args    = [ (n,d,r) for r in rs for d in ds ]
res1    = p.map(dtest_wrap, args)

# Get without covariate
p       = multiprocessing.Pool(processes=12)
args    = [ (n,d,r) for r in rs for d in ds ]
res2    = p.map(dtest_wrap_nocov, args)