03_nonlinear_regression_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "vocal-lemon",
   "metadata": {},
   "source": [
    "Nonlinear Regression Part 2\n",
    "========================="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "japanese-vaccine",
   "metadata": {},
   "source": [
    "<div class=\"overview-this-is-a-title overview\">\n",
    "<p class=\"overview-title\">Overview</p>\n",
    "<p>Questions</p>\n",
    "    <ul>\n",
    "        <li> How can I perform statistically robust non-linear regression in a Jupyter notebook?</li>\n",
    "        <li> How can I prepare accurate and attractive plots of my regression results that include the best-fit line?\n",
    "    </ul>\n",
    "<p>Objectives:</p>\n",
    "    <ul>\n",
    "        <li> Import data from a csv file\n",
    "        <li> Create a function based on the equation to be fitted\n",
    "        <li> Generate the best-fit curve with parameters and statistical descriptors\n",
    "        <li> Create plots of the data\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loaded-tomato",
   "metadata": {},
   "source": [
    "## Data Fitting\n",
    "In an earlier module, we used the `linregress` function from the scipy.stats library to perform linear regression. This is sufficient for a least squares linear fit to the data. Python contains tools that also enable us to fit data to an equation that we have defined, with two or more paramaters. So it is ideal for nonlinear curve fitting to the Michaelis-Menten equation.\n",
    "\n",
    "In this module, we will use some familiar tools: importing csv data into a pandas dataframe and plotting results using the matplotlib.pyplot library. We will also create a function based on the Michaelis-Menten equation. The exciting part is learning how to fit the data to the function using `curve_fit` from scipy.optimize.\n",
    "\n",
    "Here are the steps in our process:\n",
    "\n",
    "1. Import the libraries and functions that you need, including `curve_fit` from scipy.optimize.\n",
    "1. Set up a pandas dataframe from the MM_data.csv file you saved from the previous module.\n",
    "1. Create your independent and dependent variables for curve fitting.\n",
    "1. Create a function for the equation you want to fit. The function should accept as inputs the independent variable(s) and all the parameters to be fit.\n",
    "1. Use the function `curve_fit` from the scipy.optimize library to fit your data.\n",
    "1. Extract the fit parameters from the output of `curve_fit`.\n",
    "1. Use your function to calculate y values using your fit model to see how well your model fits the data.\n",
    "1. Graph your original data and the best-fit equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "opposite-marathon",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the libraries and functions that you need, including curve_fit from scipy.optimize.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "moderate-russell",
   "metadata": {},
   "source": [
    "## The Michaelis-Menten equation\n",
    "\n",
    "The Michaelis-Menten equation (which is probably very familiar to you) describes kinetics for enzyme-catalyzed reactions. We will use the simplest version of the equation, based on initial reaction velocities with a single substrate.\n",
    "\n",
    "In our equation,  there are two parameters, V$_{max}$ (maximum velocity) and K$_m$ (the Michaelis constant, which equals the substrate concentration at half maximal velocity), an independent variable, \\[S\\] = substrate concentration, and a dependent variable, V$_i$, the initial velocity.\n",
    "\n",
    "$$\n",
    "V_i = \\frac{V_{max} * [S]}{K_m + [S]}\n",
    "$$\n",
    "\n",
    "We will define the two parameters and two variables for our equation:\n",
    "\n",
    "- Vmax = the maximum velocity (the first parameter)\n",
    "- Km = the Michaelis constant (the second parameter)\n",
    "- SubConc = substrate concentration (the independent variable)\n",
    "- Vi = initial velocity (the dependent variable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "traditional-sitting",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up a pandas dataframe from the MM_data.csv file you saved from the previous module.\n",
    "# If you skipped nonlinear_regression_part_1, you can simply use the MM_data.csv\n",
    "# file that is provided.\n",
    "# create the filehandle\n",
    "\n",
    "\n",
    "# import the data\n",
    "\n",
    "\n",
    "# inspect the data to make sure you have substrate concentrations and initial velocities\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "medium-heating",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create your independent and dependent variables for curve fitting.\n",
    "\n",
    "\n",
    "# Inspect your data to make sure it looks like a Michaelis-Menten plot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "disturbed-lottery",
   "metadata": {},
   "source": [
    "## Creating Functions\n",
    "\n",
    "To this point we have called a number of functions that are included with python libraries such as pandas, numpy and matplotlib. Most code is organized into blocks of code which perform a particular task. These code blocks are called functions. A commercial software package likely has hundreds, thousands or millions of functions. Functions break up our code into smaller, more easily understandable statements, and also allow our code to be more modular, meaning we can take pieces and reuse them. In general, each function should perform only one computational task. \n",
    "\n",
    "In Python, the following syntax is used to declare a function:\n",
    "\n",
    "```python\n",
    "def function_name(independent_variable, parameter1, parameter2):\n",
    "    ** function body code **\n",
    "    return value_to_return\n",
    "```\n",
    "\n",
    "Functions are defined using the `def` keyword, followed by the name of the function. The function may have independent variable(s) and parameters which are passed to it; these are in parenthesis following the function name. You then add a colon after the parentheses. This causes the function body code to indent in a Jupyter notebook coding cell. In certain cases, a function may have no parameters as well. Most (though not all) functions return some type of information. It is important to note here that defining a function does not execute it.\n",
    "\n",
    "Now we will create a function for the Michaelis-Menten equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "surprising-length",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a function for the equation you want to fit. The function should accept \n",
    "# as inputs the independent variable(s) and all the parameters to be fit.\n",
    "\n",
    "\n",
    "# Test the function to see if it produces the expected output\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "insured-advocate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Common errors and how to read them"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "premium-sydney",
   "metadata": {},
   "source": [
    "In the previous cell, we received an error message that resulted when we used parentheses by mistake rather than square brackets in our attempt to assign a series from a dataframe to a variable. The error message is\n",
    "\n",
    "```python\n",
    "TypeError: 'DataFrame' object is not callable\n",
    "```\n",
    "\n",
    "The term `callable` here refers to something that happens with a function.\n",
    "\n",
    "Here are two very useful sources when you encounter errors: [Stack Overflow](https://stackoverflow.com) and [GeeksforGeeks](https://www.geeksforgeeks.org). In fact, if you just paste your error in quotes (e.g., \"TypeError: 'DataFrame' object is not callable\") into a search engine, you are likely to find hits from both sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7988cc69-96f7-482f-b39b-5ac32ba3c27c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     34.272718\n",
       "1     33.542678\n",
       "2     31.362310\n",
       "3     27.801416\n",
       "4     22.829652\n",
       "5     16.417302\n",
       "6     13.522863\n",
       "7      9.096204\n",
       "8      5.096143\n",
       "9      5.204653\n",
       "10     1.487957\n",
       "Name: initial velocities, dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Correcting the error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "powerful-cigarette",
   "metadata": {},
   "source": [
    "## Using curve_fit\n",
    "The syntax for the curve_fit function from scipy.optimize is\n",
    "```python\n",
    "popt, pcov = curve_fit(function, xdata, ydata)\n",
    "```\n",
    "| Term | Signficance |\n",
    "| :----: | :---- |\n",
    "| popt | optimized parameters |\n",
    "| pcov | a 2D matrix of the covariance for the curve fit |\n",
    "| curve_fit | the function from scipy.optimize |\n",
    "| function | the function you choose that curve_fit will use for fitting the data |\n",
    "| xdata | a data array for the independent variable |\n",
    "| ydata | a data array for the dependent variable |\n",
    "\n",
    "You can learn much more about `curve_fit` in the online documentation for [scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html). \n",
    "\n",
    "Let's review the function, `MM`, that we defined above, to clearly distinguish the data, variables and parameters associated with `curve_fit` from those associated with `MM`.\n",
    "\n",
    "```python\n",
    "def MM(SubConc, Vmax, Km):\n",
    "    Vi = (Vmax * SubConc)/(Km + SubConc)\n",
    "    return Vi\n",
    "```\n",
    "\n",
    "You need to pass three pieces of data to `MM` to predict the initial velocity for an enzyme-catalyzed reaction: the substrate concentration, and values for V$_{max}$ and K$_m$.\n",
    "\n",
    "scipy.optimize.curve_fit() can then fit the data to the `MM` function you have defined. In an earlier module, there was a discussion of different ways to linearize data so that linear regression could be used, for example, taking logs of the molecular weights of protein standards to plot versus relative migration for an SDS-PAGE gel. With curve_fit, you can do a least squares regression fit for data that demonstrate nonlinear relationships."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "asian-indonesia",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the function curve_fit to fit your data.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "automotive-fiction",
   "metadata": {},
   "source": [
    "## Extracting the output from curve_fit\n",
    "There are two outputs from curve_fit - popt and pcov. These are typical names that help to point to the contents of the variable, but you can choose any variable name you like.\n",
    "1. popt = a list of the optimized parameters. In our case, the two parameters are $V_{max}$ and $K_m$. Remember that python starts counting at 0, so popt[0] = $V_{max}$ and popt[1] = $K_m$. \n",
    "1. pcov = an n-dimensional array of the covariance for the parameters, where n is the number of parameters. In our case, n = 2. \n",
    "\n",
    "So we can extract the fitted values for $V_{max}$ and $K_m$ just by setting variables equal to them. The covariance matrix can be used to extract standard errors or standard deviations for each of the parameters. In our case, we will calculate one standard deviation for each parameter. Note that we use two numpy functions (sqrt and diag) to calculate the standard deviation values from the 2D matrix for pcov."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "oriental-primary",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the fit parameters from the output of curve_fit.\n",
    "\n",
    "\n",
    "\n",
    "# Extract one standard deviation for each parameter\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expanded-illinois",
   "metadata": {},
   "source": [
    "## Review of f string printing\n",
    "\n",
    "It is unlikely that you will want to have your values printed out to more than 10 decimal places, which is the default when python prints floats. There are a number of ways to get around this that include truncation or rounding, but we are going to look at a very useful tool - [f string printing](https://realpython.com/python-f-strings/). Here is the syntax for an f string print statement:\n",
    "\n",
    "```python\n",
    "print(f'a string you want to print that may contain {variables in curly brackets}')\n",
    "```\n",
    "Variables of type float can also be assigned precision using this syntax:\n",
    "\n",
    "```python\n",
    "print(f'Vmax for this reaction: {Vmax:.2f}')\n",
    "```\n",
    "\n",
    "Since Vmax is in curly brackets, the print statement will insert the value for the variable, Vmax, at that point in the print statement. Vmax is followed by `:.2f` which directs the print statement to print the value of Vmax to 2 decimal places. We will now use f string printing to control the precision on the values for our parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "military-camera",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fatty-bouquet",
   "metadata": {},
   "source": [
    "<div class=\"exercise-this-is-a-title exercise\">\n",
    "    <p class=\"exercise-title\"><b>Exercise</b></p>\n",
    "    <p> Print the value for K<sub>m</sub> +/- one standard deviation to three decimal places. </p>\n",
    "\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "further-point",
   "metadata": {},
   "source": [
    "## Plotting the data\n",
    "At this point, you want to compare your plotted data with the curve based on the parameters from curve_fit. We'll do this in two stages. First, we'll use the `MM` function to calculate y values using your fit model to see how well your model fits the data. \n",
    "\n",
    "We will use matplotlib.pyplot (please note the import statement at the top of this notebook) to create and annotate this plot. Our first plot will contain two sets of data: \n",
    "\n",
    "1. the actual data for initial velocity vs. substrate concentration from our pandas dataframe\n",
    "2. the fitted data using the parameters generated by curve_fit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "together-cycle",
   "metadata": {},
   "outputs": [],
   "source": [
    "# using your fit model to see how well your model fits the data.\n",
    "\n",
    "# Graph your original data and the fit equation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recreational-melissa",
   "metadata": {},
   "source": [
    "## Smoothing the plot\n",
    "You probably noticed that the curve appears a bit disjointed, because it provides values only at the selected concentrations of the substrate. A fit curve based on 1000 x-values will be much smoother. To do this, we first use the numpy linspace function to create a new list of 1000 x-values ranging from the minimum xdata value (using the numpy min function) to the maximum xdata value (using the numpy max function). Then we generate a fit curve using the MM function. We then add the experimental data and the smoothed fit data to the plot. Note that the experimental data is plotted using `'bo'`. This tells matplotlib.pyplot to use blue (`b`) circles (`o`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "retired-costume",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To generate a smoother fitted curve\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "agreed-cancellation",
   "metadata": {},
   "source": [
    "<div class=\"exercise-this-is-a-title exercise\">\n",
    "    <p class=\"exercise-title\"><b>Exercise</b></p>\n",
    "    <p> Use what you learned in previous modules to label the x-axis and y-axis and to insert the values for $V_{max}$ and $K_m$ with their standard deviations in the white space on the graph. Change the symbol for the substrate concentrations to a red triangle. </p>\n",
    "    \n",
    "  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "valid-example",
   "metadata": {},
   "source": [
    "<div class=\"exercise-this-is-a-title exercise\">\n",
    "<p class=\"exercise-title\">Take Home Exercise</p>\n",
    "    <p> If you want to dig deeper, you can complete the curve fitting analysis of the chymotrypsin data from the previous module. You may need to adjust some of the details to generate high quality output.</p>\n",
    "    \n",
    "\n",
    "\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "controversial-pointer",
   "metadata": {},
   "source": [
    "## Common errors \n",
    "\n",
    "Python is very particular about the use of parentheses (), square brackets [], and braces or curly brackets {}. \n",
    "\n",
    "Parentheses are used to pass arguments to functions. For example, \n",
    "```python\n",
    "Vi = MM(xdata, Vmax, Km)\n",
    "```\n",
    "\n",
    "Square brackets are used to create lists. In this example, the square brackets are used to call a series from the rates_df dataframe.\n",
    "```python\n",
    "xdata = rates_df['pNPP (mM)']\n",
    "```\n",
    "\n",
    "Curly brackets or braces are used to create dictionaries or sets. \n",
    "```python\n",
    "a_dict = {'Brady': 12, 'Evans': 13}\n",
    "```\n",
    "You can learn more about [dictionaries in the python online manual](https://docs.python.org/3/library/stdtypes.html#mapping-types-dict). \n",
    "\n",
    "Here are two very useful sources when you encounter errors: [Stack Overflow](https://stackoverflow.com) and [GeeksforGeeks](https://www.geeksforgeeks.org). In fact, if you just paste your error in quotes (e.g., \"TypeError: 'DataFrame' object is not callable\") into a search engine, you are likely to find hits from both sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "625a5891",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}