<jupyter_start><jupyter_text>Covid-19 diagnosis using symptoms[<NAME>](https://www.pgupta.info) 2020-05-20**Abstract:** This notebook is a tutorial on building a typical machine learning classifier. The process of doing so spans steps ranging from data wrangling to model selection. With the help of each of these steps we hope to make the reader familiar with challenges involved in building a machine learning system.  QuestionIn this hackathon, we want to build a machine learning model to predict COVID-19 infections from symptoms.It has several applications, for example, triaging patients to be attended to by a doctor or nurse, recommending self-isolation through contact tracing apps. Zoabi et al. [[1]](https://www.nature.com/articles/s41746-020-00372-6) builds a decision tree classifier using the publicly available data reported by the Israeli Ministry of Health.The paper itself dicsusses the various challenges encountered in deploying such a model. It is encouraged to read the paper and learn the challeges and ways to overcome them. However, in this hackathon, we will use their dataset and make the participant familiar with a typical pipeline of building a machine learning system.[1] [<NAME>., <NAME>. & <NAME>. Machine learning-based prediction of COVID-19 diagnosis based on symptoms. npj Digit. Med. 4, 3 (2021).]((https://www.nature.com/articles/s41746-020-00372-6))  Setup the workspaceWe will clone their Git repository to to use their dataset<jupyter_code>import pandas as pd
import math
import numpy as np

from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import LabelBinarizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import plot_roc_curve, roc_auc_score, roc_curve

import matplotlib.pyplot as plt
from matplotlib.lines import Line2D<jupyter_output><empty_output><jupyter_text>Data Let's check how the data looks like and how various features are encoded.<jupyter_code>df = pd.read_csv('data\corona_tested_individuals_ver_006.english.csv.zip')
df.info() # No nulls
df.head()
# We need to one hot encode corona_result, age_60_and_above, gender and test_indication<jupyter_output><empty_output><jupyter_text>Great! These are the features used in the paper for their prediction task. The authors also list these features in the [README.md of their Github repo](https://github.com/nshomron/covidpred).  Pandas read columns as string, so we need to convert them to the proper format before we can operate on it.<jupyter_code>df['test_date'] = pd.to_datetime(df['test_date'])
print("Start date:", min(df['test_date']))
print("End date:", max(df['test_date']))

n_days = (max(df['test_date']) - min(df['test_date'])).days
print("# of days: ", n_days)

n_obs =  df.shape[0]
print("# of observations:", df.shape[0])
print("# of features:", df.shape[1])

pos_cases = sum(df['corona_result'] == "positive")
print("# of positively diagnosed cases: {0} ({1: 2.2f}%)".format(pos_cases, 100*pos_cases / n_obs))

neg_cases = sum(df['corona_result'] == "negative")
print("# of negatively diagnosed cases: {0} ({1: 2.2f}%)".format(neg_cases, 100 * neg_cases / n_obs))

other_cases = sum(df['corona_result'] == "other") # possibly not confirmed
print("# of other cases (possibly, not confirmed): {0} ({1: 2.2f}%)".format(other_cases, 100 * other_cases / n_obs))

# Small percentage of positively diagnosed cases - we need to be careful about our accuracy metric. ROC AUC is best.<jupyter_output>Start date: 2020-03-11 00:00:00
End date: 2020-04-30 00:00:00
# of days:  50
# of observations: 278848
# of features: 10
# of positively diagnosed cases: 14729 ( 5.28%)
# of negatively diagnosed cases: 260227 ( 93.32%)
# of other cases (possibly, not confirmed): 3892 ( 1.40%)<jupyter_text>Since we do not have any information on what happened to "other" cases, we will exclude them from our exercise.<jupyter_code>df = df[df['corona_result'].isin(['positive', 'negative'])]<jupyter_output><empty_output><jupyter_text>To build out predictor, we will be splitting our dataset into **training,  validation, and test sets**. A model is trained on the training dataset while the hyperparameters are tuned on the validation dataset. Finally, a test dataset is used to report final model's performance metrics. Since we have a time dependent dataset we will split our training and test dataset based on time. Thus, we find the date before which 60% of observations are present, and use that as our training dataset.We will use next 20% of the dataset as our validation dataset, and finally, the remaining 20% will be used as a test dataset. Thus, we use 60-20-20 split.The authors use 63%-23% training-test split, and a further split of training into train-valid dataset using 80-20% split.There is no prescribed formula on how to do this split.<jupyter_code>date_counts = df.groupby(['test_date']).count()['gender'] # take count of any column. They will all be same.
date_counts = date_counts.sort_index()
cum_counts = date_counts.cumsum()
cdf = cum_counts / n_obs

max_training_date = cdf[cdf < 0.60].index.max()
training_data = df[df['test_date'] <= max_training_date]

min_test_date = cdf[cdf > 0.80].index.min()
test_data = df[df['test_date'] >= min_test_date]

valid_data = df[(max_training_date < df['test_date']) & (df['test_date'] < min_test_date)]

print("# of observations in training dataset", training_data.shape[0])
print("# of observations in validation dataset", valid_data.shape[0])
print("# of observations in test dataset", test_data.shape[0])<jupyter_output># of observations in training dataset 160717
# of observations in validation dataset 53172
# of observations in test dataset 61067<jupyter_text>**NOTE:**It is **extremely important** that you do not use the test dataset in the model building phase. While building models, it is required to tune the hyperparameter, adjust assumptions, modify features, etc. This should be done on the validation dataset. After several such iterations on the validation dataset, you will pick a model with the best performace as your final model. A test dataset is used to measure the final model's performance, which is a proxy for how it will perform (or generalize) in real life. Thus, to have a proper measure of model's genearalization, test dataset should not be part of your model building process.   Exploratory Data AnalysisIn this section, we will see the general statistics of features.In doing so, we will encounter inconsistencies in the data and address them accordingly.<jupyter_code># We want to predict 'corona_result'. 
# We will not use "test_date" as a feature. 
# So we narrow down the input features to this list 
INPUT_FEATURES = ['cough', 'fever', 'sore_throat', 'shortness_of_breath', 'head_ache', 'age_60_and_above', 'gender', 'test_indication']
TARGET_COLUMN = 'corona_result'

for col in INPUT_FEATURES:
  print("*"*25, f" {col} ", "*"*25)
  print(training_data[col].value_counts())<jupyter_output>*************************  cough  *************************
0       119702
1        25547
0        12039
1         3177
None       252
Name: cough, dtype: int64
*************************  fever  *************************
0       131874
0        13768
1        13375
1         1448
None       252
Name: fever, dtype: int64
*************************  sore_throat  *************************
0       142874
0        15973
1         1435
1          434
None         1
Name: sore_throat, dtype: int64
*************************  shortness_of_breath  *************************
0       143279
0        15906
1         1030
1          501
None         1
Name: shortness_of_breath, dtype: int64
*************************  head_ache  *************************
0       142294
0        16064
1         2015
1          343
None         1
Name: head_ache, dtype: int64
*************************  age_60_and_above  *************************
No      123867
Yes      25425
None     11425
Name: age_60_and_above, dtype: [...]<jupyter_text>Here are the questions to guide you through the process of exploring data 1.   Think about possible biases and limitations of this dataset2.   What is the format of feature values? Are there any inconsistencies? If so, how would you make them cosistent?3.   What is the statistics of these feature values? How many symptoms are reported or not?4.   Which symptoms have a reporting bias, i.e., likely to be reported when the patient is COVID positive? 5.   How will the symptoms with reporting bias affect the model’s performance?6.   Visualization: Draw the bar graph of features grouped by the target class? 7.   How does the bar graph of the symptoms with reporting bias looks like?8.   Determine if we have a class imbalance in the dataset? If so, what do you reckon will be the downstream challenges in evaluating the model? How will you overcome those challenges?9.   What does "None" value mean for feature? Should we include these features? **GOOD PRACTICE**: To make your plots accessible to everyone, it is always a good idea to use colorblind-friendly palette for your plots. Check out [this](https://medium.com/cafe-pixo/inclusive-color-palettes-for-the-web-bbfe8cf2410e) for such a palette.  Feature EngineeringIn this section, we will transform the features that models can operate upon. Note that this transformation doesn't have to be unique. It is very much dependent on the type of model you are building. Here is the list of questions to guide your feature engineering task 1.   How will you represent the features in numerical format that can be accessible by model? 2.   Are there any redundancies in your feature representation?3.   How will you represent targets in a format accessible to the model?Check out [`sklearn`'s preprocessing library](https://scikit-learn.org/stable/modules/preprocessing.html) for easy-to-use functions to do this.<jupyter_code>def normalize(x):
  """
  Normalizes the input to integer type. It maps "None" to 2.
  Args:
    x (int or str): input to be normalized 
  Returns:
    (int): normalized input
  """
  if type(x) == str:
    return int(x) if x != "None" else 2
  return x

def normalize_columns(data):
  for col in ['cough', 'fever', 'sore_throat', 'shortness_of_breath', 'head_ache']:
      data[col] = data[col].apply(normalize)
  return data
from sklearn.preprocessing import OneHotEncoder, LabelBinarizer
pd.options.mode.chained_assignment = None  # default='warn'


def preprocess(data, encoder, lb):
    """
    Transforms `data` into format required for model building
    Args:
      data (pd.DataFrame): dataframe with columns `INPUT_FEATURES` and `TARGET_COLUMN`
      encoder (sklearn.preprocessing.OneHotEncoder): A fitted OneHotEncoder to be used to transform `INPUT_FEATURES`
      lb (sklearn.preprocessing.LabelBinarizer): A fitted LabelBinarizer to be used to transform `TARGET_COLUMN`
    Returns:
      model_input (np.array): each row is an observation, columns are one-hot encoded features of `INPUT_FEATURES`
      model_target (np.array): 1D array with 1 where `TARGET_COLUMN` is "positive" and 0 otherwise.
    """

    # Clean values, turning all values into ints to consolidate duplicate values.
    normalize_columns(data)

    # One hot encode features.
    model_input = encoder.transform(data[INPUT_FEATURES]).todense()
    model_target = lb.transform(data[TARGET_COLUMN]).flatten()

    return model_input, model_target


# Clean values, turning all values into ints to consolidate duplicate values.
normalize_columns(training_data)

lb = LabelBinarizer()
lb.fit(training_data['corona_result'])
encoder = OneHotEncoder(drop='first').fit(training_data[INPUT_FEATURES]) 

X_train, y_train = preprocess(training_data, encoder, lb)
X_valid, y_valid = preprocess(valid_data, encoder, lb)
X_test, y_test = preprocess(test_data, encoder, lb)<jupyter_output><empty_output><jupyter_text>Model buildingIn this section, we will build various classifiers using `sklearn`. You do not have to restrict yourself to `sklearn`. Please feel free to use any other library.**TRY:**  Try various classifiers that you have learned so far.Here is the list of models to try :*  Logistic Regssion: [User Guide](https://scikit-learn.org/stable/modules/linear_model.htmllogistic-regression). [API](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html).*   Decision Trees: [User Guide](https://scikit-learn.org/stable/modules/tree.htmlclassification). [API](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.htmlsklearn.tree.DecisionTreeClassifier). You can head down in the User Guide to [other Tree algorithms](https://scikit-learn.org/stable/modules/tree.htmltree-algorithms-id3-c4-5-c5-0-and-cart) if you fancy. *   Categorical Naive Bayes.[User Guide](https://scikit-learn.org/stable/modules/naive_bayes.htmlcategorical-naive-bayes). [API](https://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.CategoricalNB.htmlsklearn.naive_bayes.CategoricalNB)*   Linear Discriminant Analysis. [User Guide](https://scikit-learn.org/stable/modules/lda_qda.html). [API](https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.LinearDiscriminantAnalysis.htmlsklearn.discriminant_analysis.LinearDiscriminantAnalysis)*   Quadratic Discriminant Analysis. [User Guide](https://scikit-learn.org/stable/modules/lda_qda.html). [API](https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis.htmlsklearn.discriminant_analysis.QuadraticDiscriminantAnalysis)*   Support Vector Machines. [User Guide](https://scikit-learn.org/stable/modules/svm.htmlclassification). [API](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.htmlsklearn.svm.SVC)*   Nearest neighbors classification. [User Guide](https://scikit-learn.org/stable/modules/neighbors.htmlnearest-neighbors-classification). [API](https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsClassifier.htmlsklearn.neighbors.KNeighborsClassifier)*   Neural networks - Multi-layer Perceptron (MLP). [User Guide](https://scikit-learn.org/stable/modules/neural_networks_supervised.htmlmulti-layer-perceptron). [API](https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.htmlsklearn.neural_network.MLPClassifier)<jupyter_code>from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier


# Training baseline models with no hyperparameter tweaking.

lr = LogisticRegression()
lr.fit(X_train, y_train)
dtc = DecisionTreeClassifier()
dtc.fit(X_train, y_train)
from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier()
rf.fit(X_train, y_train)
from sklearn.neighbors import KNeighborsClassifier
#knn = KNeighborsClassifier(n_neighbors=19, weights = 'uniform')
#knn_model = knn.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Evaluate ModelIn this section, we will evaluate our model's performance on the validation dataset. Here are the list of questions to think about while deciding how to evaluate your model - *   Is accuracy the right metric to evaluate the model? Are inaccuracies correctly penalized in the accuracy metric?*  Would you think that the cost of false negative is more than the false positive? Is it dependent on the application?*  Which metric will minimize false negatives and false positives?*   Which dataset should you chose to evaluate the model? Validation or Test?What other metric is relevant in our context?  For benchmarking everyone’s results we will stick to ROC AUC score as a metric. There are standard functions to compute these scores in `sklearn`, so we will use them. Specifically, we will be using [`roc_auc_score`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.htmlsklearn.metrics.roc_auc_score) and [`plot_roc_curve`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.plot_roc_curve.html).<jupyter_code>def evaluate(model, X, Y):
  """
  Returns the AUC-ROC for `model` as evaluated on (X, Y)

  Args:
    model (): Any model that has a function predict_proba and returns probability for each row in `X`.
    X (np.array): Input to the model containing feature values
    Y (np.array): 1D array containing true class i.e. 0 or 1
  
  Returns:
    (float): AUC ROC for the model
  """

  y_score = model.predict_proba(X) # (n_samples, n_clases) with each value being the probability of being in that class
  return roc_auc_score(y_true=Y, y_score=y_score[:, 1])

print('')
print(f'LR train score: {evaluate(lr, X_train, y_train)}')
print(f'LR valid score: {evaluate(lr, X_valid, y_valid)}')

print(f'DTC train score: {evaluate(dtc, X_train, y_train)}')
print(f'DTC valid score: {evaluate(dtc, X_valid, y_valid)}')

print(f'RF train score: {evaluate(rf, X_train, y_train)}')
print(f'RF valid score: {evaluate(rf, X_valid, y_valid)}')

print(f'knn train score: {evaluate(knn, X_train, y_train)}')
print(f'knn valid score: {evaluate(knn, X_valid, y_valid)}')<jupyter_output>LR train score: 0.8990602419946635
LR valid score: 0.82602918308811
DTC train score: 0.9197129777640092
DTC valid score: 0.8166051160298665
RF train score: 0.9196518046448063
RF valid score: 0.8209756694371307<jupyter_text>Hyperparameter SearchIn this section we will be searching for the best parameters to build our models. This is where we will use our validation dataset.  Hyperparameter search can become messy if you have lots of paramters. A brute force method to do such a search will be to do a grid search to fit tons of models. Thus, a smarter way to do hyperparameter search has been the subject of research. **TRY:** If interested, read [here](https://scikit-learn.org/stable/modules/grid_search.html) for more details and incorporate some of those ideas in the model building process.<jupyter_code>from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold

# Grid search for Logistic Regression
lr_params = {
    'solver': ['newton-cg', 'lbfgs', 'liblinear'],
    'penalty': ['l2'],
    'C': [100, 10, 1.0, 0.1, 0.01]
}
cv = RepeatedStratifiedKFold(n_splits=10, n_repeats=3, random_state=1)
lr_grid_search = GridSearchCV(estimator=lr, param_grid=lr_params, n_jobs=-1, cv=cv, scoring='accuracy', error_score=0)
grid_result = lr_grid_search.fit(X_train, y_train)
print(grid_result.best_score_)
print(grid_result.best_params_)
best_lr = grid_result.best_estimator_
from numpy import mean
from numpy import std
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedStratifiedKFold
from sklearn.svm import SVC
from sklearn.ensemble import VotingClassifier
from matplotlib import pyplot
from sklearn.ensemble import StackingClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier

models = list()
models.append(('lr', LogisticRegression(C= 1.0, penalty='l2', solver= 'newton-cg'))) 
models.append(('rf', RandomForestClassifier(max_depth= 5, min_samples_split= 2, n_estimators= 20)))
models.append(('knn', KNeighborsClassifier(n_neighbors=19, weights = 'uniform')))
	# define the voting ensemble
ensemble = StackingClassifier(estimators=models)
ensemble.fit(X_train, y_train)
ensemble.predict(X_test)
roc_auc_score(y_test, y_score)<jupyter_output><empty_output><jupyter_text>Report your results**NOTE:** You should use the test dataset only when you are done with hyperparameter search on your model. This is because the test dataset is not involved in the model building process, thereby making sure that the performance evaluation on the test dataset measures how well the proposed model is able to generalize.<jupyter_code>X_test, Y_test = preprocess(test_data, encoder, lb)
auc_score = evaluate(best_lr, X_test, Y_test)
print(f"LR AUC-ROC on the test dataset:{auc_score}")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1043245653143112, 'fraction_numerical': 0.02516470996185664, 'mean_word_length': 4.9887273805992285, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 47, 'https://': 26, 'lorem ipsum': 0, 'www.': 3, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '473528', 'license': 'MIT', 'n_tokens_mistral': 6428, 'n_tokens_neox': 5929, 'n_words': 2092, 'path': 'Previous versions/FINAL_COVID_diagnosis_hackathon.ipynb', 'repo_name': 'nomopix/cs-ai-l7-oct-hackathon', 'size': '20207'}