<jupyter_start><jupyter_text>Drug Allergy Prediction Author   <NAME>   <EMAIL>      https://thanaphon.dev     BackgroundAccording to the World Health Organization, Drug Allergy is unintended and harmful reactions from drug intake under normal conditions. Drug Allergy presents a serious challenge in healthcare systems as it can increase healthcare costs and deaths.A laboratory technique called Enzyme-linked immunospot (ELISpot) can test whether the patient's immune cells will respond to particular drugs. This allows a doctor to screen for drug safety for specific patients. Unfortunately, the test does not always guarantee that the patients will not be allergic to the medicine due to the occasional lack of response in ELISpot test and vice versa.There are several works published regarding drug allergy prediction. [<NAME>., <NAME>., <NAME>. et al., 2018](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2544-0) leverages Convolutional Neural Network to predict adverse drug reactions and identify associated molecular substructures. The model performs very well with over 90% accuracy and closes to 85% AUC of the ROC curve via the "Neural Fingerprint" method. On the other hand, [<NAME>, 2018](https://www.ijcmph.com/index.php/ijcmph/article/view/2471/1845) takes the traditional modeling approach and create a model which predict hospitalizations and deaths from drug allergies using supervised learning models including Logistical regression, SVM, and XGradient Boost. The model performs from 65-75% accuracy with the F1 score in the 80s. These studies have access to millions of labeled data points from United States healthcare statistics. Real-world data such as biological and medical datasets often suffer small samples with very few labeled datasets. This presents a challenge in creating an effective machine learning model. Technologies such as active learning and semi-supervised learning have been proposed to utilize the unlabeled samples for improving the model performance. [<NAME>., <NAME>., Wang, SS. et al., 2018](https://www.researchgate.net/publication/327287262_A_novel_logistic_regression_model_combining_semi-supervised_learning_and_active_learning_for_disease_classification) proposes a novel approach to this problem by using a logistic regression model with high confidence to propagate labels to the unlabeled dataset.  Research QuestionsIn this task, the goal is to develop a prediction model for drug allergy based on patient and drug information - a binary classification problem. Due to the small dataset and the small dimensionality of the input features, the experiment will leverage simpler traditional machine learning models over complicated deep neural networks. Semi-supervised learning techniques will also be used to accommodate a large number of unlabeled datasets. The research questions of this experiment can be summarized as follow:1 - Can semi-supervised learning techniques such as K-prototype clustering or the novel logistical regression model presented by [<NAME>., <NAME>., Wang, SS. et al., 2018](https://www.researchgate.net/publication/327287262_A_novel_logistic_regression_model_combining_semi-supervised_learning_and_active_learning_for_disease_classification) be used to effectively propagate the training labels?2 - Can semi-supervised learning techniques improve the final model given a large number of unlabeled samples?3 - Do machine learning models agree with the assertation of strong predictive property in ELISpot variables to drug allergy?4 - Given a lack of domain expertise to perform any feature engineering, will model with built-in feature selections such as tree-based models enhance model performance?5 - What is the best classification model for this task and its performance? Can it perform better than a dummy classifier that simply guesses the answer based on prior probability?  Performance measurementsLabeled ground truth will be split for final evaluation while preserving the posterior probability for each class. The model will be assessed based on the macro f-1 score and the area under the curve (AUC) of the receiver operating characteristic (ROC).   0 Environment Setup (Co-lab) 0.1 Install any missing packages<jupyter_code>!pip install kmodes<jupyter_output>Requirement already satisfied: kmodes in /usr/local/lib/python3.6/dist-packages (0.10.2)
Requirement already satisfied: scikit-learn>=0.19.0 in /usr/local/lib/python3.6/dist-packages (from kmodes) (0.22.2.post1)
Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.6/dist-packages (from kmodes) (0.14.1)
Requirement already satisfied: numpy>=1.10.4 in /usr/local/lib/python3.6/dist-packages (from kmodes) (1.18.4)
Requirement already satisfied: scipy>=0.13.3 in /usr/local/lib/python3.6/dist-packages (from kmodes) (1.4.1)<jupyter_text>0.2 ImportssIt's typical to include package imports at the top of the file for clarity.<jupyter_code>import os
import numpy as np 
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from kmodes.kprototypes import KPrototypes
from sklearn.cluster import KMeans
from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_transformer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import RobustScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import StratifiedKFold
from sklearn.svm import LinearSVC
import xgboost as xgb
from sklearn.metrics import f1_score, roc_curve, auc
from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression,LogisticRegressionCV
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix, log_loss, precision_recall_fscore_support
from sklearn.model_selection import cross_val_score
from sklearn.metrics import log_loss
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import SVC
from statistics import stdev, mean

%matplotlib inline

#Generate random seed for reproducability
seed = 0
rng = np.random.RandomState(seed)


# Ignore warnings for final delivery
import warnings
warnings.filterwarnings('ignore')<jupyter_output>/usr/local/lib/python3.6/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm<jupyter_text>0.3 Helper functions<jupyter_code>def plot_confusion_matrix(cm, classes=None, title='Confusion matrix'):
    """Plots a confusion matrix."""
    if classes is not None:
        sns.heatmap(cm, xticklabels=classes, yticklabels=classes, vmin=0., vmax=1., annot=True)
    else:
        sns.heatmap(cm, vmin=0., vmax=1., annot=True)
    plt.title(title)
    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    
def plot_ROC(y, y_pred):
    fpr, tpr, threshold = roc_curve(y, y_pred)
    roc_auc = auc(fpr, tpr)
    plt.title('Receiver Operating Characteristic')
    plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
    plt.legend(loc = 'lower right')
    plt.plot([0, 1], [0, 1],'r--')
    plt.ylabel('True Positive Rate')
    plt.xlabel('False Positive Rate')
    plt.show()<jupyter_output><empty_output><jupyter_text>0.4 Mount the drive for data<jupyter_code>from google.colab import drive
drive.mount('/content/drive')<jupyter_output>Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).<jupyter_text>1. Data Exploration 1.1 Quick look at the dataset<jupyter_code>raw_dataset = pd.read_csv(os.path.join(os.getcwd(), "drive/My Drive/Colab Notebooks/datasets", "DrugAllergyKaggle_v3_050720.csv"), delimiter = ',')
print("There are {} datapoints and {} attributes in the raw dataset.".format(raw_dataset.shape[0], raw_dataset.shape[1]))
raw_dataset.head(10)
print("Ground Truth Labels")
raw_dataset["Drug_Rechallenge_Result"].value_counts()
plt.figure(figsize=(5,5))
sns.countplot(x = raw_dataset['Drug_Rechallenge_Result'])
plt.xticks([0,1],['No', 'Yes'])
plt.xlabel('Drug Allergy')
plt.title('Ground Truth Data')
plt.show()
print("Training Label")
raw_dataset["Allergic_Reaction_Group"].value_counts()
raw_dataset.describe()
raw_dataset.info()
raw_dataset.hist(figsize=(20,15))
plt.title("Count plot of all features")
plt.show()<jupyter_output><empty_output><jupyter_text>1.2 Explaination of each variables ELISpot_Control The ELISpot test result for the POSITIVE CONTROL (i.e., we expect to see strong response) ELISpot_Result The ELISpot test result for SUSPECTED DRUG (i.e., this is the result that indicate whether the patient would be allergic to that drug) NARANJO_Category An ordinal data. Exposure_Time The amount of times since the patient has taken the drug until the ELISpot test date. Suspicion_ScoreThe suspicion level of the drug (1 = suspected drug, 2 = similar to suspected drug, 3 = negative control). This is ORDINAL. Allergic_Reaction_Group The severity of patient's allergic reaction. This is ORDINAL. Drug_Group A categorical data. Drug_Rechallenge_ResultThe ground truth of this dataset that we want to predict.  2. Data Analysis 2.1 Categorical, Ordinal, and Continous VariablesIt is important to separate these variables as some machine learning algorithms require categorical data to be encoded and numerical data to be standardized. Binary categorical variables are inherently encoded.<jupyter_code>continous_variable_columns = ['Age_Year', 'ELISpot_Control', 'ELISpot_Result', 'Naranjo_Score', 'Exposure_Time']
ordinal_variable_columns =  ['Naranjo_Category', 'Suspicion_Score', 'Allergic_Reaction_Group', 'Suspicion_Score']
binary_variable_columns = ['Gender', 'Steroid_Usage', 'Underlying_Condition_A', 'Underlying_Condition_B', 'Underlying_Condition_C', 'Underlying_Condition_D', 'Underlying_Condition_E']
categorical_variable_columns = ['Drug_Group']<jupyter_output><empty_output><jupyter_text>2.2 Target LabelDrug_Rechallenge_Result is the ground truth of this dataset. This is the dependent variable the model will try to predict.<jupyter_code>target_column = 'Drug_Rechallenge_Result'<jupyter_output><empty_output><jupyter_text>2.3 Patient_IDPatient_ID is a database primary key. It is to be discarded as it serves no statistical significance.<jupyter_code>raw_dataset = raw_dataset.drop(['Patient_ID'], 1)<jupyter_output><empty_output><jupyter_text>2.4 Association Analysis and Dimensionality ReductionAlthough this is a binary classification task, it may still be useful to visualize the strength of the association between the continuous variables and to gain some insight about a particular feature. Highly correlated independent variables could be compressed further using dimensionality reduction techniques such as PCAs to improve training performance.<jupyter_code>#Split data into batches for easier viewing. There are 15 columns

for row in continous_variable_columns:
    columns = []
    for column in continous_variable_columns:
        if row != column:
             columns.append(column)    
    sns.pairplot(data=raw_dataset,y_vars=[row],x_vars=columns, height=3)    
plt.show()<jupyter_output><empty_output><jupyter_text>2.4.1  Scatterplots AnalysisIt does not appear to have any strong linear relationships among the independent continuous variables. Therefore, PCA will not be an attempt in this experiment.  2.5 ELISpot_Control and OutliersELISpot_Control seems to follow a normal distribution but there appear to be some bad outliers based on  IQR, std, min, max, and mean. This outlier will be dealt with using the empirical rule (3 sigma rule).Outliers can introduce bias into the model during the training phase, but the effect of outliers may still be useful during the inference phase. Therefore, eliminating outliers during model development may be preferred over minimizing the outlier effect through log transformation.<jupyter_code>index = np.abs(raw_dataset['ELISpot_Control'] - raw_dataset['ELISpot_Control'].mean()) / raw_dataset['ELISpot_Control'].std() <= 3
raw_dataset = raw_dataset[index]
raw_dataset.shape
raw_dataset["ELISpot_Control"].hist(figsize=(5,5))
plt.title("ELISpot_Control after removing outliers")
plt.show()<jupyter_output><empty_output><jupyter_text>3. Data Preprocessing 3.1 Labeled and Unlabeled DataSplit labeled and unlabeled data.<jupyter_code># We can use NaN != NaN to split labeled data
unlabeled_data = raw_dataset.query("Drug_Rechallenge_Result != Drug_Rechallenge_Result")
labeled_data = raw_dataset.query("Drug_Rechallenge_Result == Drug_Rechallenge_Result")
print("There are {} labeled data points and {} unlabeled data points.".format(labeled_data.shape[0], unlabeled_data.shape[0]))<jupyter_output>There are 113 labeled data points and 678 unlabeled data points.<jupyter_text>Sanity ChecksLabeled data should have 1s and 0s. Unlabeled data should have all NaN.<jupyter_code>labeled_data["Drug_Rechallenge_Result"].head(5)
unlabeled_data["Drug_Rechallenge_Result"].head(5)
labeled_data_X = labeled_data.drop([target_column], 1)
labeled_data_Y = labeled_data[target_column]

unlabeled_data_X = unlabeled_data.drop([target_column], 1)
unlabeled_data_Y = unlabeled_data[target_column]<jupyter_output><empty_output><jupyter_text>3.2 Training and Validation DataBecause the dataset is partially labeled, half of the labeled data will be split into validation and training set.<jupyter_code>from sklearn.model_selection import StratifiedShuffleSplit

splitter = StratifiedShuffleSplit(n_splits=2, random_state=None, test_size=0.5)
[index, _] = splitter.split(labeled_data_X, labeled_data_Y)

train_X = labeled_data_X.iloc[index[0]]
train_Y = labeled_data_Y.iloc[index[0]]
test_X = labeled_data_X.iloc[index[1]]
test_Y = labeled_data_Y.iloc[index[1]]
print("There are {} training data points and {} validation data points".format(train_X.shape[0], test_X.shape[0]))<jupyter_output>There are 56 training data points and 57 validation data points<jupyter_text>3.2.1 Sanity CheckCheck if the split works correctly and preserves the class ratio.<jupyter_code>fig, (ax1, ax2) = plt.subplots(1, 2)
sns.countplot(x = train_Y, ax=ax1)
ax1.set_title("Train")
sns.countplot(x = test_Y, ax=ax2)
ax2.set_title("Test")
plt.show()<jupyter_output><empty_output><jupyter_text>3.3 Imputing Missing Value There are quite a bit of missing Data points in many of the input variables. This requires imputation before it can be processed by Machine Learning algorithms. There are many methods to impute data and it often depends on the context. In some circumstances, replacing missing data might be the wrong thing to do, and zero out the values might have been better. Unfortunately, due to the lack of domain expertise, this experiment will impute based on median value for continuous data and based on mode for everything else.  3.4 Standardize Continous Variable Machine Learning algorithms that use distance metrics such as Mean Square Error are generally sensitive to outliers and differences in Scale. Standardizing these features to zero mean and unit variance can lead to better performance. SKlearn's RobustScaler will be used in this task to standardize the data while minimizing the effect of outliers. 3.5 Categorical Data EncodingMany machine learning models do not accept categorical data points. They will have to be encoded using One-Hot Encoding. The exception to this is the binary categorical data, as they are inherently encoded.  3.5 Preprocessing PipelineWrapping up the preprocessing steps into a pipeline. One with one-hot encoding and the other without to accommodate certain algorithms.<jupyter_code># Transformers
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', RobustScaler())])

categorical_onehot_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent')),
     ('onehot', OneHotEncoder(handle_unknown='ignore'))])

categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent'))])

binary_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent'))])

ordinal_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent'))])

# Preprocessing Pipelines
preprocessor_with_one_hot = ColumnTransformer(
    transformers=[
        ('num', numeric_transformer, continous_variable_columns),
        ('ord', ordinal_transformer, ordinal_variable_columns),
        ('bin', binary_transformer, binary_variable_columns),
        ('cat', categorical_onehot_transformer, categorical_variable_columns)])

preprocessor = ColumnTransformer(
    transformers=[
        ('num', numeric_transformer, continous_variable_columns),
        ('ord', ordinal_transformer, ordinal_variable_columns),
        ('bin', binary_transformer, binary_variable_columns),
        ('cat', categorical_transformer, categorical_variable_columns)])


# Numerical Column Index
category_column_index = [train_X.columns.get_loc(c) for c in categorical_variable_columns if c in train_X]<jupyter_output><empty_output><jupyter_text>4. Semi-Supervised LearningDue to a large number of the unlabeled dataset, K-prototype and Logistical regression will be trained over the labeled samples then attempt to propagate the labels over the unlabeled samples.  4.1 Label Propagation using K-prototypeK-prototype is a variant of Kmean and Kmode unsupervised learning model to deal with a combination of ordinal, categorical, and numerical datasets.<jupyter_code>cluster_clf = Pipeline(steps=[('preprocessor', preprocessor),
                  ('cluster', KPrototypes(n_clusters=2, init='Huang', random_state=0))])
cluster_clf.fit(train_X, train_Y, cluster__categorical=category_column_index)
y_predict = cluster_clf.predict(train_X, categorical=category_column_index)
cm = confusion_matrix(train_Y, y_predict)
cm = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm, title="K-Prototype clustering CM")    
plt.show()<jupyter_output><empty_output><jupyter_text>4.1.1 K-prototype resultUnfortunately, based on the confusion matrix alone. K-prototype gives no meaningful clusters as it labels almost everything into one particular class. Hence, it will not be used for label propagation.  4.2 Label Propagation using Logistical RegressionBased on the method proposed by  [<NAME>., <NAME>., <NAME>. et al., 2018](https://www.researchgate.net/publication/327287262_A_novel_logistic_regression_model_combining_semi-supervised_learning_and_active_learning_for_disease_classification), a logistical regressor will be trained and samples with high confidence will see label propagation.<jupyter_code>lg_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification', LogisticRegression(random_state=0, max_iter=20000))]) 
lg_clf.fit(train_X, train_Y)
y_predict = lg_clf.predict(train_X)
cm = confusion_matrix(train_Y, y_predict)
cm = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm, title="LG CM")    
plt.show()<jupyter_output><empty_output><jupyter_text>4.2.1 Logistical regression resultThe logistical regression performs relatively well based on the Confusion Matrix.  4.2.2 Identify high confidence samplesLeverage Logistical regression to identify samples with a very high classification confidence. (Over 90% for each class).<jupyter_code>prob_Y = lg_clf.predict_proba(unlabeled_data_X)

threshold = 0.9

i = 0
for a,b in prob_Y:
    if a >= threshold:
        unlabeled_data_Y.iloc[i] = 0
    elif b >= threshold:
        unlabeled_data_Y.iloc[i] = 1
    else:
        unlabeled_data_Y.iloc[i] = -1
    i += 1

index = unlabeled_data_Y >= 0<jupyter_output><empty_output><jupyter_text>4.3 Label PropagationPropagate identified samples to the new augmented dataset.<jupyter_code>augmented_train_X = train_X.append(unlabeled_data_X[index])
augmented_train_Y = pd.Series([*train_Y, *unlabeled_data_Y[index]])
print("New augmented training set has {} data points".format( augmented_train_X.shape[0]))<jupyter_output>New augmented training set has 500 data points<jupyter_text>5. Model SelectionsIn this section, several machine learning models will be trained and evaluated against the augmented training set. Ideally, the one with a significantly higher averaged macro f1 score across all cross-validations will be selected for final evaluation.  5.0.1  Evaluate FunctionThe following function will evaluate model performance by calculating the average macro F-1 score across 5 folds. Then display a confusion matrix from one of the folds for a sanity check.<jupyter_code>def evaluate_clf(clf, clf_name):
    skf = StratifiedKFold(n_splits=5, random_state=None, shuffle=True)
    count = 0
    f1_sum = 0
    for train_index, test_index in skf.split(augmented_train_X, augmented_train_Y):
        X_tr, X_val = augmented_train_X.iloc[train_index], augmented_train_X.iloc[test_index]
        Y_tr, Y_val = augmented_train_Y.iloc[train_index], augmented_train_Y.iloc[test_index]    
        clf.fit(X_tr, Y_tr)    
        y_predict = clf.predict(X_val)
        f1_sum += f1_score(Y_val, y_predict, average='macro')
        count += 1
    # Display CM from one of the K fold
    cm = confusion_matrix(Y_val, y_predict)
    cm = cm / cm.sum(axis=1)[:, np.newaxis]
    plot_confusion_matrix(cm, title="CM for {}".format(clf_name))    
    plt.show()
    
    # Show average F-1 scores across all folds
    print("Average Macro F-1 Score: {}".format(f1_sum/count))<jupyter_output><empty_output><jupyter_text>5.1 Gaussion Naive BayesAre input features conditionally independent? Naive Bayes can give more insight into this relationship.<jupyter_code>naive_bayes_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification',GaussianNB())]) 
evaluate_clf(naive_bayes_clf, "Naive Bayes")<jupyter_output><empty_output><jupyter_text>5.1.1 DiscussionResult suggests good training performance. Input features are likely to be independent.  5.2 Random ForestRandom Forest creates boundary lines based on discriminative features and automatically perform feature selections. This can give some insights into the most discriminating feature of the dataset.<jupyter_code>random_forest_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification',RandomForestClassifier())]) 
evaluate_clf(random_forest_clf, "Random Forest")<jupyter_output><empty_output><jupyter_text>5.2.1 Random Forest - Relative Feature ImportanceRank each feature based on most discriminative features such as Information Gain, etc.<jupyter_code># Reconstruct the input column name after one hot encoding

Drug_Group = [ "Drug_Group" for c in range(1, 9)]

one_hot_encoded_features = ['Gender', 'Age_Year', 'ELISpot_Control', 'ELISpot_Result',
       'Naranjo_Score', 'Naranjo_Category', 'Exposure_Time', 'Steroid_Usage',
       'Underlying_Condition_A', 'Underlying_Condition_B',
       'Underlying_Condition_C', 'Underlying_Condition_D',
       'Underlying_Condition_E', 'Suspicion_Score', 'Allergic_Reaction_Group',
       *Drug_Group]

one_hot_encoded_importance = random_forest_clf.steps[1][1].feature_importances_


# Aggregate the importance value for each feature 
features_importances = {}  

for imp, fea in zip(one_hot_encoded_importance, one_hot_encoded_features):
    if fea not in features_importances:
        features_importances[fea] = 0
    features_importances[fea] += imp

# Display the feature importance
features = list(features_importances.keys())
importances = np.array(list(features_importances.values()))
indices = np.argsort(importances)

plt.figure(figsize=(5,5))
plt.title('Feature Importances (Random Forest)')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), [features[i] for i in indices])
plt.xlabel('Relative Importance')
plt.show()<jupyter_output><empty_output><jupyter_text>5.3 XGradientBoostXGradient Boost is a popular tree-based ensemble model that often outperforms Random Forest in Kaggle competitions.<jupyter_code>xgb_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification', xgb.XGBClassifier())]) 
evaluate_clf(xgb_clf, "XGBoost")<jupyter_output><empty_output><jupyter_text>5.3.1 XGBoost - Relative Feature ImportanceRank each feature based on the most discriminative features.<jupyter_code># Aggregate the importance value for each feature 
one_hot_encoded_importance = xgb_clf.steps[1][1].feature_importances_
features_importances = {}  

for imp, fea in zip(one_hot_encoded_importance, one_hot_encoded_features):
    if fea not in features_importances:
        features_importances[fea] = 0
    features_importances[fea] += imp

# Display the feature importance
features = list(features_importances.keys())
importances = np.array(list(features_importances.values()))
indices = np.argsort(importances)

plt.figure(figsize=(5,5))
plt.title('Feature Importances (XGBoost)')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), [features[i] for i in indices])
plt.xlabel('Relative Importance')
plt.show()<jupyter_output><empty_output><jupyter_text>5.3.1 DiscussionThe tree models performed quite well. The most discriminative features are ELISpot variables, which agree with the initial assertion.  5.4 Support Vector Machine (Linear Kernel)  Support Vector Machine generally performs well on a medium to small size dataset.<jupyter_code>linear_svc_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification', LinearSVC(penalty='l2', loss='squared_hinge',  max_iter=100000, dual = False) )]) 
evaluate_clf(linear_svc_clf, "Linear SVC")<jupyter_output><empty_output><jupyter_text>5.5 Model SelectionThe clear winner cannot be determined due to the marginal differences in the average F-1 score for each model. This score will likely have a high variance as the sample size is quite small. Despite this, the XGBoost model will be selected for final evaluation due to a slightly better average F-1 macro score in some executions. Tree-based algorithms also have much better explainability in higher-dimensional space.  6. Model OptimizationFrom the previous section, XGBoost performed the best on the training data. This step will optimize XGBoost for optimal hyperparameter using Grid search.<jupyter_code>from sklearn.model_selection import GridSearchCV
from scipy.stats import uniform


optimization_dict = {'max_depth': [2,4,6,8,10,12,14,16,18,20],
                     'n_estimators': [50,100,200,400]}

clf = GridSearchCV(xgb.XGBClassifier(), optimization_dict, scoring="f1_macro")
clf.fit(preprocessor_with_one_hot.fit_transform(augmented_train_X),augmented_train_Y)

print("Best XGBoost parameters")
print(clf.best_score_)
print(clf.best_params_)<jupyter_output>Best XGBoost parameters
0.967701754385965
{'max_depth': 2, 'n_estimators': 100}<jupyter_text>7. Model EvaluationThis section evaluates XGBoost against the Test dataset and baseline classifiers.  7.1 XGBoost Evaluation - with artificially labelled dataset<jupyter_code>precision = []
recall = []
f1 = []
for i in range(10):
    xgb_clf = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification', xgb.XGBClassifier(max_depth=2, n_estimators=100))]) 
    xgb_clf.fit(augmented_train_X, augmented_train_Y)
    y_predict = xgb_clf.predict(test_X)
    scores = precision_recall_fscore_support(test_Y, y_predict, average='macro')
    precision.append(scores[0])
    recall.append(scores[1])
    f1.append(scores[2])

# Display CM from one of tests
cm = confusion_matrix(test_Y, y_predict)
cm = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm, title="CM for XGBoost (Augmented Data)")    
plt.show()

# Plot ROC
plot_ROC(test_Y, y_predict)

# Summarize Scores
print("Precision {} std: {}".format(mean(precision), stdev(precision)))
print("Recall    {} std: {}".format(mean(recall), stdev(recall)))
print("F1        {} std: {}".format(mean(f1), stdev(f1)))<jupyter_output><empty_output><jupyter_text>7.2 XGBoost Evaluation - with only ground truth label<jupyter_code>precision = []
recall = []
f1 = []
for i in range(10):
    xgb_clf2 = Pipeline(steps=[('preprocessor', preprocessor_with_one_hot),
                  ('classification', xgb.XGBClassifier(max_depth=2, n_estimators=100))]) 
    xgb_clf2.fit(train_X, train_Y)
    y_predict = xgb_clf.predict(test_X)
    scores = precision_recall_fscore_support(test_Y, y_predict, average='macro')
    precision.append(scores[0])
    recall.append(scores[1])
    f1.append(scores[2])

# Display CM from one of tests
cm = confusion_matrix(test_Y, y_predict)
cm = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm, title="CM for XGBoost (Ground Truth only)")    
plt.show()

# Plot ROC
plot_ROC(test_Y, y_predict)

# Summarize Scores
print("Precision {} std: {}".format(mean(precision), stdev(precision)))
print("Recall    {} std: {}".format(mean(recall), stdev(recall)))
print("F1        {} std: {}".format(mean(f1), stdev(f1)))<jupyter_output><empty_output><jupyter_text>7.3 Dummy ClassifierEvaluate a classifier where the only prediction it makes is zero (No Allergy).<jupyter_code>precision = []
recall = []
f1 = []
for i in range(10):
    dummy_clf = Pipeline(steps=[('preprocessor', preprocessor),
                  ('classification', DummyClassifier(strategy='constant', constant=0))]) 
    dummy_clf.fit(train_X, train_Y)
    y_predict = dummy_clf.predict(test_X)
    scores = precision_recall_fscore_support(test_Y, y_predict, average='macro')
    precision.append(scores[0])
    recall.append(scores[1])
    f1.append(scores[2])

# Display CM from one of tests
cm = confusion_matrix(test_Y, y_predict)
cm = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm, title="CM for Dummy Classifier".format("s"))    
plt.show()

# Plot ROC
plot_ROC(test_Y, y_predict)

# Summarize Scores
print("Precision {} std: {}".format(mean(precision), stdev(precision)))
print("Recall    {} std: {}".format(mean(recall), stdev(recall)))
print("F1        {} std: {}".format(mean(f1), stdev(f1)))<jupyter_output><empty_output><jupyter_text>8. Conclusion In this experiment, I have explored the patient data and drug information in an attempt to create a predictive model that would identify vulnerable patients. Due to the lack of labeled data, a semi-supervised learning technique is deployed to increase training data size the supervised several machine learning models.Initially, a Logistical Regression model is trained from the labeled dataset, then propagate the label over to unlabeled data. Several models are evaluated afterward. Tree-based models score better average macro F-1 score in several executions with some noise and they have a better explainability compared to other models. Hence, XGBoost is selected to perform the final evaluation against the test dataset with an average F1 macro score of roughly 0.78 and AUC of  0.77. Compared to the XGBoost model using the human-labeled dataset, the augmented model sees small improvement over the model without label propagation. Dummy baseline classifier is also trained and evaluated with a 0.43 macro F1 score and 0.5 AUC, significantly underperforming against XGBoost. During the experimentation, Random Forest and XGBoost find the most discriminative feature of this dataset to be the EL5Spot variable. This agrees with the expert assertation that the EL5Spot is the best indicator of drug allergy.In conclusion, this experiment paints a promising result in using a semi-supervised learning model for predicting drug allergy in patients using only a small sample of labeled data. Nevertheless, future studies should be conducted with more human-labeled data to reaffirm this finding. Small evaluation and training dataset can introduce sampling errors with high variance.<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08547008547008547, 'fraction_numerical': 0.014887334887334888, 'mean_word_length': 5.888460715050311, 'pattern_counts': {'":': 0, '<': 132, '<?xml version=': 0, '>': 138, 'https://': 6, 'lorem ipsum': 0, 'www.': 4, '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': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '121770', 'license': 'Unlicense', 'n_tokens_mistral': 9905, 'n_tokens_neox': 9277, 'n_words': 3037, 'path': '1_Drug_Allergy_Prediction.ipynb', 'repo_name': 'As-12/Drug-Allergy-Prediction', 'size': '32223'}