week06a-multivariate-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4b1c1805-e6d7-44b0-873c-844bf74c9ac8",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Multivariate regression\n",
    "\n",
    "__Goal:__ Find a relationship that explains variable $y$ in terms of variables, $x_1, x_2, x_3$,...$x_n$\n",
    "\n",
    "![images/mult_reg_viz.png](images/mult_reg_viz.png)\n",
    "\n",
    "source: [http://www.sjsu.edu/faculty/gerstman/EpiInfo/cont-mult.htm](http://www.sjsu.edu/faculty/gerstman/EpiInfo/cont-mult.htm)\n",
    "\n",
    "This three dimensional visualization shows how linear model based on two predictor variables, $x_1$ and $x_2$ can be used to model a response variable $y$. A constant and two slopes to define a 2D plane in 3D space. The sum of squared vertical distances between the plane (model) and observations of $y$ are minimized. Like fitting a line in 2D space, this procedure assumes the validity of a linear model.\n",
    "\n",
    "## Example: Model for aragonite saturation state based on three other oceanographic variables\n",
    "\n",
    "As a motivation for multiple linear regression, we consider a model for aragonite saturation state, $\\Omega_A$. The study by Juranek et al. (2009), discussed on greater detail below, uses multiple linear regression to model $\\Omega_A$ for surveys where it is not measured, using more commonly measured parameters. We know that $\\Omega_A$ depends on both physical and biological processes, so one candidate for a model might be:\n",
    "\n",
    "$y = b_0 + b_1x_1 + b_2x_2 + b_3x_3$\n",
    "\n",
    "where $y$ = $\\Omega_A$\n",
    "\n",
    "$x_1$ = temperature\n",
    "\n",
    "$x_2$ = salinity\n",
    "\n",
    "$x_3$ = dissolved oxygen\n",
    "\n",
    "$k$ = 3 predictor variables\n",
    "\n",
    "$N$ = 1011 samples\n",
    "\n",
    "Because we have three predictor variables, it is hard to visualize in four-dimensional space. However, the same principles are involved as fitting a 1D line in 2D space, or a 2D plane in 3D space.\n",
    "\n",
    "## Equations for linear model\n",
    "\n",
    "A linear model can be represented as a system of $N$ equations.\n",
    "\n",
    "$$ b_0 + b_1x_{11} + b_2x_{12} + b_3x_{13} = \\hat{y}_1 + \\epsilon_1$$\n",
    "\n",
    "$$ b_0 + b_1x_{21} + b_2x_{22} + b_3x_{23} = \\hat{y}_2 + \\epsilon_2$$\n",
    "\n",
    "$$ b_0 + b_1x_{31} + b_2x_{32} + b_3x_{33} = \\hat{y}_3  + \\epsilon_3$$ \n",
    "\n",
    "$$ \\ldots$$\n",
    "\n",
    "$$ b_0 + b_1x_{N1} + b_2x_{N2} + b_3x_{N3} = \\hat{y}_N  + \\epsilon_4$$\n",
    "\n",
    "where $\\hat{y}_i$ is a modeled value and $\\epsilon_i$ is the difference between the modeled value, $\\hat{y}_i$ and an observation $y_i$. The least squares regression minimizes the sum of $\\epsilon_i^2$, the overall deviation between the linear model and data.\n",
    "\n",
    "## Matrix form\n",
    "\n",
    "To solve a least squares problem numerically, it helps to write the system of equations for the model in matrix form.\n",
    "\n",
    "Form a vector $Y$ of $N$ observations.\n",
    "\n",
    "$$ Y = \\begin{bmatrix}\n",
    "        y_1 \\\\\n",
    "        y_2 \\\\\n",
    "        y_3 \\\\\n",
    "        \\vdots \\\\\n",
    "        y_n \\\\\n",
    "        \\end{bmatrix}$$\n",
    "\n",
    "A vector $B$ contains $k+1$ unknown coefficients.\n",
    "\n",
    "$$ B = \\begin{bmatrix}\n",
    "        b_0 \\\\\n",
    "        b_1 \\\\\n",
    "        b_2 \\\\\n",
    "        b_3 \\\\\n",
    "        \\end{bmatrix}$$\n",
    "\n",
    "The predictor variables are stored as columns in a matrix with $N$ rows and $k+1$ columns\n",
    "\n",
    "$$\n",
    "X =     \\begin{bmatrix}\n",
    "        1 & x_{11} & x_{12} & x_{13} \\\\\n",
    "        1 & x_{21} & x_{22} & x_{23} \\\\\n",
    "        1 & x_{31} & x_{32} & x_{33} \\\\\n",
    "        \\vdots & \\vdots & \\vdots     & \\vdots \\\\\n",
    "        1 & x_{N1} & x_{N2} & x_{N3} \\\\\n",
    "        \\end{bmatrix}$$\n",
    "\n",
    "Now the system of equations for the linear model can be written as\n",
    "\n",
    "$ X B = Y$\n",
    "\n",
    "\n",
    "## Numerical solution\n",
    "\n",
    "The least squares problem is solved using a singular value decomposition method. Efficient alorithms for this procedure are typically included in scientfic computing software. In Python, create an array for the vector $Y$ and a 2D array for the matrix $X$. Then use `np.linalg,lstsq` to solve for $B$.\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "B = np.linalg.lstsq(X,Y)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11ab16cd-cb99-4e67-a39e-524ac104fb3e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Testing for significance\n",
    "\n",
    "### F test \n",
    "\n",
    "Similar to ANOVA significance calculation, which also involves ratios of squared values.\n",
    "\n",
    "$$ \\hat{y} = b_0 + b_1x_1 + b_2x_2 + b_3x_3$$\n",
    "\n",
    "* $H_0 : \\hat{y} = C_0$ (All non-constant coefficients are zero)\n",
    "\n",
    "* $H_1 :$ At least one coefficient is non-zero\n",
    "\n",
    "### Total sum of squares\n",
    "\n",
    "$$ SST =\\sum_{j=1}^N{(y_j - \\bar{y})^2}$$\n",
    "\n",
    "### Regression sum of squares\n",
    "\n",
    "$$ SSR =\\sum_{j=1}^N{(\\hat{y_j} - \\bar{y})^2}$$\n",
    "\n",
    "where $\\hat{y_j}$ are model values\n",
    "\n",
    "\n",
    "### Error sum of squares\n",
    "\n",
    "$$ SSE =\\sum_{j=1}^N{(y_i - \\hat{y})^2}$$\n",
    "\n",
    "### F-statistic\n",
    "\n",
    "Mean squares statistics, calculated by dividing by degrees of freedom:\n",
    "\n",
    "* $MST =\\frac{SST}{N-1}$\n",
    "\n",
    "* $MSR =\\frac{SSR}{k}$ , where k is the number of variables\n",
    "\n",
    "* $MSE = \\frac{SSE}{N-k-1}$\n",
    "\n",
    "$F = \\frac{MSR}{MSE}$\n",
    "\n",
    "This test statistic can be compared with a critical F value, which depends on significance level $\\alpha$ and the degrees of freedom in the numerator and denominator. If F is larger, then error is small. Find F using statistical tables, or `stats.f.ppf` in Python.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49c33a77-5d96-4d5f-a517-47823f5c3669",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Multiple regression example: modeling aragonite saturation state\n",
    "\n",
    "\n",
    "![images/arag_sat.png](images/arag_sat.png)\n",
    "\n",
    "Source: Feeley et al. (2208) Evidence for upwelling of corrosive acidification water onto the continental shelf, Science\n",
    "\n",
    "At Aragonite staturation state > 1 aragonite (calcium carbonate will dissolve in seawater)\n",
    "\n",
    "Is there a way to estimate aragonite saturation state $\\Omega_{Ar}$ based on more commonly measured parameters?\n",
    "\n",
    "$$\\Omega = \\frac{[Ca^{2+}][CO_3^{2-}]}{K'_{sp}}$$\n",
    "\n",
    "Where $K'_{sp}$ is the stoichiometric solubility product function of T,S,pr and mineral phase (aragonite, calcite)\n",
    "\n",
    "$[Ca^{2+}]$  doesn't change much<br>\n",
    "$[CO_3^{2+}]$ can be calculated from chemical measurements of DIC, $pCO_2$, total alkalinity and pH (at least two of these 4 parameters).\n",
    "\n",
    "#### Models\n",
    "\n",
    "Juranek et al. (2009) describe a set of least squares regression models for aragonite saturation state, based on more commonly measured oceangraphic variables (temperature, salinity, pressure, oxygen and nitrate).\n",
    "\n",
    "Juranek, L. W., R. A. Feely, W. T. Peterson, S. R. Alin, B. Hales, K. Lee, C. L. Sabine, and J. Peterson, 2009: A novel method for determination of aragonite saturation state on the continental shelf of central Oregon using multi-parameter relationships with hydrographic data. Geophys. Res. Lett., 36, doi:10.1029/2009GL040778.\n",
    "\n",
    "##### Model 1 \n",
    "\n",
    "$$\\Omega_{arag}^e = \\beta_0 + \\beta_1T + \\beta_2S + \\beta_3P + \\beta_4O_2 + \\beta_5NO_3^-$$\n",
    "\n",
    "* Has high $R_a^2$ (\"adjusted\" $R^2$)\n",
    "* High \"variance inflation factor\"\n",
    "* Indicates multiple collinearity\n",
    "* Coefficients are ambiguous and not meaningful - When you add more data, you get will get a different answer (this is bad!)\n",
    "\n",
    "##### Adjusted $R^2$\n",
    "\n",
    "Accounts for reduction of degrees of freedom when using multiple predictor variables.\n",
    "\n",
    "$$R_a^2 = R^2 - (1-R^2)\\frac{k}{n-k-1}$$\n",
    "\n",
    "$$= 1 - \\frac{MSE}{MST}$$\n",
    "\n",
    "If the MSE is low, the adjusted R-squared is going to be high. The more observations you have, the less this adjustment matters.\n",
    "\n",
    "\n",
    "##### Variance Inflation Factor\n",
    "\n",
    "Variance Inflation Factor\n",
    "\n",
    "$$VIF  = \\frac{1}{1 - R^2}$$ \n",
    "\n",
    "where $R^2 $ from regression of predictor variables against other predictor variables. There is no clear \"cut-off\" that defines high VIF, but greater than 5 (and definitely greater than 10) is generally considered high.\n",
    "\n",
    "![images/arag_sat_table.png](images/arag_sat_table.png)\n",
    "\n",
    "##### Final Model\n",
    "\n",
    "$$\\Omega_{arag}^e = \\alpha_0 + \\alpha_1(O_2 - O_{2,r}) + \\alpha_2(T - T_r) \\times (O_2 - O_{2,r})$$\n",
    "\n",
    "* Less variables, avoids multiple collinearity\n",
    "* Includes interaction term\n",
    "* Reference values ($T_r$ and $O_{2,r}$) keep product from getting too big\n",
    "* Using variables with differing magnitudes can lead to problems like round-off errors\n",
    "* Standardizing variables (using z-scores) another common strategy\n",
    "\n",
    "![images/arag_sat_final.png](images/arag_sat_final.png)\n",
    "\n",
    "\n",
    "__Aragonite saturation state__\n",
    "\n",
    "_Red contours:_ from measured DIC and total alkalinaity\n",
    "\n",
    "_Blue contours:_ Multiple regression model\n",
    "\n",
    "__Application of the model to time series that do not have direct observations of the carbonate system parameters__\n",
    "\n",
    "![images/arag_sat_ts.png](images/arag_sat_ts.png)\n"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}