<jupyter_start><jupyter_text>Logistic regression 1. Download data fromhttps://science.sciencemag.org/content/359/6378/926 (supplementary materials).* read the abstract of the article to get familiar with data origin* open the data in excel and get familiar with its content* load the protein level data (you need to figure out which one is that) as a pandas dataframe* handle missing values and convert features to numeric values when it is needed* get rid of the unnecessary (which does not encode protein levels or the tumor type) columns and the CancerSEEK results<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import pandas as pd


import sys 
import os
import conda

import seaborn as sns

#toc
from nbindex import floating
floating.tableofcontent()
#the needed data is Table S6
#there are 2 rows that we dont need
data0 = pd.read_excel('data/aar3247_cohen_sm_tables-s1-s11.xlsx', sheet_name='Table S6', skiprows=2)
display(data0.head())<jupyter_output><empty_output><jupyter_text>These stars! They use it to flag some discrepancy, or for different usage. Do they cause problem?<jupyter_code>data0.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 1821 entries, 0 to 1820
Data columns (total 45 columns):
 #   Column                                Non-Null Count  Dtype  
---  ------                                --------------  -----  
 0   Patient ID #                          1820 non-null   object 
 1   Sample ID #                           1817 non-null   object 
 2   Tumor type                            1817 non-null   object 
 3   AJCC Stage                            1005 non-null   object 
 4   AFP (pg/ml)                           1817 non-null   object 
 5   Angiopoietin-2 (pg/ml)                1817 non-null   object 
 6   AXL (pg/ml)                           1811 non-null   float64
 7   CA-125 (U/ml)                         1817 non-null   object 
 8   CA 15-3 (U/ml)                        1817 non-null   object 
 9   CA19-9 (U/ml)                         1817 non-null   object 
 10  CD44 (ng/ml)                          1811 non-null   object 
 11  CEA (pg/ml)      [...]<jupyter_text>Almost all columns have one or two starred numbers!<jupyter_code>import regex as re
#https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.replace.html
data0 = data0[:-4].replace(to_replace=r'[\*,\*\*]', value='', regex=True)
display(data0)
#facepalm, yes, this that simple... now 'to_numeric'
data0.info()
for column in data0.columns[4:-1]:
    data0[column] = pd.to_numeric(data0[column])
data0.info()
for i in data0.columns:
    print("{0} NaNs in \t{1}".format(sum(data0[i].isna()),i))
nyani = pd.unique(data0['AJCC Stage'])[3]
data0['AJCC Stage'] = data0['AJCC Stage'].replace('I','1').replace('II','2').replace('III','3').replace(nyani,'0')
pd.unique(data0['AJCC Stage'])
data0 = data0.dropna()
data0.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 1803 entries, 0 to 1816
Data columns (total 45 columns):
 #   Column                                Non-Null Count  Dtype  
---  ------                                --------------  -----  
 0   Patient ID #                          1803 non-null   object 
 1   Sample ID #                           1803 non-null   object 
 2   Tumor type                            1803 non-null   object 
 3   AJCC Stage                            1803 non-null   object 
 4   AFP (pg/ml)                           1803 non-null   float64
 5   Angiopoietin-2 (pg/ml)                1803 non-null   float64
 6   AXL (pg/ml)                           1803 non-null   float64
 7   CA-125 (U/ml)                         1803 non-null   float64
 8   CA 15-3 (U/ml)                        1803 non-null   float64
 9   CA19-9 (U/ml)                         1803 non-null   float64
 10  CD44 (ng/ml)                          1803 non-null   float64
 11  CEA (pg/ml)      [...]<jupyter_text>VERI NICE, WE CLEANED IT!<jupyter_code>cancerSeeker_outcome = data0[data0.columns[-2:]]
display(cancerSeeker_outcome) 
tumor_type = data0[data0.columns[2]]
display(tumor_type)
data0 = data0.drop(columns = ['CancerSEEK Logistic Regression Score','CancerSEEK Test Result'])
data0.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 1803 entries, 0 to 1816
Data columns (total 43 columns):
 #   Column                       Non-Null Count  Dtype  
---  ------                       --------------  -----  
 0   Patient ID #                 1803 non-null   object 
 1   Sample ID #                  1803 non-null   object 
 2   Tumor type                   1803 non-null   object 
 3   AJCC Stage                   1803 non-null   object 
 4   AFP (pg/ml)                  1803 non-null   float64
 5   Angiopoietin-2 (pg/ml)       1803 non-null   float64
 6   AXL (pg/ml)                  1803 non-null   float64
 7   CA-125 (U/ml)                1803 non-null   float64
 8   CA 15-3 (U/ml)               1803 non-null   float64
 9   CA19-9 (U/ml)                1803 non-null   float64
 10  CD44 (ng/ml)                 1803 non-null   float64
 11  CEA (pg/ml)                  1803 non-null   float64
 12  CYFRA 21-1 (pg/ml)           1803 non-null   float64
 13  DKK1 (ng/ml)      [...]<jupyter_text>2. Predict if a sample is cancerous or not* you need to build a classifier that predicts the probability of a sample coming from a cancerous (tumor type is normal or not) person based on the measured protein levels* train a logistic regression (sklearn API) on every second sample (not first 50% of the data (!), use every second line)* generate prediction for the samples that were not used during the trainingI didn't find any kind of Logistic Classifier, so I took it as only Logistic Regression is needed to be done. I imported others before reafng through the exercise description.<jupyter_code>#https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
from sklearn.neighbors import KNeighborsRegressor, KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
data0[data0.columns[2]].value_counts()<jupyter_output><empty_output><jupyter_text>If normal -> 0Everything else is cancerous -> 1<jupyter_code>data0[data0.columns[2]] = data0[data0.columns[2]].replace(['Normal','Colorectum','Breast','Lung', 'Pancreas', 'Stomach', 'Ovary', 'Esophagus', 'Liver'], 
                                ['0','1','1','1','1','1','1','1','1'])
data0[data0.columns[2]] = pd.to_numeric(data0[data0.columns[2]])
data0[data0.columns[3]] = pd.to_numeric(data0[data0.columns[3]])
data0[data0.columns[2]].value_counts()  #strange one-hot encoding hahaha!
data0.info()
#this is for testing if i know it right!
#lista = np.array([10,20,30,40])
#display(lista[1::2])

#lets separate stuff
X_ = data0[data0.columns[4:]]   #X
Y_ = data0[data0.columns[2]]     #target Y

X_train_2 = X_[1::2]
X_test_2 = X_[::2]

Y_train_2 = Y_[1::2]
Y_test_2 = Y_[::2]

#‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’, ‘saga’
LogResPred = LogisticRegression(max_iter=100000,solver='lbfgs').fit(X_train_2,Y_train_2)
Y_pred_2 = LogResPred.predict(X_test_2)
print("My prediction is {0:.3f}% accurate".format(LogResPred.score(X_test_2, Y_test_2)*100))<jupyter_output>My prediction is 83.481% accurate<jupyter_text>3. Comparision to CancerSEEK* plot the ROC curve and calculate the confusion matrix for the predictions* do the same for the CancerSEEK predictions* compare your model's performance to CancerSEEK performance<jupyter_code>from sklearn import metrics
from sklearn.metrics import roc_curve as roc_curve
#Predict confidence scores for samples.
LogResPred.decision_function(X_test_2)
fpr, tpr, thresholds = roc_curve(Y_test_2, LogResPred.decision_function(X_test_2))
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.plot(fpr,tpr)

ax.set_title("Receiver Operating Characteristic \n My cancer prediction score is {0:.3f}%".format(LogResPred.score(X_test_2, Y_test_2)*100), fontsize=25)
ax.set_xlabel("False positive rates", fontsize=22)
ax.set_ylabel("True positive rates", fontsize=22)
ax.grid()
ax.tick_params(axis="both", labelsize=15)

fig.tight_layout()

plt.show()
cancerSeeker_outcome[cancerSeeker_outcome.columns[1]] = pd.to_numeric(cancerSeeker_outcome[cancerSeeker_outcome.columns[1]].replace(['Positive', 'Negative'], ['1', '0']))
cancerSeeker_outcome
def accurace_score(arr1, arr2):
    """
    INPUT:
        2 arrays to compare
    OUTPUT:
        Percentage of similarity by counting equal values
    NOTE:
        1. Doesn't check for same lenght
        2. Doesn't check for non-zero lenght
    """
    return sum(np.equal(arr1,arr2))/len(arr1)
fpr_2, tpr_2, thresholds = roc_curve(Y_test_2, cancerSeeker_outcome[cancerSeeker_outcome.columns[0]][0::2])

fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.plot(fpr_2,tpr_2)

ax.set_title("Receiver Operating Characteristic \n CancerSEEK {0:.3f}%".format(accurace_score(np.array(Y_test_2),np.array(cancerSeeker_outcome[cancerSeeker_outcome.columns[1]][0::2]))*100), fontsize=25)
ax.set_xlabel("False positive rates", fontsize=22)
ax.set_ylabel("True positive rates", fontsize=22)
ax.grid()
ax.tick_params(axis="both", labelsize=15)

fig.tight_layout()

plt.show()
fig,ax = plt.subplots(1,1,figsize=(8,8))

ax.plot(fpr,tpr, label="My prediction")
ax.plot(fpr_2,tpr_2, label="cancerSEEK")
ax.plot(np.linspace(-1.1,1.1,100),np.linspace(-1.1,1.1,100),"--c")

ax.set_title("Receiver Operating Characteristic \n Both".format(accurace_score(np.array(Y_test_2),np.array(cancerSeeker_outcome[cancerSeeker_outcome.columns[1]][0::2]))*100), fontsize=25)
ax.set_xlabel("False positive rates", fontsize=22)
ax.set_ylabel("True positive rates", fontsize=22)
ax.grid()
ax.tick_params(axis="both", labelsize=15)

ax.set_xlim(-0.1,1.1)
ax.set_ylim(-0.1,1.1)

ax.legend(loc="best", fontsize=16)

fig.tight_layout()

plt.show()<jupyter_output><empty_output><jupyter_text>Confusion matric<jupyter_code>C_MY = metrics.confusion_matrix(Y_test_2, Y_pred_2)
C_CK = metrics.confusion_matrix(Y_test_2, np.array(cancerSeeker_outcome[cancerSeeker_outcome.columns[1]][0::2]))

display(C_MY)
display(C_CK)
fig,ax = plt.subplots(1,2,figsize=(14,7))

for i in range(2):
    ax[i].set_xticks(np.arange(0,len(C_MY[0,:])))            #class ticks place
    ax[i].set_yticks(np.arange(0,len(C_MY[:,0])))            #survived ticks place
    ax_x = ["Normal", "Cancer"]                           #class tick labels
    ax_y = ["Normal", "Cancer"]                           #survived tick labels
    ax[i].set_xticklabels(ax_x, fontsize=18)
    ax[i].set_yticklabels(ax_y, fontsize=18)

    
ax[0].set_title("My predictions's Confusion Matrix", fontsize=20)
ax[1].set_title("CancerSEEK's Confusion Matrix", fontsize=20)
ax[0].imshow(C_MY, cmap="gray")
ax[1].imshow(C_CK, cmap="gray")

for i in range(0,len(C_MY)):
    for j in range(0,2):
        color0="black"
        if(int(round(C_MY[i,j])) < int(round(np.max(C_MY)))*0.7):
            color0 = "white"
        text = ax[0].text(j, i, str(C_MY[i, j]),
                       ha="center", va="center", color=color0, fontsize=21)    
        
for i in range(0,len(C_CK)):
    for j in range(0,2):
        color0="black"
        if(int(round(C_CK[i,j])) < int(round(np.max(C_CK)))*0.7):
            color0 = "white"
        text = ax[1].text(j, i, str(C_CK[i, j]),
                       ha="center", va="center", color=color0, fontsize=21)    

fig.tight_layout()

plt.show()<jupyter_output><empty_output><jupyter_text>So, my prediction said that someone is normal rather then cancerous compared to cancerSEEK, which called normal people cancerous, while accurately predicting normal people.`CONCLUSION`: cancerSEEK is missing a lot of cancerous people in the process, which is not safe. It may be more wasteful but safer if we call more people cancerous, rather then missing a lot.  4. Hepatocellular carcinoma* fit a logistic regression (using statsmodels API this time) to predict if a sample has Hepatocellular carcinoma (liver cancer) or not. You need to keep only the liver and the normal samples for this exercise! For fitting use only the first 25 features and all the rows (which are liver or normal)* select the 5 best predictor based on P values.* Write down the most important features (based on P value) and compare them to the tumor markers that you find on wikipeida https://en.wikipedia.org/wiki/Hepatocellular_carcinoma or other sources!<jupyter_code>import statsmodels.api as sm
#ar0 = np.array([1,2,3,4,5,3,2,1])*10
#m1 = ar0 == 30
#m2 = ar0 == 20
#ar0[m1+m2]
m1 = tumor_type == 'Liver'
m2 = tumor_type == 'Normal'

t_coded = tumor_type[m1 + m2].replace(['Liver', 'Normal'], ['1', '0'])
t_coded<jupyter_output><empty_output><jupyter_text>Liver is for cancer, so it stays 1. Otherwise its normal which is 0 with my conventions! Now, let's slice!<jupyter_code>#lets separate stuff
X_2_ = data0[data0.columns[4:(25+4)]][m1+m2]   #X
Y_2_ = pd.to_numeric(t_coded)                  #target Y
 
display(X_2_.info())

#X_2_ = np.array(X_2_)
#Y_2_ = np.array(Y_2_)
#fitting with statsmodel.api LOGIT

log_reg = sm.Logit(Y_2_, X_2_).fit()
display(log_reg.summary())
highest_p = np.sort(log_reg.pvalues)[::-1][0:5]
display(highest_p)
locs = data0[data0.columns[4:(25+4)]][m1+m2].columns == "nyani"
for i in range(5):
    mask = highest_p[i] == log_reg.pvalues 
    locs = locs + mask
for i in range(5):
    print(data0[data0.columns[4:(25+4)]][m1+m2].columns[locs][i])<jupyter_output><empty_output><jupyter_text>These are the mos impactful predictors. Let's check if they align with the experiments results!`CHECKS SITE`It seems like AFP should be in the 5 ostimpactful feature. Maybe scaling data solves it?<jupyter_code>from sklearn.preprocessing import StandardScaler

def norm_data(x):
    #I DONT LIKE BURNING
    x = pd.DataFrame(StandardScaler().fit_transform(x), index=x.index, columns=x.columns)
    return x
X_2_ = norm_data(data0[data0.columns[4:(25+4)]][m1+m2])   #normed X
Y_2_ = pd.to_numeric(t_coded)                             #target Y

log_reg = sm.Logit(Y_2_, X_2_).fit()

display(log_reg.summary())

highest_p = np.sort(log_reg.pvalues)[::-1][0:5]
display(highest_p)
locs = data0[data0.columns[4:(25+4)]][m1+m2].columns == "nyani"
for i in range(5):
    mask = highest_p[i] == log_reg.pvalues 
    locs = locs + mask
for i in range(5):
    print(data0[data0.columns[4:(25+4)]][m1+m2].columns[locs][i])<jupyter_output>Optimization terminated successfully.
         Current function value: 0.090309
         Iterations 12<jupyter_text>Now they seem to be better as Leptin is connected to type-2 diabetes, which can be a common metabollic risk factor. Midkine is a great hit too as it seem to be connected to the development of the tumor. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5464806/On the other hand, G-CSF seems to be not conencted at all. [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6498295/]And all others seem to b connected to a different type of cancer...  5. Multiclass classification* Again, using every second datapoint train a logistic regression (sklearn API) to predict the tumor type. It is a multiclass classification problem.* Generate prediction for the rest of the dataset and show the confution matrix for the predictions!* Plot the ROC curves for the different cancer types on the same plot!* Intepret your results. Which cancer type can be predicted the most reliably?<jupyter_code>#tumor type
display(tumor_type.value_counts())
t_type_2 = tumor_type.replace(['Normal','Colorectum','Breast','Lung', 'Pancreas', 'Stomach', 'Ovary', 'Esophagus', 'Liver'], 
                                ['0','1','2','3','4','5','6','7','8'])
display(t_type_2.value_counts())
#lets separate stuff
X_ = data0[data0.columns[4:]]    #X
Y_ = t_type_2     #target Y

X_train_5 = X_[1::2]
X_test_5 = X_[::2]

Y_train_5 = Y_[1::2]
Y_test_5 = Y_[::2]

#‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’, ‘saga’
LogResPred = LogisticRegression(max_iter=10000000,solver='saga').fit(X_train_5,Y_train_5)
Y_pred_5 = LogResPred.predict(X_test_5)
print("My multiclass prediction has an accuracy of {0:.3f}% ".format(LogResPred.score(X_test_5, Y_test_5)*100))
C_MUL = metrics.confusion_matrix(Y_test_5, Y_pred_5)

fig,ax = plt.subplots(1,1,figsize=(14,14))

for i in range(0,len(C_MUL)):
    for j in range(0,len(C_MUL)):
        color0="black"
        if(int(round(C_MUL[i,j])) < int(round(np.max(C_MUL)))*0.7):
            color0 = "white"
        text = ax.text(j, i, str(C_MUL[i, j]),
                       ha="center", va="center", color=color0, fontsize=21)    


ax.set_xticks(np.arange(0,len(C_MUL[0,:])))            #class ticks place
ax.set_yticks(np.arange(0,len(C_MUL[:,0])))            #survived ticks place
ax_x = ['Normal','Colorectum','Breast','Lung', 'Pancreas', 'Stomach', 'Ovary', 'Esophagus', 'Liver']                           #class tick labels
ax_y = ['Normal','Colorectum','Breast','Lung', 'Pancreas', 'Stomach', 'Ovary', 'Esophagus', 'Liver']                           #survived tick labels
ax.set_xticklabels(ax_x, fontsize=18)
ax.set_yticklabels(ax_y, fontsize=18)

    
ax.set_title("Multiclass' Confusion Matrix", fontsize=20)
ax.imshow(C_MUL, cmap="gray")

fig.tight_layout()

plt.show()
#'key of type tuple not found and not a MultiIndex'
one_hot_encoded = np.array(pd.get_dummies(Y_test_5))

labels = ['Normal','Colorectum','Breast','Lung', 'Pancreas', 'Stomach', 'Ovary', 'Esophagus', 'Liver'] 

fig,ax=plt.subplots(1,1,figsize=(10,10))

for i in range(9):
    fpr, tpr, thresholds = roc_curve(one_hot_encoded[:,i], LogResPred.decision_function(X_test_5)[:,i])    
    ax.plot(fpr,tpr, label=labels[i])
    
    
ax.plot(np.linspace(-0.1,1.1,100),np.linspace(-0.1,1.1,100),"--c")


ax.set_title("Multiclass ROC", fontsize=26)
ax.set_xlabel("False positive rates", fontsize=22)
ax.set_ylabel("True positive rates", fontsize=22)
ax.grid()
ax.tick_params(axis="both", labelsize=15)
ax.legend(loc="best", fontsize=22)

ax.set_xlim(-0.1,1.1)
ax.set_ylim(-0.1,1.1)

fig.tight_layout()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13466126615990678, 'fraction_numerical': 0.046995505742662154, 'mean_word_length': 3.431767887878043, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 49, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '218692', 'license': 'MIT', 'n_tokens_mistral': 6888, 'n_tokens_neox': 6295, 'n_words': 1559, 'path': 'hw5/HW_05.ipynb', 'repo_name': 'AdamGTaylor/DataMining_2021', 'size': '18023'}