Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "38cf12af-4bf7-42c7-851e-f89a0062139b",
   "metadata": {},
   "source": [
    "---\n",
    "title: Regression\n",
    "project:\n",
    "  type: website\n",
    "format:\n",
    "  html:\n",
    "    code-fold: true\n",
    "    code-tools: true\n",
    "jupyter: python 3\n",
    "number-sections: false\n",
    "filters:\n",
    "    - pyodide\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff6bd615-5678-4a1c-bb21-c5e05b17d776",
   "metadata": {},
   "source": [
    "In this notebook, we'll implement an example application of linear regression applied to behavioral and electrophysiologicial data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04aa2ddd-5bd5-42d3-9543-f747df122830",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Load modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy.stats import pearsonr\n",
    "import statsmodels.api as sm\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf3b269b-e554-4451-88bd-1c7670fc3ec5",
   "metadata": {},
   "source": [
    "# Load the data\n",
    "\n",
    "The example behavioral and electrophysiologicial data is [regression_example_data.csv](https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/Data/regression_example_data.csv). Get these data, and load them:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ecb4bbf-a932-4ea5-b2c1-af514beb8dbb",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "df = pd.read_csv(\"https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/Data/regression_example_data.csv\")\n",
    "\n",
    "# Extract the variables from the loaded data\n",
    "task_performance = np.array(df.iloc[:,0])  #Get the values associated with the first column of the dataframe\n",
    "firing_rate = np.array(df.iloc[:,1])  #Get the values associated with the second column of the dataframe\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d74a68de-1492-43fc-a84d-38bf546d74d7",
   "metadata": {},
   "source": [
    "# Visualize the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f20129e5-89c6-403b-8862-7c7977320285",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Plot it ...\n",
    "plt.figure()\n",
    "plt.plot(firing_rate, task_performance, '.')\n",
    "plt.xlabel('Firing rate [Hz]')\n",
    "plt.ylabel('Task Performance [a.u.]')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "854d3e14-5334-476f-9086-36a404489489",
   "metadata": {},
   "source": [
    "# Correlation\n",
    "\n",
    "Compute the [correlation](https://numpy.org/doc/stable/reference/generated/numpy.correlate.html) between $x$ and $y$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "919171a8-1186-4585-b07b-0121222bafe4",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "N       = np.size(firing_rate)\n",
    "x       = firing_rate - np.mean(firing_rate)\n",
    "y       = task_performance - np.mean(task_performance)\n",
    "sigma_x = 'SOMETHING'    #Standard deviation of x\n",
    "sigma_y = 'SOMETHING'    #Standard deviation of y\n",
    "\n",
    "correlation = 'SOMETHING'\n",
    "print(correlation)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "925ca081-0a27-4092-8886-cf28795297d6",
   "metadata": {},
   "source": [
    "# Regression (compute it)\n",
    "\n",
    "Model the data using regression."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd0cd49c-50a5-4b4d-be34-031116030ac0",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "data = {\"x\": firing_rate, \"y\": task_performance}\n",
    "\n",
    "res1 = ols(\"y ~1 + x\", data=data).fit()\n",
    "res1.summary()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a376d89e-9939-4f39-9455-994eb58a2874",
   "metadata": {},
   "source": [
    "# Regression (plot it)\n",
    "\n",
    "Plot the estimated regression model with the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "426aeb05-7f58-41f0-a885-35e25ff9c2f3",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Get model prediction.\n",
    "fitted_values = res1.fittedvalues\n",
    "\n",
    "# Sort x values for better plotting of the regression line\n",
    "x_sorted = np.sort(firing_rate)\n",
    "fitted_sorted = np.sort(fitted_values)\n",
    "\n",
    "# Plot the regression line (fitted model)\n",
    "plt.figure()\n",
    "plt.scatter(firing_rate,task_performance)\n",
    "plt.plot(x_sorted, fitted_sorted, label=\"Fitted Model\", color=\"red\")\n",
    "\n",
    "plt.xlabel('Firing rate [Hz]')\n",
    "plt.ylabel('Task Performance [a.u.]')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "961dacc7-0a12-4568-91f6-63f72fac5110",
   "metadata": {},
   "source": [
    "---\n",
    "# Regression example (Part 2)\n",
    "\n",
    "We learn that an additional predictor - age - impacts task performance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b580395d-64fc-467f-a615-3103aeef424c",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "df = pd.read_csv(\"https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/Data/regression_example_data.csv\")\n",
    "\n",
    "# Extract the variables from the loaded data\n",
    "task_performance = np.array(df.iloc[:,0])\n",
    "firing_rate = np.array(df.iloc[:,1])\n",
    "age = np.array(df.iloc[:,2])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aa87b33-cc9a-4b05-b24d-28e017f086ee",
   "metadata": {},
   "source": [
    "# Visualize the new data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f55b361-8a13-4b04-8638-ece113f2cf6c",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Plot it ...\n",
    "plt.figure()\n",
    "plt.plot(age, task_performance, '.')\n",
    "plt.xlabel('Age [months]')\n",
    "plt.ylabel('Task Performance [a.u.]')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "778f0827-6d6f-4812-9b17-0273ef64f40d",
   "metadata": {},
   "source": [
    "# Correlation (between task performance and age)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e884254-06b3-4ef0-8ae4-b9c0d1e76cf4",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Compute the correlation between task performance and age\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29006d51-cf15-408a-a167-36f1544da0f5",
   "metadata": {},
   "source": [
    "# Visualize all data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f878ce-3614-48fc-812a-fb3d70685564",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "fig = plt.figure(figsize=(12, 12))\n",
    "ax  = fig.add_subplot(projection='3d')\n",
    "ax.scatter(age, firing_rate, task_performance)\n",
    "ax.set_xlabel('Age [months]')\n",
    "ax.set_ylabel('Firing Rate [Hz]')\n",
    "ax.set_zlabel('Task Performance');\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04d9c052-2811-48b7-a88d-1f2e7dc6b65c",
   "metadata": {},
   "source": [
    "# Regression (compute it with all data)\n",
    "\n",
    "Model all data using regression."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c63f625-6f02-4f3c-a61d-98d899eb72b6",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "from statsmodels.formula.api import ols\n",
    "data = {\"firing_rate\": firing_rate, \"age\": age, \"y\": task_performance}\n",
    "# Write the model and print out the summary\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb1ba49e-db3b-4b7a-b216-e568556d174b",
   "metadata": {},
   "source": [
    "# Regression (plot it with all data)\n",
    "\n",
    "Plot the estimated regression model with two predictors and all data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abd064b7-07d7-4133-bc51-dca781209433",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "# Plot the mean model fit.\n",
    "\n",
    "# First, plot the data.\n",
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax  = fig.add_subplot(projection='3d')\n",
    "ax.set_xlabel('Age [months]')\n",
    "ax.set_ylabel('Firing Rate [Hz]')\n",
    "ax.set_zlabel('Task Performance')\n",
    "ax.scatter(age, firing_rate, task_performance)\n",
    "\n",
    "# Then, define model parameter estimates.\n",
    "# REPLACE THESE VALUES WITH YOUR PARAMETER ESTIMATES\n",
    "alpha  = 1\n",
    "beta_1 = 1\n",
    "beta_2 = 1\n",
    "\n",
    "# Finally, plot the model fit.\n",
    "x      = np.arange(8, 12, 0.1)          # Firing rate\n",
    "y      = np.arange(10,20, 0.1)          # Age\n",
    "xx, yy = np.meshgrid(x, y)              # Two dim coordinates\n",
    "zz     = alpha + beta_1*xx + beta_2*yy  # Model predictions\n",
    "ax.plot_surface(yy,xx,zz);\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5644779f-0f7d-4293-87e4-ebdd60d10478",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}