DSV Lecture 20
2014-04-01
James Bagrow, james.bagrow@uvm.edu, http://bagrow.com
Linear regression
We discussed linear regression way back in Lecture 9.
- Why are we talking about it again?
First, let's recap:
The goal with linear regression was to find the line that best fit some given XY-data (and test if that fit was significant). For the \(i\)th data point:
\[y_i = \beta_0 + \beta_1 x_i + \epsilon_i\]
We want to find the pair of values (\(\beta_0\),\(\beta_1\)) that minimize the sum of square errors \(S\):
\[S(\beta) = \sum_i \epsilon_i^2 = \sum_i \left[y_i - \left(\beta_0 + \beta_1 x_i\right)\right]^2\]
Of course, this holds for all the data points together:
\[y_1 = \beta_0 + \beta_1 x_1 + \epsilon_1\\
y_2 = \beta_0 + \beta_1 x_2 + \epsilon_2\\
y_3 = \beta_0 + \beta_1 x_3 + \epsilon_3\\
\vdots \\
y_n = \beta_0 + \beta_1 x_n + \epsilon_n
\]
We can write this in matrix notation:
Now a single matrix equation captures all individual linear regression equations. Furthermore, if we want to estimate the \(\mathbf{\beta}\)'s with Ordinary Least Squares (OLS), we have a (relatively) simple way to write down the best \(\mathbf{beta}\)'s.
In matrix form the residuals are:
\[\mathbf{\epsilon} = \mathbf{y} - \mathbf{X}\mathbf{\beta}\]
The sum of the squared residuals is then (using the inner product):
\[\mathbf{\epsilon}^\mathrm{T}\mathbf{\epsilon} = \left(\mathbf{y} - \mathbf{X}\mathbf{\beta}\right)^\mathrm{T}\left(\mathbf{y} - \mathbf{X}\mathbf{\beta}\right)\]
where \(\mathbf{\epsilon}^\mathrm{T}\) is the transpose of \(\mathbf{\epsilon}\).
With a little calculus and elbow grease we can find the minimum of this equation, which can be solved to give us the best estimates for \(\mathbf{\beta}\), called \(\mathbf{\hat{\beta}}\):
\[\mathbf{\hat{\beta}} = \left(\mathbf{X}^\mathrm{T}\mathbf{X}\right)^{-1}\mathbf{X}^\mathrm{T} \mathbf{y}\]
assuming \(\mathbf{X}^\mathrm{T}\mathbf{X}\) is nonsingular, meaning we can compute the matrix inverse.
We already computed the values of \(\beta_0\) and \(\beta_1\), so why are we now getting into all this matrix nonsense?
Because the above derivation is not limited to two coefficients and a linear regression of the form \(y = \beta_0 + \beta_1 x\).
Until now, we've been doing simple linear regression.
Multiple linear regression.
There is no reason why we need to stop our equation at two coefficients. What about this:
\[\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} =
\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\
1 & x_{21} & x_{22} & \cdots & x_{2p} \\
1 & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n1} & x_{n2} & \cdots & x_{np}
\end{bmatrix}
\begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix} +
\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}
\]
In matrix form that is exactly the same equation:
\[\mathbf{y} = \mathbf{X}\mathbf{\beta} + \mathbf{\epsilon}\]
except that \(\mathbf{X}\) has \(p+1\) columns and \(\mathbf{\beta}\) has \(p+1\) rows.
- This means we can solve exactly the same linear regression but we are not limited to a single column of \(x\) data!
*** For example, we may have reason to study a linear model relating a person's height and weight:
\[\mathrm{weight} = \beta_0 + \beta_1 \mathrm{height}\]
Now we can also incorporate (probably unrealistically) age data as well:
\[\mathrm{weight} = \beta_0 + \beta_1 \mathrm{height} + \beta_2 \mathrm{age}\]
We have a new coefficient, and we have a new column (the age column) in our data matrix. But estimating \(\mathbf{\hat{\beta}}\) remains the same. ***
A note on jargon
There's a lot of names for the quantities in these equations, it's hard to keep straight.
The \(y_i\)'s are called:
- Regressand
- endogenous variable
- response variable
- dependent variable
The \(x_i\)'s (including the constant \(1\)) are called:
- Regressors
- exogenous variables
- explanatory variables
- independent variables
The matrix of data \(\mathbf{X}\) is called the design matrix.
The \(\beta\)'s are called:
- effects
- regression coefficients
- regression parameters
The vector \(\mathbf{\beta}\) is called the paramter or coefficient vector.
The \(\epsilon\)'s are called errors, error term, or noise.
Handy mnemonic for the \(x\)'s and \(y\)'s:
The \(x\)'s are the exogenous variable because "exogenous" contains "x".
Let's see this in action:
We're going to use a new library called statsmodels.
First let's generate some data:
Let's even print the first few lines:
Looks good.
If you're comfortable with matrix operations, you can also build the data using matrix multiplication:
But matrix operations are not very comfortable in python (at least compared to MATLAB), so you might just want to avoid them. It's up to you.
Anyway, now that we've got the data generated, let's use statsmodels OLS function to find the \(\mathbf{\hat{\beta}}\):
Wow, that was pretty easy, but we now we need to understand what results is...
model.fit returns an object that contains a bunch of data about the fit, residuals, coefficient values, and other useful statistics. It also has a few methods attached to it. One very handy method it has is summary:
This type of display is very common across statistical packages. There's three main panels.
The top panel contains information about the size of the data and the overall accuracy of the fit, as well as some nice records like time of the calculation.
The middle panel contains information about each coefficient of the fit.
The bottom panel contains information about the residuals \(\epsilon_i\) of the fit, are they normally distributed (omnibus test), do they possess serial correlations (Durbin-Watson), etc.
We can also access the information shown in the summary directly:
Here are all the things result has:
There's lots of data we can interact with regarding the fit...
Here's another example. In this case we are going to use sm.add_constant to more quickly add the column of ones to the design matrix.
- This uses real data as well so we don't need to add fake \(\epsilon\)'s ourselves.
- But let's add a random variable and see if the model rejects it...
Let's make a quickie plot of the solution and the data:
OK, so we see that, since the "Nonsense" axis isn't meaningful, our regression is essentially a flat plane. Let's do an example with real data, where both variables matter:
This data measures the time to recovery from anaesthetic as a function of blood pressure and (the logarithm) of the dosage.
Fit it up!!!