{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Titanic Submission Using XGBoost With Model Diagnostic\n\nThis is my submission for the widely known Titanic survival classification dataset. The project will be divided into several steps starting from:\n\n- Data Import\n- Data Cleaning & Exploration\n- Creating Baseline Model\n- Hyperparameter Tuning\n- Prediction","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport re\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 01 Data Import\n\nThis section we will be working with data import and a quick glance into our data structure.\n\nBefore we import the data into our workspace, let's have a general overview of the data dictionary provided by Kaggle.","metadata":{}},{"cell_type":"markdown","source":"| Variable | Definition | Key\n| --- | --- | ---\n|survival  |Survival |0 = No, 1 = Yes\n|pclass | Ticket class | 1 = 1st, 2 = 2nd, 3 = 3rd\n|sex |Sex |\n|Age | Age in years |\n|sibsp | # of siblings / spouses aboard the Titanic\n|parch | of parents / children aboard the Titanic \n|ticket | Ticket number \n|fare | Passenger fare\n|cabin | Cabin number\n|embarked | Port of Embarkation | C = Cherbourg, Q = Queenstown, S = Southampton","metadata":{}},{"cell_type":"code","source":"# Import train and test dataset\ntrain = pd.read_csv('/kaggle/input/titanic/train.csv')\ntest = pd.read_csv('/kaggle/input/titanic/test.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Observing the Data Structure","metadata":{}},{"cell_type":"code","source":"train.info()\nprint(f\" Our training dataset consist of {train.shape[0]} training examples across {train.shape[1]} initial training features\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Seems like we have a mix between numerical, categorical, and string data. Some of our data contain missing values with `Cabin` containing the highest number of missing values.","metadata":{}},{"cell_type":"markdown","source":"### 02 Data Cleaning and Exploration\n\nThis section we will be working toward understanding our data, feature by feature, as well as exploring relationship between variables. Along the way we will also be cleaning our data of missing values, redundant columns, and perform feature engineering where we think is necessary.","metadata":{"tags":[]}},{"cell_type":"markdown","source":"#### Observe Distribution of Null Values","metadata":{"tags":[]}},{"cell_type":"code","source":"# Create series for count of null and percentage of null\ntrain_null = train.isnull().sum().sort_values(ascending=False)\ntrain_null_pct =(train_null/train.shape[0]).apply(lambda x: '{:.2%}'.format(x))\nd_type = train.dtypes\n\n# Combine null values into dataframe\nnull_df = pd.DataFrame({'null_count':train_null, \n                        'null_percentage': train_null_pct,\n                        'd_type':d_type})\n\n# Observe only columns with null values\nnull_df[null_df.null_count>0]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Only 3 columns have missing values. We can probably impute missing values on `Age` and `Embarked` but `Cabin` might not be a very useful feature consider over 77% of the data is missing. We should be exploring the distribution of these 3 columns to identify best way for imputation. Let's observe the distribution of our columns starting with `Age`.","metadata":{}},{"cell_type":"code","source":"# Create figure and subplots\nfig, axes = plt.subplots(1, 3, figsize=(20, 5))\n# Plot distribution of Age\naxes[0].set_title('Distribution of Age')\nsns.histplot(data=train, x='Age',kde=True, ax=axes[0])\n\n# Plot Cabin Distribution\naxes[1].set_title('Distribution of Cabin')\ng=sns.countplot(data=train, x=\"Cabin\", ax=axes[1])\ng.set_xticks(np.arange(0, sns.countplot(data=train, x=\"Cabin\", ax=axes[1]).get_xticks()[-1], 15))\n\n# Plot Embarked Distribution\naxes[2].set_title('Distribution of Embarked')\nsns.countplot(data=train, x=\"Embarked\", ax=axes[2], palette='Set2')","metadata":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"`Age` have a fairly normal distribution, we can impute the missing values by the mean value for now. `Cabin` have too many missing values, we can try to extract the room initial as a separate feature but very unlikely that the feature will tell us something meaningful. `Embarked` has only 2 missing values, we will impute those with majority class 'S'.","metadata":{}},{"cell_type":"markdown","source":"#### Impute Missing Values","metadata":{}},{"cell_type":"code","source":"# Impute missing values\ntrain['Age'] = train['Age'].fillna(value = np.mean(train['Age']))\ntrain['Embarked'] = train['Embarked'].fillna(value='S')\n\n# Drop cabin column\ntrain.drop('Cabin',axis=1,inplace=True)\n\n# Check data structure\ntrain.info()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Observe Distribution of Categorical Features","metadata":{}},{"cell_type":"code","source":"# Parse categorical features into categorical dtype\ntrain = train.assign(Survived = train['Survived'].astype('category'),\n             Pclass = train['Pclass'].astype('category'),\n             Sex = train['Sex'].astype('category'),\n             Embarked = train['Embarked'].astype('category'))\n\n# Create new dataframe that only includes categorical variables\ntrain_cat = train.select_dtypes(include = 'category')\n\n# Check data structure\ntrain_cat.info()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create figure and subplots\nfig, axes = plt.subplots(2, 2, figsize=(16, 7))\n\n# Loop through axes and plot each column\nfor i, ax in enumerate(fig.axes):\n    ax.set_title(f\"Distribution of {train_cat.columns[i]}\")\n    sns.countplot(data=train, x=train_cat.columns[i], order = train_cat.iloc[:,i].value_counts().index, ax=ax, palette='Set2')\n\nfig.tight_layout()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Initial observation among the categorical variables show that:\n\n1. Less than half of passengers survived.\n2. Most passengers have 3rd class tickets with 2nd class tickets having the lowest frequency.\n3. Most passengers are male.\n4. Majority of the passengers embarked on Southampton.","metadata":{}},{"cell_type":"markdown","source":"#### Observation of Numerical Variables","metadata":{}},{"cell_type":"code","source":"# Create new dataframe that only includes numerical variables\ntrain_num = train.select_dtypes(include = np.number)\ntrain_num.drop('PassengerId', axis = 1, inplace=True)\n\n# Check data structure\ntrain_num.info()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create figure and subplots\nfig, axes = plt.subplots(2, 2, figsize=(16, 7))\n\n# Loop through axes and plot each column\nfor i, ax in enumerate(fig.axes):\n    ax.set_title(f\"Distribution of {train_num.columns[i]}\")\n    sns.histplot(data=train, x=train_num.columns[i], kde=True, ax=ax)\n\nfig.tight_layout()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Observing the numerical distribution aboves, these observations can be made:\n\n1. Age have a fairly normal distribution, even more so after using mean imputation on missing values.\n2. The other variables have a rather heavy right-skewed distribution.\n3. Both `SibSp` and `Parch` can be categorical variable with simpler categorization (e.g. binary class of low values against high values), but we will be keeping the data as is for now.","metadata":{}},{"cell_type":"code","source":"# Correlation between numerical values\nplt.figure(figsize=(10,10))\nsns.heatmap(train_num.corr(),cbar=True,annot=True,cmap='viridis')\nplt.title('Correlation Between Numerical Variables')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Correlation between variables are rather weak. The one that are somewhat correlated is `Parch` and `SibSp`, both representing the size of the family in different aspect.","metadata":{}},{"cell_type":"markdown","source":"#### What factors likely to influence survival rate the most?","metadata":{}},{"cell_type":"code","source":"# Crosstabluation of categorical variables against `Survived`\ndisplay(np.round(pd.crosstab(train['Survived'], train['Sex'], normalize = 'columns')*100,2))\ndisplay(np.round(pd.crosstab(train['Survived'], train['Pclass'], normalize = 'columns')*100,2))\ndisplay(np.round(pd.crosstab(train['Survived'], train['Embarked'], normalize = 'columns')*100,2))\n\n# Observe Embarked against Pclass\ndisplay(np.round(pd.crosstab(train['Pclass'], train['Embarked'], normalize = 'columns')*100,2))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Female passengers have a drastically higher survival rate compared to male passengers, probably due to lady and children first policy.\n- Ticket class also have an impact on survival rate with passengers from the upper class have a higher chance of survival.\n- Where passenger Embarked have to do with their ticket class as well. Majority of passenger from Queenstown and Southamptom are 3rd class passengers which result in lower survival rate compare to Cherbourg where higher proportion of 1st class passengers are observed.","metadata":{}},{"cell_type":"code","source":"# Create figure and subplots\nfig, axes = plt.subplots(2, 2, figsize=(16, 8), dpi = 320)\n\n# Loop through axes and plot each column\nfor i, ax in enumerate(fig.axes):\n    ax.set_title(f\"Survived X {train_num.columns[i]}\")\n    sns.boxplot(data=train, y=train_num.columns[i], x = 'Survived', orient='v',ax=ax, palette='Set2')\n\nfig.tight_layout()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Impact of `Parch`, `SibSp`, and `Age` on `Survived` is unclear.\n- Those who survived tends to pay higher fare on average than those who don't.","metadata":{}},{"cell_type":"markdown","source":"#### Creating title feature from Name variable","metadata":{}},{"cell_type":"markdown","source":"The `Name` variable, by itself is not very useful feature for ML model to use. However, the column also contain potentially useful information about the passenger social status, hence we will be creating a new feature based on their title(s).","metadata":{}},{"cell_type":"code","source":"# Split name by space\nname_list = train.Name.apply(lambda x: x.split())\n\n# Extract title\ntitle_match = []\nfor name in name_list:\n    for word in name:\n        loc = word.find('.')\n        if loc != -1:\n            title_match.append(word)\n            break\n\n# Observe distribution\npd.Series(title_match).shape\n\n# Add new col to train df\ntrain['Title'] = title_match\ntrain.Title.value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that there's still too many level of titles, many of which belong to one specific passenger, we can further group them up to aid our model picking up patterns from the feature.","metadata":{}},{"cell_type":"code","source":"# Create function to group title\ndef group_title(x):\n    # Create list of titles to be grouped\n    default_list = ['Mr.','Miss.','Mrs.','Ms.']\n    military_list = ['Major.','Col.']\n    \n    # Group titles\n    if x in default_list:\n        if x == 'Ms.':\n            x = 'Miss.'\n        else:\n            x = x\n    elif x in military_list:\n        x = 'Military'\n    else:\n        x = 'Others'\n        \n    return x\n\n# Update the `Title` column with new grouping\ntrain['Title'] = train['Title'].apply(lambda x:group_title(x)).astype('category')\n\n# Observe distribution\ntrain['Title'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Observe the relationship between the new titles with survival rate\ndisplay(np.round(pd.crosstab(train['Survived'], train['Title'], normalize = 'columns')*100,2))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems that titles do have relatively big impact on survival rate, most notable of which is the high survival rate among married women and nobles. Majority of males did not survived as inferred from the `Sex` variable.","metadata":{}},{"cell_type":"markdown","source":"#### Creating family related variables","metadata":{}},{"cell_type":"code","source":"train['FamilySize'] = train['Parch'] + train['SibSp'] + 1\ntrain['AdultFlg'] = train['Age'].apply(lambda x:1 if x >= 21 else 0).astype('category')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Drop redundant columns","metadata":{}},{"cell_type":"code","source":"train.drop(['Ticket','Name','PassengerId'], axis=1, inplace=True)\ntrain.info()","metadata":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 03 Creating Baseline Model\n\nThis section we will be creating some baseline models to be used as the benchmark for model prediction power on passenger survival chance. The models that we will be testing are:\n\n- Naive Bayes (baseline)\n- Logistic Regression (reference)\n- XGBoost","metadata":{}},{"cell_type":"code","source":"# Get dummies for non-ensemble models\nlog_x = train.drop('Survived', axis=1)\nlog_y = train.Survived\n\n# Get dummies for ensemble models\ntree_train = pd.get_dummies(train, columns=train.select_dtypes(include='category').columns)\ntree_x = tree_train.drop(['Survived_1','Survived_0'], axis = 1)\ntree_y = tree_train.Survived_1","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import cross_validate\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.preprocessing import StandardScaler, OneHotEncoder\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.compose import ColumnTransformer\n\n# Scale numerical features and apply one hot encoding to categorical\nnumeric_features = log_x.select_dtypes(np.number).columns\nnumeric_transformer = StandardScaler()\ncategorical_features = log_x.select_dtypes('category').columns\ncategorical_transformer = OneHotEncoder(drop='first')\n\n# Create preprocessor steps to process columns\npreprocessor = ColumnTransformer(transformers = [(\"num\", numeric_transformer, numeric_features), (\"cat\", categorical_transformer, categorical_features)])\n\n# Create pipeline for preprocessing and fitting logistic regression model\nclf = Pipeline(steps=[(\"preprocessor\", preprocessor), (\"classifier\", LogisticRegression(penalty='none', max_iter=500))])\nscores = cross_validate(clf, log_x, log_y, scoring = ['accuracy','f1'], cv=5)\n\n# Print performance\nscores = pd.DataFrame(scores)\nprint(f\" Average {scores.columns[-2]} is {scores.iloc[:,-2].mean():.4f}, Average {scores.columns[-1]} is {scores.iloc[:,-1].mean():.4f}\")\nscores","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.naive_bayes import GaussianNB\n\ngnb_model = GaussianNB()\nscores = cross_validate(gnb_model, tree_x, tree_y, scoring = ['accuracy','f1'], cv=5)\nscores = pd.DataFrame(scores)\nprint(f\" Average {scores.columns[-2]} is {scores.iloc[:,-2].mean():.4f}, Average {scores.columns[-1]} is {scores.iloc[:,-1].mean():.4f}\")\nscores","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from xgboost import XGBClassifier\nxgb = XGBClassifier(random_state = 123)\n\n# Create pipeline for preprocessing and fitting logistic regression model\nclf = Pipeline(steps=[(\"preprocessor\", preprocessor), (\"classifier\", xgb)])\nscores = cross_validate(clf, log_x, log_y, scoring = ['accuracy','f1'], cv=5)\n\n# Print performance\nscores = pd.DataFrame(scores)\nprint(f\" Average {scores.columns[-2]} is {scores.iloc[:,-2].mean():.4f}, Average {scores.columns[-1]} is {scores.iloc[:,-1].mean():.4f}\")\nscores","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Comparing across 3 models with base parameters we have:\n\n- 79.1% CV Accuracy on Naive Bayes\n- 81.7% CV Accuracy for Logistic Regression\n- 82.0% CV Accuracy for XGBoost.","metadata":{}},{"cell_type":"markdown","source":"### 04 Hyperparameter Tuning","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import validation_curve\n\n# Define function to create validation score for specific parameter\ndef validation_score(model, x, y, param_name, param_range, cv):\n    train_scores, valid_scores = validation_curve(model, x, y, param_name =param_name, param_range = param_range, n_jobs=-1,cv=cv)\n    train_mean_param = np.mean(train_scores,axis=1)\n    valid_mean_param = np.mean(valid_scores,axis=1)\n    train_score = pd.DataFrame({param_name :  param_range, 'score' : train_mean_param, 'label' : 'train'})\n    valid_score = pd.DataFrame({param_name :  param_range, 'score' : valid_mean_param, 'label' : 'valid'})\n    score = pd.concat([train_score, valid_score], ignore_index=True)\n    return score","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Metric in consideration, all of which is for adjusting bias/variance\nmetrics_eval = ['max_depth','lambda','colsample_bytree','gamma', 'subsample','eta']\nmetrics_param_range = [np.arange(1,25,1), np.arange(0, 100,2), np.linspace(0,1,10), np.arange(0,25,1), np.linspace(0,1,5), np.linspace(0,1,5)]\nfig, axes = plt.subplots(2, 3, figsize=(16, 6), dpi = 320)\n\nfor i, ax in enumerate(fig.axes):\n    ax.set_title(f\"{metrics_eval[i]}\")\n    score=validation_score(xgb,tree_x,tree_y,metrics_eval[i],metrics_param_range[i],5)\n    sns.lineplot(data=score, x=metrics_eval[i], y=\"score\", hue = 'label',ax=ax)\n\nfig.tight_layout()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Observing the distribution above suggests that our model might slightly overfit the training data as seen with the initial increase in validation performance when we increase the `lambda` parameter. Other metrics such as `max_depth` which allows for a more complex model and `min_child_weight` that controls model complexity do not show big improvement in validation score as we adjust it forward.","metadata":{}},{"cell_type":"code","source":"#from sklearn.model_selection import GridSearchCV \n\n#xgb = XGBClassifier(random_state = 123)\n\n# Grid search\n#param_grid = {\n    #'max_depth': np.arange(1,10,2),\n    #'lambda': np.linspace(0,40,10),\n    #'colsample_bytree':np.linspace(0,1,5),\n    #'gamma': np.arange(0,10,2),\n    #'subsample': np.linspace(0,1,5)\n#}\n\n#gs_xgb = GridSearchCV(xgb, param_grid = param_grid, verbose = 1, n_jobs=-1, scoring = 'accuracy')\n#xgb_model = gs_xgb.fit(tree_x, tree_y)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import GridSearchCV \n\nxgb = XGBClassifier(random_state = 123)\n\n# Grid search\nparam_grid = {\n    'max_depth': [7],\n    'lambda': [4.444444444444445],\n    'colsample_bytree':[0.25],\n    'gamma':[0],\n    'subsample': [0.5]\n}\n\ngs_xgb = GridSearchCV(xgb, param_grid = param_grid, verbose = 1, n_jobs=-1, scoring = 'accuracy')\nxgb_model = gs_xgb.fit(tree_x, tree_y)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Best performance seen is: {xgb_model.best_score_}\")\nprint(f\"Parameters that achieved best score are: {xgb_model.best_params_}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from xgboost import plot_importance\nplot_importance(xgb_model.best_estimator_, importance_type = 'gain')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 05 Model Prediction","metadata":{}},{"cell_type":"code","source":"x_test = test.assign(Pclass = test['Pclass'].astype('category'),\n             Sex = test['Sex'].astype('category'),\n             Embarked = test['Embarked'].astype('category'))\n\n# Split name by space\nname_list = x_test.Name.apply(lambda x: x.split())\n# Extract title\ntitle_match = []\nfor name in name_list:\n    for word in name:\n        loc = word.find('.')\n        if loc != -1:\n            title_match.append(word)\n            break\n# Observe distribution\npd.Series(title_match).shape\n# Add new col to train df\nx_test['Title'] = title_match\nx_test['Title'] = x_test['Title'].apply(lambda x:group_title(x)).astype('category')\nx_test.drop(['Cabin','Ticket','Name','PassengerId'], axis=1, inplace=True)\nx_test['FamilySize'] = x_test['Parch'] + train['SibSp'] + 1\nx_test['AdultFlg'] = x_test['Age'].apply(lambda x:1 if x >= 21 else 0).astype('category')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_test = pd.get_dummies(x_test, columns=x_test.select_dtypes(include='category').columns)\nx_test.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Survived = xgb_model.predict(x_test)\nfinal_output = pd.DataFrame({'PassengerId' : test.PassengerId, 'Survived' : Survived})\nfinal_output.to_csv('final_submission_2.csv',index=False)","metadata":{},"execution_count":null,"outputs":[]}]}