<jupyter_start><jupyter_text>Opioid DataHW 2 Part 1 - Dimensionality Reduction.  Use one summary vector per patient, specifically the monthly average per patient.    PreprocessingPatient files are in one of two directories: R or NR.  Each patient is represented by one CSV file.  Each row of each CSV contains readings from one day.    Here, we load each patient average across all days.<jupyter_code>from os import listdir
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA, KernelPCA
from sklearn.preprocessing import StandardScaler
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.cluster import FeatureAgglomeration
from sklearn.random_projection import GaussianRandomProjection
from sklearn.manifold import TSNE
from sklearn.model_selection import cross_validate
import matplotlib.pyplot as plt
from matplotlib import colors
mycmap = colors.ListedColormap(['red','blue'])  # list color for label 0 then 1
np.set_printoptions(precision=2)
pathR='data/ChunkedData_R/'
pathN='data/ChunkedData_NR/'
filesR = listdir(pathR)
filesN = listdir(pathN)
print(len(filesR)," files of type R")
print(len(filesN)," files of type N")
CLASS_SEPARATOR=13  # data[:13] vs data[13:]
WITH_VARIANCE_COLUMNS=True   # Use mean and variance per patient
# Read one CSV file. 
# Create a Pandas data frame.
# Drop the date column.
def file_mean (filepath):
    mydata = pd.read_csv(filepath)
    # Drop the date column.
    mydata = mydata.drop('Date',axis=1) 
    mcols = [name+'M' for name in mydata.columns]
    vcols = [name+'V' for name in mydata.columns]
    # Transpose column of mean values into a row.
    mymean = mydata.mean(axis=0).to_frame().T
    mymean.columns=mcols
    myvar = mydata.var(axis=0).to_frame().T
    myvar.columns=vcols
    return mymean,myvar
# Read directory of CSV files (R or NR). 
# Create one dataframe representing all files.
# Retain only one row per file = column averages.
def mean_per_file (directory,variance):
    files = listdir(directory)
    means = pd.DataFrame()
    for fp in files:
        dfp = directory+fp
        m,v = file_mean(dfp)
        m_v = pd.concat([m,v],axis=1)
        # Let Pandas number the rows sequentially.
        means = means.append(m_v,ignore_index=True)
    return means
meansR = mean_per_file(pathR,WITH_VARIANCE_COLUMNS)
meansR.describe()
meansN = mean_per_file(pathN,WITH_VARIANCE_COLUMNS)
meansN.describe()
# Combine all data (R and NR) into one data frame.
# Label R = positive = 1 = blue.
# Label NR = negative = 0 = red.
def make_labels(positives,negatives):
    rows = positives.shape[0]
    labelsP = pd.DataFrame(np.ones(rows,dtype=np.int8))  # one = blue
    rows = negatives.shape[0]
    labelsN = pd.DataFrame(np.zeros(rows,dtype=np.int8))  # zero = red
    labelsAll = pd.concat((labelsP,labelsN),ignore_index=True)
    return labelsAll
# If we had all the data and their labels in one place, we could do this:
#allData.groupby('R_NR').mean().T
# Instead, we have unlabeled data in meansAll orgranized by position (R then NR)
meansAll = pd.concat((meansR,meansN),ignore_index=True)
print("(rows,columns)=",meansAll.shape)
labels_df = make_labels(meansR,meansN)
#print(labels_df.T)<jupyter_output>(rows,columns)= (40, 518)<jupyter_text>Scaling and covarianceNormalize by subtracting the column mean from every column value.  Since columns have widely different numerical ranges,   also normalize by making each column have unit variance.  Note: without normalization, the covariance plot would be all black except for the few features with large absolute values.<jupyter_code># Standardize features by shifting the mean to zero and scaling to unit variance.
# Subtract the mean and divide by the std.dev: z = (x - u) / s
def scale_features(X):
    s = StandardScaler()
    z = s.fit_transform(X)
    return z
scaledMeans = scale_features(meansAll)
column_names = meansAll.columns
print("Column Names:")
print(column_names)
print("Scaled Means & Variances:")
scaledMeans
def plot_covariance(X):
    cv=X.T.dot(X)/len(X)
    plt.imshow(cv, cmap='hot', interpolation='nearest')
    plt.show()
plot_covariance(scaledMeans)
pathR = None
pathN = None
filesR = None
filesN = None
meansR = None
meansN = None
meansAll = None<jupyter_output><empty_output><jupyter_text>Principal Component Analysis (PCA)PCA is an example of unsupervised learning.PCA does dimensionality reduction by a linear transformationto orthogonal axes where each successive axis captures most of the remaining variance.  It is important to normalize first.Otherwise, most variance will be explained by those columns with large absoulte values.That was our mistake in HW 1, when 67% of variance was explained by PC1.  We will examine prinicpal components 1, 2, 3.<jupyter_code>def show_PC_variance_explained(model):
    # Assumes at least 3 principal components, or else crashes.
    e1,e2,e3=model.explained_variance_ratio_[:3]*100.0
    print("Variance explained by PC1=%.2f%% PC2=%.2f%% PC3=%.2f%%"%(e1,e2,e3))
def show_PC_eigenvalues(model):
    e1,e2,e3=model.lambdas_[:3]
    print("Eigenvalues EV1=%.2f EV2=%.2f EV3=%.2f"%(e1,e2,e3))
def construct_PCA():
    # Use the PCA class from sklearn.
    # Linear dimensionality reduction using Singular Value Decomposition.
    # This is unsupervised learning (but we'll use labels for visualization).
    # Every transform returns a COPY of the data; see copy parameter.
    # Does centering (setting mean = 0); no way to disable this.
    # Does not do scaling (setting variance = 1); no way to enable this.
    # Does NOT do whitening (setting variance = covariance = 1); see whitening parameter.
    # When n_components = None, PCA uses min(features,instances).
    # Can also be set to 'mle' or min percent of variance to explain.
    # For approximation that computes faster, use svd_solver='randomized'.
    model = PCA()
    return model
print("PCA of total feature space:")
pca=construct_PCA()
pca.fit(scaledMeans)
show_PC_variance_explained(pca)
P_transformed=pca.fit_transform(scaledMeans)
print("Shape of transformed data",P_transformed.shape)<jupyter_output>PCA of total feature space:
Variance explained by PC1=27.42% PC2=11.37% PC3=6.55%
Shape of transformed data (40, 40)<jupyter_text>Linear Discriminant Analysis (LDA)LDA is a classification algorithm. LDA is an example of supervised learning.LDA finds a linear decision boundary.LDA can be used for dimensionality reduction.It projects the data onto some number of axes that most discriminate between the classes.The maximum number of dimensions is n_classes-1.<jupyter_code>def construct_LDA():
    # A classifier with a linear decision boundary, 
    # generated by fitting class conditional densities 
    # to the data and using Bayes’ rule.
    # Assumes a Gaussian density for each class.
    # Assumes that all classes share the same covariance matrix.
    # Solver = 'svd' (default), 'lsqr', or 'eigen'.
    # By default, priors are inferred from inputs.
    # By default, n_components = None, and LDA uses min(features,classes-1).  
    # It is possible to ask LDA for the means, classes, priors,
    # variance explained, decision boundary line, and within-class covariance.
    return LinearDiscriminantAnalysis()
# Must use ravel() to convert shape (40,1) to (40,).
labels_ravel=np.ravel(labels_df)
lda=construct_LDA()
training = cross_validate(lda, scaledMeans, labels_ravel, cv=5)
training['test_score']
print("LDA of total feature space:")
lda.fit(scaledMeans,labels_ravel)
# The linear decision boundary is a (n_features)-dimensional vector.
print("Shape of decision boundary",lda.coef_.shape)
L_transformed = lda.transform(scaledMeans)
# The transformed data is (n_classes-1)-dimensional i.e. 1D.
print("Shape of transformed data",L_transformed.shape)<jupyter_output>LDA of total feature space:
Shape of decision boundary (1, 518)
Shape of transformed data (40, 1)<jupyter_text>HW2.1.a: Comparison of PCA vs LDAThe LDA has more information, specifically, the labels.Also, since we have only two classes, the LDA must perform maximum separation on one axis.Thus, it is not surprising that the LDA achieves greater separation compared to the first principal component.<jupyter_code>fig1 = plt.figure(figsize=(10,2)) 
ax1 = fig1.add_subplot(121)
ax1.set_title('PCA first axis')
# Extract list of PC1 coordinate for each of 40 intances.
pc1_coords = [c[0] for c in P_transformed]
ax1.hist(pc1_coords[:CLASS_SEPARATOR],histtype='step',color='blue') 
ax1.hist(pc1_coords[CLASS_SEPARATOR:],histtype='step',color='red')  
ax2 = fig1.add_subplot(122)
ax2.set_title('LDA first axis')
ax2.hist(L_transformed[:CLASS_SEPARATOR],histtype='step',color='blue') 
ax2.hist(L_transformed[CLASS_SEPARATOR:],histtype='step',color='red')  
plt.show()<jupyter_output><empty_output><jupyter_text>HW2.1.b: Discussion of PCA<jupyter_code>def plot_PCA (transformed,labels):  
    lims=[-5,5]  #  Use this if defaults don't work:   ax1d.set_xlim(lims)
    fig = plt.figure(figsize=(10,3))
    #pc1,pc2,pc3 = np.split(transformed,3,axis=1)  # crashes if n_components != 3
    pc1 = [c[0] for c in transformed]
    pc2 = [c[1] for c in transformed]
    pc3 = [c[2] for c in transformed]
    # Plot first PC as histogram
    ax1d = fig.add_subplot(131)
    ax1d.set_title('First PC')
    ax1d.hist(pc1[:CLASS_SEPARATOR],histtype='step',color='blue')  
    ax1d.hist(pc1[CLASS_SEPARATOR:],histtype='step',color='red')  
    # Plot first 2 PCs in 2D. 
    ax2d = fig.add_subplot(132)
    ax2d.set_title('First 2 PCs')
    ax2d.scatter(pc1,pc2,c=labels,cmap=mycmap)
    # Plot first 3 PCs in 3D.
    ax3d = fig.add_subplot(133,projection='3d')
    ax3d.set_title('First 3 PCs')
    ax3d.scatter(pc1,pc2,pc3,c=labels,cmap=mycmap)
    # Output to screen.
    plt.show()
plot_PCA(P_transformed,labels_df)
def plot_scree(model):
    pc_values = np.arange(model.n_components_) + 1
    plt.rcParams["figure.figsize"] = (10,5)
    plt.plot(pc_values, model.explained_variance_ratio_, 'ro-', linewidth=2)
    plt.title('Scree Plot')
    plt.xlabel('Principal Component')
    plt.ylabel('Proportion of Variance Explained')
    plt.grid()
    plt.show()
plot_scree(pca)
print("Principal components (count,dimensions):",pca.components_.shape)

print("Additional variance explained, per principal component (first 10):")
for pc in range(0,10):
    print("%d=%.3f "%(pc+1,pca.explained_variance_ratio_[pc]),end=" ")
print()

print("Cumulative variance explained, per principal component (first 10):")
sum = 0.0  # alternately, use np.cumsum()
for pc in range(0,10):
    sum += pca.explained_variance_ratio_[pc]
    print("%d=%.3f "%(pc+1,sum),end=" ")
print()

def print_top_ten_features(pc,cols,cmps,TOP=10):
    print("PC%d top ten features"%pc)
    n_features=len(cols)
    mylist=[]
    for i in range(0,n_features):
        name = cols[i]
        weight = cmps[pc-1][i]
        triple = (i,name,np.abs(weight))   
        mylist.append(triple)
    myary = np.array(mylist,dtype=[('index',int),('name','S30'),('weight',float)])
    sorted = np.sort(myary,order='weight')
    for i in range(1,TOP+1):
        val = sorted[-i]  # start at end and work to front
        print("  # %3i (%30s) %7.4f"%(val[0],val[1],val[2]))
print_top_ten_features(1,column_names,pca.components_)
print_top_ten_features(2,column_names,pca.components_)
print_top_ten_features(3,column_names,pca.components_)<jupyter_output>Principal components (count,dimensions): (40, 518)
Additional variance explained, per principal component (first 10):
1=0.274  2=0.114  3=0.065  4=0.059  5=0.057  6=0.045  7=0.040  8=0.038  9=0.030  10=0.027  
Cumulative variance explained, per principal component (first 10):
1=0.274  2=0.388  3=0.453  4=0.512  5=0.569  6=0.614  7=0.654  8=0.692  9=0.722  10=0.750  
PC1 top ten features
  # 127 (      b'NBack_RT_variance_PMM')  0.0766
  # 477 (            b'NBack0_FA_AMPMV')  0.0743
  # 130 (           b'NBack_RT_std_PMM')  0.0743
  #  97 (b'Flanker(Neu)_RT_variance_PMM')  0.0739
  # 144 (     b'NBack0_RT_variance_AMM')  0.0739
  # 126 (      b'NBack_RT_variance_AMM')  0.0738
  # 480 (      b'NBack0_Accuracy_AMPMV')  0.0729
  # 147 (          b'NBack0_RT_std_AMM')  0.0726
  #  82 (    b'Flanker(Con)_RT_std_PMM')  0.0724
  # 488 (              b'NBack1_FA_PMV')  0.0722
PC2 top ten features
  # 120 (          b'NBack_RT_mean_AMM')  0.1114
  # 123 (        b'NBack_RT_median_AMM')  0.1100
[...]<jupyter_text>Non-Linear Dimensionality Reduction<jupyter_code>def construct_FAC(n_clusters):
    # Feature-wise Agglomerative Clustering (can also be applied sample-wise).
    # Agglomerative means bottom-up.
    # Agglomerative hierarchical clustering (like nearest neighbor).
    # Forms a binary tree - every node has two children - stored as array of pairs.
    # Application of StandardScalar is a prerequisite.
    # Default n_clusters is 2.
    # By default, "ward" objective function minimizes variance within each cluster.
    # transform() returns same num samples withnum features = n_clusters.
    # Transformed feature values are "pooled" for that cluster (default=np.mean()).
    # This seems to be an unsupervised decision tree.
    model = FeatureAgglomeration(n_clusters)
    # model.n_clusters_ (matches given)
    # model.labels_ (for each feature, the cluster to which it is assigned)
    # fa.children_ (array representation of the feature tree)
    # In single-linkage mode (not default), this is similar to Kruskal's algorithm!
    return model
fac = construct_FAC(3)
fac.fit(scaledMeans)
F_transformed = fac.transform(scaledMeans)
print("Example of one data instance in the reduced feature space.")
F_transformed[0]
plot_PCA(F_transformed,labels_df)<jupyter_output><empty_output><jupyter_text>The separation is not very good even at 3 levels.For binary classification (R vs NR), we would have to rely on the first level.  Random Projections<jupyter_code>def construct_GRP(n_components):
    model = GaussianRandomProjection(n_components)
    return model
grp = construct_GRP(3)
grp.fit(scaledMeans)
G_transformed = grp.transform(scaledMeans)
print("Example of one data instance in the reduced feature space.")
G_transformed[0]
plot_PCA(G_transformed,labels_df)<jupyter_output><empty_output><jupyter_text>t-SNE<jupyter_code>def construct_TSNE(n_components):
    model = TSNE(n_components)
    return model
tsne = construct_TSNE(3)
T_transformed = tsne.fit_transform(scaledMeans)  # doesn't have separate fit() and transform()
print("Example of one data instance in the reduced feature space.")
T_transformed[0]
plot_PCA(T_transformed,labels_df)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10627356873500171, 'fraction_numerical': 0.030442235173123073, 'mean_word_length': 4.607843137254902, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '307620', 'license': 'MIT', 'n_tokens_mistral': 5201, 'n_tokens_neox': 4851, 'n_words': 1395, 'path': 'PR677/Opioid_HW2_averaged.ipynb', 'repo_name': 'ShepherdCode/ShepherdML', 'size': '14585'}