w9_training_validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biol 359  |  Cross-Validation\n",
    "### Spring 2021, Week 9\n",
    "\n",
    "<hr style=\"border:2px solid gray\"> </hr>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns \n",
    "from sklearn import datasets\n",
    "from sklearn import linear_model\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load microarray dataset\n",
    "\n",
    "dataset = pd.read_csv('microarray_data_clean2.csv')\n",
    "input_data = dataset.iloc[:, 4:]\n",
    "output_data = dataset.iloc[:, 3]\n",
    "conditions_data = dataset.iloc[:,0]\n",
    "\n",
    "#input_data\n",
    "#output_data\n",
    "#conditions_data\n",
    "dataset.head(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate scatter plot of (first 25) independent vs dependent variable\n",
    "plt.style.use('ggplot')\n",
    "fig = plt.figure(figsize = (20, 20))\n",
    " \n",
    "for index, feature_name in enumerate(input_data.columns):\n",
    "    if index >= 25:\n",
    "        break\n",
    "    ax = fig.add_subplot(5, 5, index + 1)\n",
    "    ax.scatter(input_data.iloc[:, index], output_data) \n",
    "    ax.set_ylabel('viability [24hr]', size = 16)\n",
    "    ax.set_xlabel(feature_name, size = 16)\n",
    " \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define X & Y variables\n",
    "X = input_data\n",
    "Y = output_data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split data into training and validation set\n",
    "x_train, x_validation, y_train, y_validation = train_test_split(X, Y, test_size = 0.20)\n",
    " \n",
    "print(\"Shape of training data X = % s and Y = % s : \"%(x_train.shape, y_train.shape)) \n",
    "print(\"Shape of validation data X = % s and Y = % s : \"%(x_validation.shape, y_validation.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (1) Multiple Linear Regression (Ordinary Least Squares)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train multiple linear regression model\n",
    "lreg = LinearRegression()\n",
    "lreg.fit(x_train, y_train)\n",
    "\n",
    "# Predict validation data\n",
    "lreg_y_pred = lreg.predict(x_validation)\n",
    "\n",
    "print(\"Linear regression R2 (training data):\",f\"{lreg.score(x_train, y_train):.3f}\")\n",
    "print(\"Linear regression Q2 (validation data):\",f\"{lreg.score(x_validation, y_validation):.3f}\" \"\\n\")\n",
    "\n",
    "sum_squared_error = round(np.sum((lreg.predict(x_train) - y_train)**2),3)\n",
    "print(\"SSE on training data: \", sum_squared_error)\n",
    " \n",
    "sum_squared_error = round(np.sum((lreg.predict(x_validation) - y_validation)**2),3)\n",
    "print(\"SSE on validation data: \", sum_squared_error)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parity plot\n",
    "\n",
    "plt.plot(lreg_y_pred, y_validation, '.', label='validation')\n",
    "plt.plot(lreg.predict(x_train), y_train, '.', label='training')\n",
    "\n",
    "min_value = min(min(lreg_y_pred), min(y_validation), min(lreg.predict(x_train)), min(y_train))\n",
    "max_value = max(max(lreg_y_pred), max(y_validation), max(lreg.predict(x_train)), max(y_train))\n",
    "plt.plot([min_value, max_value],[min_value, max_value], '--', label=\"parity\")\n",
    "plt.xlabel('Predicted Values')\n",
    "plt.ylabel('True Values')\n",
    "ax.set_box_aspect(1)\n",
    "sns.despine()\n",
    "plt.title('Parity Plot')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine coefficient and corresponding variable names\n",
    "lreg_coefficient = pd.DataFrame()\n",
    "lreg_coefficient[\"Columns\"] = x_train.columns\n",
    "lreg_coefficient['Coefficient Estimate'] = pd.Series(lreg.coef_)\n",
    "#lreg_coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the magnitude of coefficients\n",
    "fig, ax = plt.subplots(figsize =(7, 7))\n",
    "color =['midnightblue', 'navy', 'darkblue', 'mediumblue', 'blue', \n",
    "        'darkslateblue','mediumslateblue', 'slateblue', \n",
    "        'mediumpurple', 'blueviolet', 'darkviolet']\n",
    "\n",
    "ax.bar(lreg_coefficient[\"Columns\"], lreg_coefficient['Coefficient Estimate'], color = color)\n",
    "ax.spines['bottom'].set_position('zero')\n",
    " \n",
    "plt.style.use('ggplot')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (2) Ridge Regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train a Ridge regression model\n",
    "from sklearn.linear_model import Ridge\n",
    " \n",
    "ridgeR = Ridge(alpha = 0.2)\n",
    "ridgeR.fit(x_train, y_train)\n",
    "\n",
    "# Predict validation data\n",
    "ridge_y_pred = ridgeR.predict(x_validation)\n",
    " \n",
    "print(\"Ridge regression R2 (training data):\",f\"{ridgeR.score(x_train, y_train):.3f}\")\n",
    "print(\"Ridge regression Q2 (validation data):\",f\"{ridgeR.score(x_validation, y_validation):.3f}\" \"\\n\")\n",
    "\n",
    "sum_squared_error = round(np.sum((ridgeR.predict(x_train) - y_train)**2),3)\n",
    "print(\"SSE on training data: \", sum_squared_error)\n",
    " \n",
    "sum_squared_error = round(np.sum((ridgeR.predict(x_validation) - y_validation)**2),3)\n",
    "print(\"SSE on validation data: \", sum_squared_error)\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parity plot\n",
    "\n",
    "plt.plot(ridgeR.predict(x_validation), y_validation, '.', label='validation')\n",
    "plt.plot(ridgeR.predict(x_train), y_train, '.', label='train')\n",
    "\n",
    "min_value = min(min(ridgeR.predict(x_validation)), min(y_validation), min(ridgeR.predict(x_train)), min(y_train))\n",
    "max_value = max(max(ridgeR.predict(x_validation)), max(y_validation), max(ridgeR.predict(x_train)), max(y_train))\n",
    "plt.plot([min_value, max_value],[min_value, max_value], '--', label=\"parity\")\n",
    "plt.xlabel('Predicted Values')\n",
    "plt.ylabel('True Values')\n",
    "ax.set_box_aspect(1)\n",
    "sns.despine()\n",
    "plt.title('Parity Plot')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine coefficient and corresponding variable names\n",
    "ridge_coefficient = pd.DataFrame()\n",
    "ridge_coefficient[\"Columns\"]= x_train.columns\n",
    "ridge_coefficient['Coefficient Estimate'] = pd.Series(ridgeR.coef_)\n",
    "#ridge_coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the magnitude of coefficients\n",
    "fig, ax = plt.subplots(figsize =(7, 7))\n",
    "color =['midnightblue', 'navy', 'darkblue', 'mediumblue', 'blue', \n",
    "        'darkslateblue','mediumslateblue', 'slateblue', \n",
    "        'mediumpurple', 'blueviolet', 'darkviolet']\n",
    "\n",
    "ax.bar(ridge_coefficient[\"Columns\"], ridge_coefficient['Coefficient Estimate'], color = color)\n",
    " \n",
    "ax.spines['bottom'].set_position('zero')\n",
    " \n",
    "plt.style.use('ggplot')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3) LASSO Regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train a LASSO regression model\n",
    "from sklearn.linear_model import Lasso\n",
    " \n",
    "lasso = Lasso(alpha = 0.2)\n",
    "lasso.fit(x_train, y_train)\n",
    "\n",
    "# Predict validation data\n",
    "lasso_y_pred = lasso.predict(x_validation)\n",
    " \n",
    "print(\"LASSO regression R2 (training data):\",f\"{lasso.score(x_train, y_train):.3f}\")\n",
    "print(\"LASSO regression Q2 (validation data):\",f\"{lasso.score(x_validation, y_validation):.3f}\" \"\\n\")\n",
    "\n",
    "sum_squared_error = round(np.sum((lasso.predict(x_train) - y_train)**2),3)\n",
    "print(\"SSE on training data: \", sum_squared_error)\n",
    " \n",
    "sum_squared_error = round(np.sum((lasso.predict(x_validation) - y_validation)**2),3)\n",
    "print(\"SSE on validation data: \", sum_squared_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parity plot\n",
    "\n",
    "plt.plot(lasso.predict(x_validation), y_validation, '.', label='validation')\n",
    "plt.plot(lasso.predict(x_train), y_train, '.', label='train')\n",
    "\n",
    "min_value = min(min(lasso.predict(x_validation)), min(y_validation), min(lasso.predict(x_train)), min(y_train))\n",
    "max_value = max(max(lasso.predict(x_validation)), max(y_validation), max(lasso.predict(x_train)), max(y_train))\n",
    "plt.plot([min_value, max_value],[min_value, max_value], '--', label=\"parity\")\n",
    "plt.xlabel('Predicted Values')\n",
    "plt.ylabel('True Values')\n",
    "ax.set_box_aspect(1)\n",
    "sns.despine()\n",
    "plt.title('Parity Plot')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine coefficient and corresponding variable names\n",
    "lasso_coeff = pd.DataFrame()\n",
    "lasso_coeff[\"Columns\"] = x_train.columns\n",
    "lasso_coeff['Coefficient Estimate'] = pd.Series(lasso.coef_)\n",
    " \n",
    "#lasso_coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the magnitude of coefficients\n",
    "fig, ax = plt.subplots(figsize =(7, 7))\n",
    "color =['midnightblue', 'navy', 'darkblue', 'mediumblue', 'blue', \n",
    "        'darkslateblue','mediumslateblue', 'slateblue', \n",
    "        'mediumpurple', 'blueviolet', 'darkviolet']\n",
    "\n",
    "ax.bar(lasso_coeff[\"Columns\"], lasso_coeff['Coefficient Estimate'], color = color)\n",
    " \n",
    "ax.spines['bottom'].set_position('zero')\n",
    " \n",
    "plt.style.use('ggplot')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (4) Elastic Net Regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train and Elastic Net model\n",
    "from sklearn.linear_model import ElasticNet\n",
    " \n",
    "e_net = ElasticNet(alpha = 0.2)\n",
    "e_net.fit(x_train, y_train)\n",
    " \n",
    "# Predict validation data\n",
    "elastic_y_pred = e_net.predict(x_validation)\n",
    "\n",
    "print(\"Elastic net regression R2 (training data):\",f\"{e_net.score(x_train, y_train):.3f}\")\n",
    "print(\"Elastic net regression Q2 (validation data):\",f\"{e_net.score(x_validation, y_validation):.3f}\" \"\\n\")\n",
    "\n",
    "sum_squared_error = round(np.sum((e_net.predict(x_train) - y_train)**2),3)\n",
    "print(\"SSE on training data: \", sum_squared_error)\n",
    " \n",
    "sum_squared_error = round(np.sum((e_net.predict(x_validation) - y_validation)**2),3)\n",
    "print(\"SSE on validation data: \", sum_squared_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parity plot\n",
    "\n",
    "plt.plot(e_net.predict(x_validation), y_validation, '.', label='validation')\n",
    "plt.plot(e_net.predict(x_train), y_train, '.', label='train')\n",
    "\n",
    "min_value = min(min(e_net.predict(x_validation)), min(y_validation), min(e_net.predict(x_train)), min(y_train))\n",
    "max_value = max(max(e_net.predict(x_validation)), max(y_validation), max(e_net.predict(x_train)), max(y_train))\n",
    "plt.plot([min_value, max_value],[min_value, max_value], '--', label=\"parity\")\n",
    "plt.xlabel('Predicted Values')\n",
    "plt.ylabel('True Values')\n",
    "ax.set_box_aspect(1)\n",
    "sns.despine()\n",
    "plt.title('Parity Plot')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine coefficient and corresponding variable names\n",
    "e_net_coeff = pd.DataFrame()\n",
    "e_net_coeff[\"Columns\"] = x_train.columns\n",
    "e_net_coeff['Coefficient Estimate'] = pd.Series(e_net.coef_)\n",
    "#e_net_coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the magnitude of coefficients\n",
    "fig, ax = plt.subplots(figsize =(7, 7))\n",
    "color =['midnightblue', 'navy', 'darkblue', 'mediumblue', 'blue', \n",
    "        'darkslateblue','mediumslateblue', 'slateblue', \n",
    "        'mediumpurple', 'blueviolet', 'darkviolet']\n",
    "\n",
    "ax.bar(e_net_coeff[\"Columns\"], e_net_coeff['Coefficient Estimate'], color = color)\n",
    " \n",
    "ax.spines['bottom'].set_position('zero')\n",
    " \n",
    "plt.style.use('ggplot')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (5) Compare model architectures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "k = 5\n",
    "reg_alpha = 0.2\n",
    "\n",
    "lm_lreg = LinearRegression()\n",
    "scores_lreg = cross_val_score(lm_lreg, X, Y, cv=k)\n",
    "np.mean(scores_lreg)\n",
    "print(\"MLR k-fold cross-validation Q2:\",f\"{np.mean(scores_lreg):.3f}\" \"\\n\")\n",
    "\n",
    "lm_ridge = linear_model.Ridge(alpha=reg_alpha)\n",
    "scores_ridge = cross_val_score(lm_ridge, X, Y, cv=k)\n",
    "np.mean(scores_ridge)\n",
    "print(\"Ridge k-fold cross-validation Q2:\",f\"{np.mean(scores_ridge):.3f}\" \"\\n\")\n",
    "\n",
    "lm_lasso = linear_model.Lasso(alpha=reg_alpha)\n",
    "scores_lasso = cross_val_score(lm_lasso, X, Y, cv=k)\n",
    "np.mean(scores_lasso)\n",
    "print(\"LASSO k-fold cross-validation Q2:\",f\"{np.mean(scores_lasso):.3f}\" \"\\n\")\n",
    "\n",
    "lm_en = linear_model.ElasticNet(alpha=reg_alpha)\n",
    "scores_en = cross_val_score(lm_en, X, Y, cv=k)\n",
    "np.mean(scores_en)\n",
    "print(\"Elastic Net k-fold cross-validation Q2:\",f\"{np.mean(scores_en):.3f}\" \"\\n\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (6) Observe impact of regularization strength "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for alpha in [0,0.01,0.02,0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 6, 8]:\n",
    "    lm_ridge = linear_model.Ridge(alpha=alpha)\n",
    "    scores_ridge = cross_val_score(lm_ridge, X, Y, cv=k)\n",
    "    mean_ridge = round(np.mean(scores_ridge),2)\n",
    "    print(f\"Alpha = {alpha} results in average Ridge Q2 = {mean_ridge}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for alpha in [0,0.01,0.02,0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 6, 8]:\n",
    "    lm_lasso = linear_model.Lasso(alpha=alpha)\n",
    "    scores_lasso = cross_val_score(lm_lasso, X, Y, cv=k)\n",
    "    mean_lasso = round(np.mean(scores_lasso),2)\n",
    "    print(f\"Alpha = {alpha} results in average LASSO Q2 = {mean_lasso}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for alpha in [0,0.01,0.02,0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 6, 8]:\n",
    "    lm_en = linear_model.ElasticNet(alpha=alpha)\n",
    "    scores_en = cross_val_score(lm_en, X, Y, cv=k)\n",
    "    mean_en = round(np.mean(scores_en),2)\n",
    "    print(f\"Alpha = {alpha} results in average EN Q2 = {mean_en}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}