fMRI example


In [39]:
import numpy as np
from scipy.io import loadmat
from os.path import join
import glob
import pandas as pd
import tensorflow as tf 
from sklearn import svm
from sklearn.metrics import accuracy_score

In [59]:
def center_normalize(x):
    x_mean_rows = np.mean(x,1).reshape(x.shape[0],1)
    x_std_rows = np.std(x,1).reshape(x.shape[0],1)
    return (x - x_mean_rows) / x_std_rows

def onehot(y):
    ynp=np.array(y)
    y_onehot=[0]*len(ynp)
    for i,j in enumerate(ynp):
        y_onehot[i]=[0]*ynp.max()
        y_onehot[i][j-1]=1
        
    return y_onehot

def model_svm(x_train,y_train,x_test,y_test):
    C = 1  # SVM regularization parameter
    svc = svm.LinearSVC(C=C).fit(x_train, y_train)
    y_pred = svc.predict(x_test)
    return accuracy_score(y_test, y_pred)

def model_tf_regression(x_train,y_train,x_test,y_test):
    tf.reset_default_graph()
    sess = tf.InteractiveSession()

    x = tf.placeholder(tf.float32, shape=[None, x_train.shape[1]])
    y_ = tf.placeholder(tf.float32, shape=[None, 3])


    W = tf.Variable(tf.zeros([x_train.shape[1],3]))
    b = tf.Variable(tf.zeros([3]))

    sess.run(tf.global_variables_initializer())

    y = tf.matmul(x,W) + b

    cross_entropy = tf.reduce_mean(
        tf.nn.softmax_cross_entropy_with_logits(labels=y_, logits=y))
    
    ##train
    train_step = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy)
    for i in range(1000):
        train_step.run(feed_dict={x: x_train.tolist(), y_: onehot(y_train)})
    
    #tf.argmax gives an index of the highest entry in a tensor along some axis
    correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))

    #we can take this list of booleans and calculate the fraction correct
    accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
    
    return accuracy.eval(feed_dict={x: x_test.tolist(), y_: onehot(y_test)})
    sess.close()

In [60]:
def model_tf_nn(x_train,y_train,x_test,y_test):
    tf.reset_default_graph()

    n_inputs = x_train.shape[1]
    n_hidden1 = 10
    n_outputs = 3
    learning_rate = 0.01


    def neuron_layer(X, n_neurons, name, activation=None):
        with tf.name_scope(name):
            n_inputs = int(X.get_shape()[1])
            stddev = 1 / np.sqrt(n_inputs)
            init = tf.truncated_normal((n_inputs, n_neurons), stddev=stddev)
            W = tf.Variable(init, name="weights")
            b = tf.Variable(tf.zeros([n_neurons]), name="biases")
            Z = tf.matmul(X, W) + b
            if activation=="relu":
                return tf.nn.relu(Z)
            else:
                return Z
    X = tf.placeholder(tf.float32, shape=(None, n_inputs), name="X")
    y = tf.placeholder(tf.int64, shape=(None), name="y")

    with tf.name_scope("dnn"):
        hidden1 = neuron_layer(X, n_hidden1, "hidden1", activation="relu")
        logits = neuron_layer(hidden1, n_outputs, "output")

    with tf.name_scope("loss"):
        xentropy = tf.nn.sparse_softmax_cross_entropy_with_logits(labels=y, logits=logits)
        loss = tf.reduce_mean(xentropy, name="loss")

    with tf.name_scope("train"):
        optimizer = tf.train.GradientDescentOptimizer(learning_rate)
        training_op = optimizer.minimize(loss)

    with tf.name_scope("eval"):
        correct = tf.nn.in_top_k(logits, y, 1)
        accuracy = tf.reduce_mean(tf.cast(correct, tf.float32))

    init = tf.global_variables_initializer()
    saver = tf.train.Saver()
    
    n_epochs = 50
    acc_test_high = 0.0

    sess = tf.InteractiveSession()

    sess.run(init)
    for epoch in range(n_epochs):
        for i in range(x_train.shape[0]):
            x_data = x_train[i,:].reshape([1,x_train.shape[1]])
            y_data = np.array([y_train[i]])
            sess.run(training_op, feed_dict={X: x_data, y: y_data })
        acc_train = accuracy.eval(feed_dict={X: x_train, y: y_train})
        acc_test = accuracy.eval(feed_dict={X: x_test, y: y_test})
        if acc_test > acc_test_high:
            acc_test_high = acc_test
#         print(epoch, "Train accuracy:", acc_train, "Test accuracy:", acc_test)

    # save_path = saver.save(sess, "./one_layer.ckpt")
    return acc_test_high
    sess.close()

In [61]:
#data dir test 
data_dir = '../data/greco_mri'
behav_dir = '../data/greco_behav'

#rois
rois = ['left_CA1','right_CA1','left_DG','right_DG'];

#subjec list
subjects = ['S1_A','S2_B','S3_A','S4_A','S5_A','S6_A','S7_A','S8_B',
            'S9_A','S10_B','S11_B','S12_B','S13_B','S14_B','S15_B',
            'S16_A','S21_B','S22_B','S24_A'];
subjects = np.array(subjects);

In [62]:
scores_svm = np.empty((subjects.shape[0],4))
scores_reg = np.empty((subjects.shape[0],4))
scores_nn = np.empty((subjects.shape[0],4))
for i,subj in enumerate(subjects):
    for run in [0,1,2,3]:
       
        behav = pd.read_table(join(behav_dir,subj,subj + '.txt'))

        fname = join(data_dir,subj + '_right_DG.csv')
        betas= pd.read_csv(fname,header=None)

        test_ind   = behav['run_num'] == run
        train_ind  = behav['run_num'] != run
        y_test     = behav['currCity'][test_ind].as_matrix()
        
        y_test_r = np.empty(y_test.shape)
        y_test_r[y_test == 1] = 0
        y_test_r[y_test == 2] = 1
        y_test_r[y_test == 3] = 2
        y_train    = behav['currCity'][train_ind].as_matrix()
        
        y_train_r = np.empty(y_train.shape)
        y_train_r[y_train == 1] = 0
        y_train_r[y_train == 2] = 1
        y_train_r[y_train == 3] = 2
        x_test = center_normalize(betas[test_ind].as_matrix())
        x_train = center_normalize(betas[train_ind].as_matrix())
        scores_svm[i,run]=model_svm(x_train,y_train,x_test,y_test)
        scores_reg[i,run] = model_tf_regression(x_train,y_train,x_test,y_test)
        scores_nn[i,run] = model_tf_nn(x_train,y_train_r,x_test,y_test_r)

In [72]:
print('Logistic regression (mean accuracy): ' + str(np.mean(scores_reg)))
print('SVM (mean accuracy): ' + str(np.mean(scores_svm)))
print('Neural network (mean accuracy): ' + str(np.mean(scores_nn)))


Logistic regression (mean accuracy): 0.324736842788
SVM (mean accuracy): 0.324736842105
Neural network (mean accuracy): 0.395789474837

In [76]:
from scipy import stats
print('Logistic regression (1samp ttest):')
print(stats.ttest_1samp(np.mean(scores_reg,1),.33))
print
print('SVM (1samp ttest):')     
print(stats.ttest_1samp(np.mean(scores_svm,1),.33))
print
print('Neural network (1samp ttest)')
print(stats.ttest_1samp(np.mean(scores_nn,1),.33))


Logistic regression (1samp ttest):
Ttest_1sampResult(statistic=-0.30706705741261137, pvalue=0.76231644869561455)

SVM (1samp ttest):
Ttest_1sampResult(statistic=-0.35242936701678962, pvalue=0.72860912450487225)

Neural network (1samp ttest)
Ttest_1sampResult(statistic=5.0932099220274942, pvalue=7.5913412853404732e-05)

In [ ]: