linear_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply Jupyter notebook style\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "from custom.styles import style_string\n",
    "\n",
    "HTML(style_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align:center;\">\n",
    "  <img src=\"custom/molssi_main_horizontal.png\" style=\"display: block; margin: 0 auto; max-height:200px;\">\n",
    "</div>\n",
    "\n",
    "# Linear Fitting\n",
    "\n",
    "<div class=\"overview admonition\"> \n",
    "<p class=\"admonition-title\">Overview</p>\n",
    "\n",
    "Questions:\n",
    "\n",
    "* How can I use Python to fit a linear model?\n",
    "\n",
    "Objectives:\n",
    "\n",
    "* Fit a linear model using statsmodels\n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "A common exercise in introductory chemistry is fitting a linear model to absorbance measurements at a particular wavelength vs. concentration of an analyte using spectrophotometry. \n",
    "In these cases, the absorbance of a sample is described using [Beer-Lambert Law](https://en.wikipedia.org/wiki/Beer%E2%80%93Lambert_law)\n",
    "\n",
    "$$\n",
    "A = \\varepsilon \\cdot l \\cdot c\n",
    "$$\n",
    "\n",
    "where $A$ is the absorbance, $\\varepsilon$ is the absorptivity constant, $l$ is path length, and $c$ is the concentration of the sample.\n",
    "\n",
    "\n",
    "In this notebook, we will read in absorbance data, fit a linear model using a Python library called [statsmodels](https://www.statsmodels.org/stable/index.html), and calculate the concentrations of unknowns using our model.\n",
    "We will utilize skills learned in our previous lessons, including reading, accessing, and analyzing data using `pandas` and visualizing data using `plotly`.\n",
    "We will add a new skill - fitting a linear model with a Python library called `statsmodels`. \n",
    "\n",
    "## Importing Libraries and Reading Data\n",
    "\n",
    "Note that when fitting a model in Python, there are a number of options for the library you might pick including NumPy, SciPy, SciKit-Learn, and statsmodels.\n",
    "The library that you pick might be based on personal preference or features that a particular library offers.\n",
    "In this notebook, we show fitting with `statsmodels` because of the ease of seeing fit statics and defining a formula for fitting.\n",
    "\n",
    "For `statsmodels`, we will use the formula API. \n",
    "This will allow us to define the relationship we'd like to fit using a formula as a string.\n",
    "\n",
    "To start this notebook, we will import the libraries we need. \n",
    "We will import the following:\n",
    "\n",
    "Library Name | Purpose\n",
    "-------------|--------\n",
    "[pandas](https://pandas.pydata.org/docs/)       | reading and processing data\n",
    "[plotly](https://plotly.com/python/getting-started/)       | interactive plots and visualization\n",
    "[statsmodels](https://www.statsmodels.org/stable/index.html)  | data fitting\n",
    "\n",
    "\n",
    "The cell below imports the libraries we will use for our analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports \n",
    "import pandas as pd # for reading data from a file and putting in a table\n",
    "import plotly.express as px # for plotting\n",
    "\n",
    "import statsmodels.formula.api as smf # for fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the example, we will use data stored in the file `data/protein_assay.csv`.\n",
    "This data represents absorbance data recorded in a Bradford Assay for determining protein concentration.\n",
    "\n",
    "After we have imported our libraries, we will next use `pandas` to read in our data.\n",
    "Our data is stored in a comma separated value (CSV) file, though pandas can also read from Excel files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in file\n",
    "df = pd.read_csv(\"data/protein_assay.csv\")\n",
    "\n",
    "# View the first five rows.\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing Data\n",
    "\n",
    "After reading in our data, we might wish to inspect it visually. \n",
    "We can do this using plotly express (imported as `px`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure using px.scatter.\n",
    "fig = px.scatter(df, x='concentration', y=\"absorbance\")\n",
    "\n",
    "# Show the figure.\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you can change the axis labels using the following syntax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure using px.scatter.\n",
    "fig = px.scatter(df, x='concentration', y=\"absorbance\", \n",
    "                 labels={\"concentration\":\"Concentration (mg/mL)\", \"absorbance\":\"Absorbance at 595 nm\"})\n",
    "\n",
    "# Show the figure.\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a Linear Model\n",
    "\n",
    "To fit a linear equation to our data, we can use a library called [statsmodels](https://www.statsmodels.org/stable/index.html). \n",
    "`statsmodels` is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration.\n",
    "\n",
    "While there are numerous options for fitting data in Python, `statsmodels` is particularly beneficial when dealing with linear models. \n",
    "Note that this encompasses more than just \"linear equations\". \n",
    "This library is especially handy when you need straightforward access to fit parameters and in-depth statistical metrics related to the fit.\n",
    "\n",
    "In particular, we are using a part of `statsmodels` called the formula API.\n",
    "The formula API lets you define a formula as a string for fitting and is specifically designed to work with dataframes.\n",
    "When defining a formula, you use the column names in a string to define the relationship.\n",
    "Note that the column names should not have spaces, or entering the relationship is a bit more complicated.\n",
    "For example, if we had data representing pressure and temperature at constant volume and expected them to follow a linear relationship such as for an ideal gas\n",
    "\n",
    "$$\n",
    "P = \\frac{n R T}{V}\r",
    "$$\n",
    "\n",
    "we would write `\"P ~ T\"` when using the `statsmodels` formula API. \n",
    "This is assuming that our dataframe has columns named `P` and `T`\n",
    "\n",
    "\n",
    "As a slightly more complicatd eexample lone could also fit something like a polynomial using `\"y ~ np.power(x, 2) + x\"` if you had imported NumPy (`import numpy as n`). \n",
    "\n",
    "To use the formula API, we will use `smf` (imported in first cell). \n",
    "We will use ordinary least squares (`ols`) for our fit, though [a number of other options are offered](https://www.statsmodels.org/dev/api.html#statsmodels-formula-api).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression = smf.ols(\"absorbance ~ concentration\", data=df).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see a summary of the fit including the `R-squared` by using the `.summary()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we would like to force the model to not have an intercept, we use a special `statsmodel` syntax. \n",
    "Adding `-1` to our formula forces the intercept to be 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# force the intercept to be 0\n",
    "regression = smf.ols(\"absorbance ~ concentration - 1\", data=df).fit()\n",
    "regression.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model parameters are in `.params` of the fit variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the slope, we will get the coefficient in front of the concentration variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = regression.params[\"concentration\"]\n",
    "print(slope)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see what our model predicts for our input concentration values by using the `predict` method.\n",
    "In the cell below, we save the results in a new column in our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"predicted\"] = regression.predict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter(x= df[\"concentration\"], y=df[\"absorbance\"])\n",
    "fig.add_scatter(x=df[\"concentration\"], y=df[\"predicted\"], mode=\"lines\", name=\"model\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting unknowns\n",
    "\n",
    "Now that we have our model and slope, we can use it to calculate the protein concentrations for a set of unknowns.\n",
    "\n",
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Use pandas to read in the file `data/protein_samples.csv`.\n",
    "\n",
    "Next, use our calculated slope to predict concentration based on measured absorbance.\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise - Perform a Linear Fit\n",
    "\n",
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Now that you have completed a simple linear regression exercise with protein assay data, here is a problem with a slightly larger dataset, taken from a ground water survey of wells in Texas kindly provided by Houghton-Mifflin. \n",
    "The data for this exercise is in the file `data/ground_water.csv`.\n",
    "Using the skills you have learned with pandas and `statsmodels`, get the linear regression statistics for the relationship between pH (dependent variable) and bicarbonate levels (ppm in well water in Texas; independent variable).\n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}