CS109B Data Science 2: Advanced Topics in Data Science

Lecture 3 - Coding Environment Setup and review of statsmodels¶

Notebook B¶

Harvard University
Spring 2021
Instructors: Mark Glickman, Pavlos Protopapas, and Chris Tanner
Additional Instructor: Eleni Kaxiras

Content: Eleni Kaxiras and Will Claybaugh

---

Learning Goals¶

By the end of this lab, you should be able to:

• use np.linalg.vander
• use the weird R-style formulas in statsmodels
• practice least-squares regression in statsmodels

Basis Functions¶

In our models we can use various types of functions as basis functions. Strictly speaking, in linear algebra where a basis for a subspace S of $\mathbb{R}^n$ is a set of vectors that spans S and is linearly independent. As a reminder, a set of vectors $\textbf{v}_1, \textbf{v}_2, ..., \textbf{v}_k$ are considered linearly independent if they cannot be written as a linear combination of each other, such that, if: $c_1\textbf{v}_1+c_2\textbf{v}_2+ ...+ c_k\textbf{v}_k = \textbf{0}$ then $c_1,c_2,...,c_k$ are all zero.

In data science where we have lots of imperfect data (with errors), as well as imperfect computers (with round-off errors), when we substitute their values into the matrices we almost always get column degeneracy, meaning, some of our columns become linear combinations of each other. Especially so if we use the monomial basis and go beyond ~5,6 degree of the polynomial.

Examples are:

• Monomials such as $x,x^2,x^4,x^5$
• Sigmoid/ReLU functions (neural networks)
• Fourier functions
• Wavelets
• Splines

The matrix produced when we substitute the values of our data into the basis functions is called the design matrix.

Linear/Polynomial Regression¶

We will use the diabetes dataset.

Variables are:

• subject: subject ID number
• age: age diagnosed with diabetes
• acidity: a measure of acidity called base deficit Response:
• y: natural log of serum C-peptide concentration

Original source is Sockett et al. (1987) mentioned in Hastie and Tibshirani's book "Generalized Additive Models".

Reading the data in Pandas:

Create the design matrix for a fictitious dataset¶

Let's keep just the age feature and create some columns of our own. Let's see how good this matrix is before we create the design matrix.

Let's check if the columns of A are linearly independent by using some linear algebra methods from numpy.linalg and sympy.

Create the design matrix for age using a polynomial basis¶

Let's keep just the age feature again and create the design matrix using a polynomial of degree n. First we will use the basic numpy formula vander().

Change the degree of the polynomial and comment on what happens to the condition and rank of the matrix.

Linear/Polynomial regression with statsmodels.¶

As you remember from 109a, we have two tools for Linear Regression:

• statsmodels https://www.statsmodels.org/stable/regression.html, and
• sklearnhttps://scikit-learn.org/stable/index.html

Previously, in this notebook, we worked from a vector of target values and a design matrix we built ourself. In 109a we used e.g. sklearn's PolynomialFeatures to build the matrix. Now we will look at statsmodels which allows users to fit statistical models using R-style formulas. They build the target value and design matrix for you.

Note: Categorical features (e.g. let's say we had a categorical feature called Region, are designated by C(Region)), polynomial features (e.g. age) are entered as np.power(age, n) where n is the degree of the polynomial OR np.vander(age, n, increasing=True).

# Example: if our target variable is 'Lottery', while 'Region' is a categorical predictor and all the others are numerical:
df = dta.data[['Lottery', 'Literacy', 'Wealth', 'Region']]

formula='Lottery ~ Literacy + Wealth + C(Region) + Literacy * Wealth'

For more on these formulas see:

• https://www.statsmodels.org/stable/examples/notebooks/generated/formulas.html
• https://patsy.readthedocs.io/en/latest/overview.html

Let's build a dataframe to predict values on (sometimes this is just the test or validation set). Very useful for making pretty plots of the model predictions - predict for TONS of values, not just whatever's in the training set.

Use get_prediction( ).summary_frame() to get the model's prediction (and error bars!)

Plot the data, the fitted model, the confidence intervals, and the prediction intervals. For more on how statsmodels calculates these intervals see: https://www.statsmodels.org/stable/_modules/statsmodels/regression/_prediction.html

• Fit a 3rd degree polynomial model to predict y using only age and
• Plot the model and its confidence intervals.
• Change the degree of your polynomial and see what happens to the fitted curve.
• Does our model have an intercept? Note: we can discover the existence or not of an intercept in our model by running:

  model_name.params

It's easier to build higher order polynomials using np.vandrer().

QR decomposition (Beyond the scope of this class)¶

As you know, to find the parameters of our model, we may try to solve the so-called normal equations, which, written in matrix form, are:
\begin{equation} (\boldsymbol{A^T}\cdot \boldsymbol{A}) \cdot \boldsymbol{b} = \boldsymbol{A} \cdot \boldsymbol{y} \end{equation}

The direct solution is $\hat{\boldsymbol{b}}=(\boldsymbol{A}^T\cdot \boldsymbol{A})^{-1}\cdot \boldsymbol{A}^T \cdot \boldsymbol{y}$

Solving the least-squares problem directly via the normal equations is susceptible to roundoff error when the condition of the matrix $\boldsymbol{A}$ is large. An alternative technique involves QR decomposition (details in any good linear algebra book). statsmodels lets you use this technique via a parameter in the .fit:

   .fit(method='qr')

Let's try with QR now¶