multivariate_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5530528d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import scipy.io\n",
    "import h5py\n",
    "import numpy as np\n",
    "import pickle\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.model_selection import train_test_split\n",
    "from scipy import stats\n",
    "import antropy as ant\n",
    "filepath = '/Users/yutingzhang/Downloads/subInfo_316_attention_210226.mat'\n",
    "arrays = {}\n",
    "f = h5py.File(filepath)\n",
    "for k, v in f.items():\n",
    "    arrays[k] = np.array(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0024625a",
   "metadata": {},
   "source": [
    "***THIS IS THE MULTIVARIATE METHOD***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f8ef040",
   "metadata": {},
   "outputs": [],
   "source": [
    "#split the mat_rest and label into training and testing groups\n",
    "mat_rest = arrays['mat_rest']\n",
    "subList = arrays['subList']\n",
    "X = mat_rest\n",
    "Y = subList[1]\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.33, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5005642a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#standardize the data, fit on the training data and apply the fit to both training and testing datasets. \n",
    "scaler = StandardScaler()\n",
    "scaler.fit(X_train)\n",
    "X_train = scaler.transform(X_train)\n",
    "X_test = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35caf52c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#initiate PCA model\n",
    "#either pick the variance to be used\n",
    "pca_var = PCA(0.2) # usually use 5 components to retain 20% variance\n",
    "#or pick the number of componenets to be used\n",
    "pca_com = PCA(n_components=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dfea6f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fit the data using pca_var\n",
    "pca_var.fit(X_train)\n",
    "#look for how many components are retained to achieve 95% variance\n",
    "pca_var.n_components_\n",
    "#apply the transform to both \n",
    "X_train = pca_var.transform(X_train)\n",
    "X_test = pca_var.transform(X_test)\n",
    "#35778"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d6ff197",
   "metadata": {},
   "outputs": [],
   "source": [
    "ran = round(len(mat_rest[0])/1000)\n",
    "n_sub = len(mat_rest)\n",
    "sub_result = []\n",
    "for i in range(n_sub):\n",
    "    sub_data = mat_rest[i]\n",
    "    for j in range(ran):\n",
    "        sphere_prior = j*1000\n",
    "        if j != ran-1:\n",
    "            sphere_posterior = (j+1)*1000\n",
    "        else:\n",
    "            sphere_posterior = len(mat_rest[0])-1\n",
    "        for k in range(ran):\n",
    "            seed_prior = k*1000\n",
    "            if k != ran-1:\n",
    "                seed_posterior = (k+1)*1000\n",
    "            else:\n",
    "                seed_posterior = len(mat_rest[0])-1\n",
    "            X = sub_data[sphere_prior:sphere_posterior]\n",
    "            Y = sub_data[seed_prior:seed_posterior]\n",
    "            X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=42)\n",
    "            pca_com = PCA(n_components=3)\n",
    "            #fit the data using pca_var\n",
    "            pca_var.fit(X_train)\n",
    "            #look for how many components are retained to achieve 95% variance\n",
    "            pca_var.n_components_\n",
    "            #apply the transform to both \n",
    "            X_train = pca_var.transform(X_train)\n",
    "            X_test = pca_var.transform(X_test)\n",
    "            logisticRegr = LogisticRegression(solver = 'lbfgs')\n",
    "            logisticRegr.fit(X_train,Y_train)\n",
    "            score = logisticRegr.score(X_test,Y_test)\n",
    "            sub_result.append(score)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c1c96bf",
   "metadata": {},
   "source": [
    "****CPM METHOD STARTS FROM HERE!!!****"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e694e9",
   "metadata": {},
   "source": [
    "**construct connectivity matrix and select the edge**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba5cc104",
   "metadata": {},
   "source": [
    "*1. read in data, a matrix of 268*268 for each subject*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cba7182a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def read_file(path):\n",
    "    filepath = path\n",
    "    arrays_conmat = {}\n",
    "    f = h5py.File(filepath)\n",
    "    \n",
    "    for k, v in f.items():\n",
    "        arrays_conmat[k] = np.array(v)\n",
    "        \n",
    "    return arrays_conmat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7411c2c3",
   "metadata": {},
   "source": [
    "*2. Extract the relevant data into numpy array*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05b6fe11",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_data(arrays_conmat):\n",
    "    ##extract the phenotype data\n",
    "    pheno = arrays_conmat['subList'][1]\n",
    "    ##calculate the mean of the resting data\n",
    "    matrix = arrays_conmat['conMatRun268']\n",
    "    rest = matrix[0:4] #first four columns are resting state. \n",
    "    shape = np.shape(matrix)\n",
    "    new = []\n",
    "    \n",
    "    for sub in range(shape[1]):\n",
    "        sub_data = rest[:,sub,:,:]\n",
    "        sub_mean = np.mean(sub_data,axis=0) #Take the average between the four runs. \n",
    "        new.append(sub_mean) #the size of new is (316,268,268)\n",
    "        \n",
    "    return pheno, new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50689ff7",
   "metadata": {},
   "source": [
    "*Now, we have pheno, which is the pheno data (size:316, one for each subject), and new, which is the matrix (size: (316,268,268), one matrix for each of 316 subject)*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5da0bac5",
   "metadata": {},
   "source": [
    "*3. Calculate the correlation matrix and p-value matrix*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0e883f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_gen(pheno,new):\n",
    "    new_new = np.asarray(new)\n",
    "    cor_matrix = []\n",
    "    p_matrix = []\n",
    "    \n",
    "    for i in range(268):\n",
    "        for j in range(268):\n",
    "            data = new_new[:,i,j]\n",
    "            cor = stats.pearsonr(data,pheno)\n",
    "            cor_matrix.append(cor[0])\n",
    "            p_matrix.append(cor[1])\n",
    "            \n",
    "    cor_matrix = np.asarray(cor_matrix)\n",
    "    p_matrix = np.asarray(p_matrix)\n",
    "    \n",
    "    cor_matrix_resize = np.resize(cor_matrix,(268,268))\n",
    "    cor_matrix_resize = np.nan_to_num(cor_matrix_resize)\n",
    "    \n",
    "    p_matrix_resize = np.resize(p_matrix,(268,268))\n",
    "    p_matrix_resize = np.nan_to_num(p_matrix_resize)\n",
    "    \n",
    "    return cor_matrix_resize,p_matrix_resize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89cee44a",
   "metadata": {},
   "source": [
    "*4. With the correlation matrix and p-value matrix, we need to select positive and negative edges. If specific edge is greater than 0 and have p-value less than 0.05, it will be marked with 1, otherwise 0, in the pos_matrix. Likewise, if specific edge is less than 0 and have p-value less than 0.05, it will be marked with 1, otherwise 0, in the neg_matrix.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ddd17e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def edge_selection(cor_matrix_resize,p_matrix_resize,target_size):\n",
    "    \n",
    "    #create positive edge matrix and negative edge matrix\n",
    "    pos_matrix = []\n",
    "    neg_matrix = []\n",
    "    \n",
    "    assert np.shape(cor_matrix_resize) == target_size, \"Wrong size!\"\n",
    "    assert np.shape(p_matrix_resize) == target_size, \"Wrong size!\"\n",
    "    \n",
    "    for i in range(np.shape(cor_matrix_resize)[0]):\n",
    "        for j in range(np.shape(cor_matrix_resize)[0]):\n",
    "            cor = cor_matrix_resize[i,j]\n",
    "            p = p_matrix_resize[i,j]\n",
    "            if cor>0 and p<0.05:\n",
    "                pos_matrix.append(1)\n",
    "                neg_matrix.append(0)\n",
    "            elif cor<0 and p<0.05: \n",
    "                pos_matrix.append(0)\n",
    "                neg_matrix.append(1)\n",
    "            else:\n",
    "                pos_matrix.append(0)\n",
    "                neg_matrix.append(0)\n",
    "                \n",
    "    pos_matrix = np.reshape(pos_matrix,target_size)\n",
    "    neg_matrix = np.reshape(neg_matrix,target_size)\n",
    "    \n",
    "    return pos_matrix,neg_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "715a327e",
   "metadata": {},
   "source": [
    "*5. Now, we have two matrices, pos_mat and neg_mat, composed of 1 and 0. We will do an elementwise multiplication to select the positive edge and negative edge on the matrix*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99b20464",
   "metadata": {},
   "outputs": [],
   "source": [
    "def element_multiply(all_sub_data,pos_matrix,neg_matrix):\n",
    "    all_pos = []\n",
    "    all_neg = []\n",
    "    for sub in all_sub_data:\n",
    "        sub_data_pos = np.multiply(sub,pos_matrix)\n",
    "        sub_data_neg = np.multiply(sub,neg_matrix)\n",
    "        sub_data_pos = np.reshape(sub_data_pos,(268,268))\n",
    "        sub_data_neg = np.reshape(sub_data_neg,(268,268))\n",
    "        all_pos.append(sub_data_pos)\n",
    "        all_neg.append(sub_data_neg)\n",
    "    return all_pos,all_neg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf106da",
   "metadata": {},
   "source": [
    "*6. Now we have the matrix of positive edge and matrix of negative edge. We want to calculate the summary statistics. Therefore, for each subject, there will be two numbers, with one corresponding to positive edges' summary stats and the other corresponding to negative edges' summary stats.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cabcf98",
   "metadata": {},
   "outputs": [],
   "source": [
    "def summary_stats(all_pos,all_neg):\n",
    "    pos_sum = []\n",
    "    neg_sum = []\n",
    "    for pos_mat in all_pos:\n",
    "        pos = np.sum(pos_mat)\n",
    "        pos_sum.append(pos)\n",
    "    for neg_mat in all_neg:\n",
    "        neg = np.sum(neg_mat)\n",
    "        neg_sum.append(neg)\n",
    "    return pos_sum, neg_sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33726f04",
   "metadata": {},
   "outputs": [],
   "source": [
    "def main():\n",
    "    size = (268,268) #set up the size of each matrix for each subject\n",
    "    filepath = '/Users/yutingzhang/Downloads/conmat_HCP_run_316_201216.mat'\n",
    "    print('finished with filepath')\n",
    "    array = read_file(filepath) #read in the file\n",
    "    pheno, sub_mat = extract_data(array) #extract the matrix and phenotype data\n",
    "    print('finished with a')\n",
    "    cor_mat, p_mat = matrix_gen(pheno,sub_mat) #extract the correlation matrix and p-value matrix for all subjects\n",
    "    pos_mat, neg_mat = edge_selection(cor_mat,p_mat,size) #edge selection with two matrices, one for pos, one for neg\n",
    "    print('finished with b')\n",
    "    sub_pos, sub_neg = element_multiply(sub_mat,pos_mat,neg_mat) #using the matrix to apply to the actual subject data to select that\n",
    "    sub_sum_pos,sub_sum_neg = summary_stats(sub_pos,sub_neg) #calculate the summary statistics\n",
    "    print('finished with c')\n",
    "    return sub_sum_pos,sub_sum_neg,pheno"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b44b8c06",
   "metadata": {},
   "source": [
    "****Implementation****"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d8c74ea",
   "metadata": {},
   "source": [
    "Note: Try leave-one-out ordinary least square regression, ridge regression, lasso regression, Elastic-Net, LARS-lasso, stochastic gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85e2bd70",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_sum_pos,sub_sum_neg,pheno = main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1347e5fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(sub_sum_neg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd3e4f6",
   "metadata": {},
   "source": [
    "1. Ordinary least sqaure regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2c763c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "def OLS(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = LinearRegression().fit(X_train, Y_train)\n",
    "        prediction = reg.predict(np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg))))\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    return original,predict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3385f668",
   "metadata": {},
   "source": [
    "Evaluate the result using correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e313662",
   "metadata": {},
   "outputs": [],
   "source": [
    "original, predict = OLS(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc89ae37",
   "metadata": {},
   "source": [
    "2. Ridge regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d324567",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "def ridgeRegression(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    original = []\n",
    "    predict = []\n",
    "\n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = Ridge(alpha = 0.0).fit(X_train, Y_train)\n",
    "        prediction = reg.predict(np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg))))\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    return original,predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c9da79a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "original, predict = ridgeRegression(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c0a506",
   "metadata": {},
   "source": [
    "3. Lasso regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d2b7a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "def lassoRegression(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    \n",
    "    test_image = []\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = Lasso(alpha = 1.0).fit(X_train, Y_train)\n",
    "        test_x = np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg)))\n",
    "        prediction = reg.predict(test_x)\n",
    "        test_image.append(test_x)\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    \n",
    "    return original,predict,test_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb5ef25e",
   "metadata": {},
   "outputs": [],
   "source": [
    "original, predict, test_image = lassoRegression(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06b4a1b3",
   "metadata": {},
   "source": [
    "4. Elastic-Net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b458bee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import ElasticNet\n",
    "def elasticNet(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    \n",
    "    test_image = []\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        print(np.shape(loo_sub_sum_pos))\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = ElasticNet(alpha = 1.0).fit(X_train, Y_train)\n",
    "        test_x = np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg)))\n",
    "        prediction = reg.predict(test_x)\n",
    "        test_image.append(test_x)\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    \n",
    "    return original,predict,test_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0b6904",
   "metadata": {},
   "outputs": [],
   "source": [
    "original, predict, test_image = elasticNet(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0b08c78",
   "metadata": {},
   "source": [
    "5. LARS-lasso\n",
    "\n",
    "Note: It is a Linear Model trained with an L1 prior as regularizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5184a947",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoLars\n",
    "\n",
    "def lassoLars(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    \n",
    "    test_image = []\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = LassoLars(alpha = 1.0).fit(X_train, Y_train)\n",
    "        test_x = np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg)))\n",
    "        prediction = reg.predict(test_x)\n",
    "        test_image.append(test_x)\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    \n",
    "    return original,predict,test_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db1fea88",
   "metadata": {},
   "outputs": [],
   "source": [
    "original, predict, test_image = lassoLars(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1ae16a",
   "metadata": {},
   "source": [
    "*plot to see whether the prediction and original are perfectly negatively correlated*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd66154",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set()\n",
    "dataFrame = pd.DataFrame(\n",
    "    {'original': original,\n",
    "     'prediction': predict\n",
    "    }) #create a dataframe for the data\n",
    "ax = sns.scatterplot(data=dataFrame, x=\"original\", y=\"prediction\")\n",
    "ax.set_title('Relationship between original and prediction using LassoLars',\n",
    "             fontdict= { 'fontsize': 15, 'fontweight':'bold'});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28fcacde",
   "metadata": {},
   "source": [
    "6. Stochastic gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a519a761",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import SGDRegressor\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "def SGDRegress(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    \n",
    "    test_image = []\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = make_pipeline(StandardScaler(), SGDRegressor(alpha=1.0)).fit(X_train, Y_train)\n",
    "        '''\n",
    "        Args\n",
    "        \n",
    "        1. loss\n",
    "        The ‘squared_loss’ refers to the ordinary least squares fit. \n",
    "        ‘huber’ modifies ‘squared_loss’ to focus less on getting outliers correct by switching from \n",
    "        squared to linear loss past a distance of epsilon. ‘epsilon_insensitive’ ignores errors less \n",
    "        than epsilon and is linear past that; this is the loss function used in SVR. \n",
    "        ‘squared_epsilon_insensitive’ is the same but becomes squared loss past a tolerance of epsilon.\n",
    "        \n",
    "        2. penalty\n",
    "        The penalty (aka regularization term) to be used. Defaults to ‘l2’ which is the standard regularizer \n",
    "        for linear SVM models. ‘l1’ and ‘elasticnet’ might bring sparsity to the model (feature selection) not \n",
    "        achievable with ‘l2’.\n",
    "        \n",
    "        3. alphafloat, default=0.0001\n",
    "        Constant that multiplies the regularization term. The higher the value, the stronger the regularization. \n",
    "        Also used to compute the learning rate when set to learning_rate is set to ‘optimal’.\n",
    "        \n",
    "        \n",
    "        '''\n",
    "        test_x = np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg)))\n",
    "        prediction = reg.predict(test_x)\n",
    "        test_image.append(test_x)\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    \n",
    "    return original,predict,test_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f42e971a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "original, predict, test_image = SGDRegress(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2637c17",
   "metadata": {},
   "source": [
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c794e19",
   "metadata": {},
   "source": [
    "**Use correlation to evaluate different models with different alpha values (the value that control the regularization/penalty)**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8d8417f",
   "metadata": {},
   "source": [
    "*Correlation of each model with different alpha*\n",
    "\n",
    "|  | LSR | Ridge | Lasso | ElasticNet | LARS-Lasso | Stochastic gradient descent | \n",
    "| --- | --- | --- | --- | --- | --- | --- |\n",
    "| alpha = 0.0 | 0.5894942746983071 |0.5894942746983067 | 0.5894942746983067 | 0.5894942746983067 | 0.5894942746983068 | 0.5894944595848238 |\n",
    "| alpha = 0.1 | / | 0.5894942657346615 | 0.5894519058277794 | 0.5894719470529471 | 0.5169165935569177 | 0.5891520624645177 | \n",
    "| alpha = 0.2 | / | 0.5894942567660058 | 0.5894032548024233 | 0.5894471287093305 | -1.0 | 0.5895695339551591 |\n",
    "| alpha = 0.3 | / | 0.5894942477923397 | 0.589348130112169 | 0.5894197471228471 | -1.0 | 0.5887090551304505 |\n",
    "| alpha = 0.4 | / | 0.589494238813663 | 0.5892863911737424 | 0.5893897830994119 | -1.0 | 0.5884137119056182 |\n",
    "| alpha = 0.5 | / | 0.5894942298299759 | 0.5892178948780954 | 0.5893572173765593 | -1.0 | 0.5881458294354422 |\n",
    "| alpha = 0.6 | / | 0.5894942208412788 | 0.5891424955577386 | 0.5893220306242982 | -1.0 | 0.5874995443343022 |\n",
    "| alpha = 0.7 | / | 0.5894942118475708 | 0.5890600449542791 | 0.5892842034459883 | -1.0 | 0.5882114874616413 |\n",
    "| alpha = 0.8 | / | 0.5894942028488538 | 0.5889703921861977 | 0.5892437163792422 | -1.0 | 0.5873394456538039 |\n",
    "| alpha = 0.9 | / | 0.5894941938451262 | 0.5888733837169156 | 0.5892005498968458 | -1.0 | 0.5876962520402159 |\n",
    "| alpha = 1.0 | / | 0.589494184836389 | 0.588768863323186 | 0.5891546844077098 | -1.0 | 0.5865179031388511 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74156ddb",
   "metadata": {},
   "source": [
    "Question: \n",
    "\n",
    "1. what if the summary statics use entropy rather than simple summary?\n",
    "\n",
    "2. What if the original coding scheme is changed to sparse encoding? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f9eaeae",
   "metadata": {},
   "source": [
    "**Using entropy as summary statistics**\n",
    "\n",
    "The following steps are divided into\n",
    "\n",
    "1) calculate the entropy stats as the summary stats\n",
    "\n",
    "2) call the main function to calculate the steps up to the summary stats \n",
    "\n",
    "3) Use least sqaure multiple regression to do the prediction\n",
    "\n",
    "4) Use correlation to evaluate whether the model is a good fit. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb19d142",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyitlib import discrete_random_variable as drv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0159d54c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy_stats(all_pos,all_neg):\n",
    "    pos_sum = []\n",
    "    neg_sum = []\n",
    "    size = 268\n",
    "    for pos_mat in all_pos:\n",
    "        reshaped_pos = np.reshape(pos_mat, (268 * 268, 1))\n",
    "        droped_pos = reshaped_pos[reshaped_pos!=0.]\n",
    "        droped_pos = droped_pos[droped_pos!=-0.] #it doesn't matter whether it's zscored or not. \n",
    "        pos = ant.sample_entropy(droped_pos)\n",
    "        '''\n",
    "        Entropy candidate to consider: app_entropy, perm_entropy, sample_entropy\n",
    "        '''\n",
    "        pos_sum.append(pos)\n",
    "        \n",
    "    for neg_mat in all_neg:\n",
    "        reshaped_neg = np.reshape(neg_mat, (268 * 268, 1))\n",
    "        droped_neg = reshaped_neg[reshaped_neg!=0.]\n",
    "        droped_neg = droped_neg[droped_neg!=-0.]\n",
    "        neg = ant.sample_entropy(droped_neg)\n",
    "        neg_sum.append(neg)\n",
    "        \n",
    "    return pos_sum, neg_sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4ccc17a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy_main():\n",
    "    size = (268,268) #set up the size of each matrix for each subject\n",
    "    filepath = '/Users/yutingzhang/Downloads/conmat_HCP_run_316_201216.mat'\n",
    "    print('finished with filepath')\n",
    "    array = read_file(filepath) #read in the file\n",
    "    pheno, sub_mat = extract_data(array) #extract the matrix and phenotype data\n",
    "    print('finished with a')\n",
    "    cor_mat, p_mat = matrix_gen(pheno,sub_mat) #extract the correlation matrix and p-value matrix for all subjects\n",
    "    pos_mat, neg_mat = edge_selection(cor_mat,p_mat,size) #edge selection with two matrices, one for pos, one for neg\n",
    "    print('finished with b')\n",
    "    sub_pos, sub_neg = element_multiply(sub_mat,pos_mat,neg_mat) #using the matrix to apply to the actual subject data to select that\n",
    "    sub_sum_pos,sub_sum_neg = entropy_stats(sub_pos,sub_neg) #calculate the summary statistics\n",
    "    print('finished with c')\n",
    "    return sub_sum_pos,sub_sum_neg,pheno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dd80402",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sub_sum_pos,sub_sum_neg,pheno = entropy_main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80fdc1f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "def OLS(sub_sum_pos,sub_sum_neg,pheno):\n",
    "    original = []\n",
    "    predict = []\n",
    "    \n",
    "    for i in range(np.shape(sub_sum_pos)[0]):\n",
    "        print(\"On subject %s\" %i)\n",
    "        loo_sub_sum_pos = np.delete(sub_sum_pos,[i])\n",
    "        single_sub_sum_pos = sub_sum_pos[i]\n",
    "        loo_sub_sum_neg = np.delete(sub_sum_neg,[i])\n",
    "        single_sub_sum_neg = sub_sum_neg[i]\n",
    "        loo_pheno = np.delete(pheno,[i])\n",
    "        single_sub_pheno = pheno[i]\n",
    "        X_train = np.vstack((np.reshape(loo_sub_sum_pos,(1,315)), np.reshape(loo_sub_sum_neg,(1,315))))\n",
    "        X_train = np.transpose(X_train)\n",
    "        Y_train = loo_pheno\n",
    "        reg = LinearRegression().fit(X_train, Y_train)\n",
    "        prediction = reg.predict(np.transpose(np.vstack((single_sub_sum_pos,single_sub_sum_neg))))\n",
    "        original.append(single_sub_pheno)\n",
    "        predict.append(prediction[0])\n",
    "    return original,predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fecc4ff",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "original, predict = OLS(sub_sum_pos,sub_sum_neg,pheno)\n",
    "cor = np.corrcoef(np.asarray(original),np.asarray(predict))[0][1]\n",
    "print(\"Correlation of the predicted and original score is %s.\" %cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4294bfa",
   "metadata": {},
   "source": [
    "*Correlation of least sqaure model with different entropy as the summary stats*\n",
    "\n",
    "|  | app_entropy | perm_entropy | sample_entropy | \n",
    "| --- | --- | --- | --- | \n",
    "| cor | 0.01641536189675583 | 0.030878640150066024 | -0.06188538442233747 |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ca51f6a",
   "metadata": {},
   "source": [
    "*Ignore the following*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "093dbc2c",
   "metadata": {},
   "source": [
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c244c09a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#apply logistic regression\n",
    "logisticRegr = LogisticRegression(solver = 'lbfgs')\n",
    "logisticRegr.fit(X_train,Y_train)\n",
    "logisticRegr.predict(X_test[0:100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b286b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#test the results for logistic regression\n",
    "logisticRegr.score(X_test,Y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87e07fc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#apply linear model, ordinary least squares\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42af47a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#apply "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28888ffe",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fit the data using pca_com\n",
    "pca_com.fit(X_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0250667d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(n_components = 2)\n",
    "Principle_Component = pca.fit_transform(mat_rest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ccbd98b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataframe = pd.DataFrame(data = Principle_Component, columns = ['Principal Component 1','Principal Component 2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32a699e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d0df864",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5a3ed92",
   "metadata": {},
   "outputs": [],
   "source": [
    "#save pickle file\n",
    "\n",
    "mat_rest = arrays['mat_rest']\n",
    "pickle_out = open(\"/Users/yutingzhang/Downloads/attention_array.pickle\",\"wb\")\n",
    "pickle.dump(mat_rest, pickle_out)\n",
    "pickle_out.close()\n",
    "\n",
    "sub_list = arrays['subList']\n",
    "pickle_out = open(\"/Users/yutingzhang/Downloads/attention_sublist.pickle\",\"wb\")\n",
    "pickle.dump(sub_list, pickle_out)\n",
    "pickle_out.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34b64854",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load in pickle file as a list\n",
    "pickle_file = open('/Users/yutingzhang/Downloads/attention_array.pickle','rb')\n",
    "objects = []\n",
    "while True:\n",
    "    try:\n",
    "        objects.append(pickle.load(pickle_file))\n",
    "    except EOFError:\n",
    "        break\n",
    "pickle_file.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4915c272",
   "metadata": {},
   "outputs": [],
   "source": [
    "pointer = arrays['id_name'][0][0]\n",
    "arrays['id_name'][pointer]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}