These are some notes on simple regression, multiple regression, and the general linear model.
For a more background, please see The general linear model and fMRI: Does love last forever?.
This page starts by setting up a simple regression. Then we show how the simple regression gets expressed in a design matrix. Once we have that, it's easy to extend simple regression to multiple regression. By adding some specially formed regressors, we can also express group membership, and therefore do analysis of variance. This last step is where multiple regression becomes the general linear model.
The example regression problem¶
Let's imagine that we have measured scores for a "psychopathy" personality trait in 12 students. We also have some other information about these students. For example, we measured how much sweat each student had on their palms, and we call this a "clammy" score. We first try and work out whether the "clammy" score predicts the "psychopathy" score. We'll do this with simple linear regression.
Simple linear regression¶
We first need to get our environment set up to run the code and plots we need.
# Set up notebook to show data plots in the notebook
%matplotlib inline
# Import numerical and plotting libraries
import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt
# Get t and gamma distribution code from scipy library
from scipy.stats import t as t_dist, gamma
# Get symbolic mathematics routines
from sympy import (Matrix, Symbol, Eq, MatAdd, MatMul, ones,
init_printing)
# Set up notebook to output symbolic math expressions nicely
init_printing(use_latex=True)
Here are our scores of "psychopathy" from the 12 students:
psychopathy = [11.416, 4.514, 12.204, 14.835,
8.416, 6.563, 17.343, 13.02,
15.19 , 11.902, 22.721, 22.324]
These are the skin-conductance scores to get a measure of clamminess for the handshakes of each student:
clammy = [0.389, 0.2 , 0.241, 0.463,
4.585, 1.097, 1.642, 4.972,
7.957, 5.585, 5.527, 6.964]
We happen to believe that there is some relationship between clammy and psychopathy. Plotting them together we get:
plt.plot(clammy, psychopathy, '+')
plt.xlabel('Clamminess of handshake')
plt.ylabel('Psychopathy score')
It looks like there may be some sort of straight line relationship. For example, if we draw some sort of straight line through the points, it looks kind of OK. The line we chose was just a guess by eye. It has intercept $10$ and slope $0.9$:
plt.plot(clammy, psychopathy, '+')
def my_line(x):
# My prediction for psychopathy given clamminess
return 10 + 0.9 * x
x_vals = [0, max(clammy)]
y_vals = [my_line(0), my_line(max(clammy))]
plt.plot(x_vals, y_vals)
plt.xlabel('Clamminess of handshake')
plt.ylabel('Psychopathy score')
plt.title('Clammy vs psychopathy with guessed line')
What does our straight line relationship mean?
We are saying that the values of psychopathy can be partly predicted by a straight line of formula 10 + clammy * 0.9.
To make this more general, let's call our psychopathy data $y$ - a vector with 12 values, one for each student. $y_1$ is the value for the first student (= 11.416) and $y_i$ is the value for student $i$ where $i \in 1 .. 12$.
Our clammy score is a predictor. Lets call the clammy scores $x$ - another vector with 12 values. $x_1$ is the value for the first student (= 0.389) and $x_i$ is the value for student $i$ where $i \in 1 .. 12$.
Our straight line model says:
$y_i \approx c + bx_i$
where $c$ is the intercept and $b$ is the slope. For the guessed line above:
$y_i \approx 10 + 0.9 x_i$
With the $\approx$ above, we are accepting that we will not succeed in explaining all the variation in our psychopathy data. We can rephrase this by saying that each observation is equal to the predicted value (from the formula above) plus some error for each observation:
$y_i = c + bx_i + e_i$
Simple regression in matrix form¶
We're going to rephrase the simple regression model in matrix form. Let's make the data and predictor and errors into vectors.
$y$ is the vector of values $y_1 .. y_{12}$:
y = Matrix(psychopathy)
y
$x$ is the vector of values $x_1 .. x_{12}$:
x = Matrix(clammy)
x
$\epsilon$ is the vector of errors $e_1 .. e_{12}$:
epsilons = [Symbol('e_{0}'.format(i)) for i in range(1, 13)]
epsilon = Matrix(epsilons)
epsilon
Now we can rephrase our model as:
$$y = c + b x + \epsilon$$Bear with with us for a little trick. If $o$ is a vector of ones, then we can rewrite the formula as:
$y = co + bx + \epsilon$
because $o_i = 1$ and so $co_i = c$
c, b = Symbol('c'), Symbol('b')
o = ones((12, 1))
# Compile expression without any simplification or reordering
rhs = MatAdd(MatMul(c, o, evaluate=False),
MatMul(b, x, evaluate=False),
epsilon)
Eq(y, rhs)
We can do the same calculation with matrix multiplication.
Call $X$ the matrix of two columns, where the first column is the column of ones ($o$ above) and the second column is $x$. Call $B$ the column vector:
$$ \left[ \begin{array}{B} c \\ b \\ \end{array} \right] $$X = o.row_join(x)
B = Matrix([c, b])
Eq(y, MatAdd(MatMul(X, B), epsilon))
In symbols:
$$y = X B + \epsilon$$We still haven't found our best fitting line. But before we go further, it might be obvious that we can easily add a new predictor here.
Multiple regression¶
Let's say we think that psychopathy increases with age. We add the student's age as another predictor:
age = Matrix([22.5, 25.3, 24.6, 21.4,
20.7, 23.3, 23.8, 21.7,
21.3, 25.2, 24.6, 21.8])
age
Now rename the clammy predictor from $x$ to $x_1$. Call the age predictor $x_2$. Call the slope for clammy $b_1$ (slope for $x_1$). Call the slope for age $b_2$ (slope for $x_2$). Our model is:
In this model $X$ has three columns (ones, $x_1$, and $x_2$), and the $B$ vector has three values $c, b_1, b_2$. This gives the same matrix formulation, with our new $X$ and $B$: $y = X B + \epsilon$.
This is a linear model because our model says that the data $y_i$ comes from the sum of some components ($c, b_1 x_{1, i}, b_2 x_{2, i}, e_i$).
We can keep doing this by adding more and more regressors. In general, a linear model with $p$ predictors looks like this:
$$y_i = b_1 x_{1, i} + b_2 x_{2, i} + ... b_p x_{p, i} + e_i$$In the case of the models above, the first predictor $x_1$ would be a column of ones, to express the intercept in the model.
Any model of the form above can still be phrased in the matrix form:
$$y = X B + \epsilon$$Solving the model with matrix algebra¶
The reason to formulate our problem this way is so we can use some basic matrix algebra to find the "best" line.
Let's assume that we want to find the line (intercept and slope) that gives the smallest "distance" between the estimated values (from the line), and the actual values (the data).
We'll define 'distance' as the squared difference of the predicted value from the actual value. These are the squared error terms $e_1^2, e_2^2 ... e_{n}^2$, where $n$ is the number of observations - 12 in our case.
Revising: our model is:
$$ y = \bf{X} B + \epsilon $$Where $y$ is the data vector $y_1, y_2, ... y_n$, $\bf{X}$ is the design matrix of shape $n, p$, $B$ is the parameter vector, $\beta_1, \beta_2 ... \beta_p$, and $\epsilon$ is the error vector giving errors for each observation $\epsilon_1, \epsilon_2 ... \epsilon_n$.
Each column of $\bf{X}$ is a regressor vector, so $\bf{X}$ can be thought of as the column concatenation of $p$ vectors $x_1, x_2 ... x_p$, where $x_1$ is the first regressor vector, and so on.
In our case, we want the vector $B$ such that the errors $\epsilon = y - \bf{X} B$ have the smallest sum of squares $\sum_{i=1}^n{e_i^2}$. $\sum_{i=1}^n{e_i^2}$ is called the residual sum of squares.
It might or might not be obvious that this also means that the design matrix $\bf{X}$ should be orthogonal to the errors - meaning that $\bf{X}^T \epsilon$ should be a vector of zeros.
If that is the case then we can multiply $y = {\bf X} \boldsymbol{\beta} + \epsilon$ through by $\bf{X}^T$:
$$ \bf{X}^T y = \bf{X}^T X B + \bf{X}^T \epsilon $$The last term now disappears because it is zero and:
$$ \bf{X}^T y = \bf{X}^T \bf{X} B $$If $\bf{X}^T \bf{X}$ is invertible then there is a unique solution:
$$ B = (\bf{X}^T \bf{X})^{-1} \bf{X} y $$It turns out that, if $\bf{X}^T \bf{X}$ is not invertible, then are an infinite number of solutions, and we have to choose one solution, taking into account that the parameters $B$ will depend on which solution we chose. The pseudoinverse operator gives us one particular solution. If $\bf{A}^-$ is the pseudoinverse of matrix $\bf{A}$ then the general solution for $B$, even when $\bf{X}^T \bf{X}$ is not invertible, is:
$$ B = \bf{X}^-y $$So what line do we estimate for psychopathy and clammy?
X = np.column_stack((np.ones(12), clammy))
X
B = npl.pinv(X).dot(psychopathy)
B
plt.plot(clammy, psychopathy, '+')
def my_best_line(x):
# Best prediction for psychopathy given clamminess
return B[0] + B[1] * x
x_vals = [0, max(clammy)]
y_vals = [my_best_line(0), my_best_line(max(clammy))]
plt.plot(x_vals, y_vals)
plt.xlabel('Clamminess of handshake')
plt.ylabel('Psychopathy score')
Our claim was that this slope and intercept minimize the sum of squared error:
fitted = X.dot(B)
errors = psychopathy - fitted
np.sum(errors **2)
Is this sum of squared errors smaller than our earlier guess of an intercept of 10 and a slope of 0.9?
fitted = X.dot([10, 0.9])
errors = psychopathy - fitted
np.sum(errors **2)
Contrasts¶
We can combine the values in the $B$ vector in different ways by using a contrast vector. A contrast vector $C$ is a vector of weights $c_1, c_2 ... c_p$ for each value in the $B$ vector. Assume that all vectors we've defined up until now are column vectors. Then a scalar value that is a linear combination of the $B$ values can be written:
$$ C^T B $$We now return to our original question, whether clamminess of handshake predicts psychopathy score.
If clamminess does predict psychopathy, then we would expect the slope of the best fit line between clammy and psychopathy would be different from zero.
In our model, we have two predictors, the column of ones and clammy. $p = 2$ and $B$ is length 2. We could choose just the second of the values in $B$ (therefore $b_1$) with a contrast:
This is the slope relating clammy to psychopathy. Now we might be interested if our estimate of this slope is different from zero.
To test whether the estimate is different from zero, we can divide the estimate by the variability of the estimate. This gives us an idea of how far the estimate is from zero, in terms of the variability of the estimate. We won't go into the estimate of the variability here though, we'll just assume it (the formula is in the code below). The estimate divided by the variability of the estimate gives us a t statistic.
With that introduction, here's how to do the estimation and a t statistic given the data $y$, the design $\bf{X}$, and a contrast vector $C$.
def t_stat(Y, X, C):
""" betas, t statistic and significance test given data, design matrix, contrast
This is OLS estimation; we assume the errors to have independent
and identical normal distributions around zero for each $i$ in
$\epsilon_i$ (i.i.d).
"""
# Make sure X, Y, C are all arrays
Y = np.asarray(Y)
X = np.asarray(X)
C = np.atleast_2d(C)
# Calculate the parameters
B = npl.pinv(X).dot(Y)
# The fitted values
fitted = X.dot(B)
# Residual sum of squares
RSS = ((Y - fitted)**2).sum(axis=0)
# Degrees for freedom is the number of observations n minus the number
# of independent regressors we have used. If all the regressor columns
# in X are independent then the (matrix rank of X) == p
# (where p the number of columns in X). If there is one column that can
# expressed as a linear sum of of the other columns then (matrix rank of X)
# will be p - 1 - and so on.
df = X.shape[0] - npl.matrix_rank(X)
# Mean residual sum of squares
MRSS = RSS / df
# calculate bottom half of t statistic
SE = np.sqrt(MRSS * C.dot(npl.pinv(X.T.dot(X)).dot(C.T)))
t = C.dot(B)/SE
# Get p value for t value using t distribution function
ltp = t_dist(df).cdf(t) # lower tail p
p = 1 - ltp # upper tail p
return B, t, df, p
So, does clammy predict psychopathy? If it does not, then our estimate of the slope will not be convincingly different from 0. The t test divides our estimate of the slope by the error in the estimate; large values mean that the slope is large compared to the error in the estimate.
X = np.column_stack((np.ones(12), clammy))
Y = np.asarray(psychopathy)
B, t, df, p = t_stat(Y, X, [0, 1])
t, p
Dummy coding and the general linear model¶
So far we have been doing multiple regression. That is, all the columns (except the column of ones) are continuous vectors of numbers predicting our outcome data psychopathy. These type of predictors are often called covariates.
It turns out we can use this same framework to express the fact that different observations come from different groups.
That's the same as saying that we can express analysis of variance designs using this same notation.
To do this, we use columns of dummy variables.
Let's say we get some new and interesting information. The first 4 students come from Berkeley, the second set of 4 come from Stanford, and the last set of 4 come from MIT. Maybe the student's college predicts if they are a psychopath?
How do we express this information? Let's forget about the clamminess score for now and just use the school information. Our model might be that we can can best predict the psychopathy scores by taking a representative psychopathy score for the relevant school:
We can code this with predictors in our design using indicator variables. The "Berkeley" indicator variable vector is 1 when the student is from Berkeley and zero otherwise. Similarly for the other two schools:
berkeley_indicator = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
stanford_indicator = [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0]
mit_indicator = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]
X = np.column_stack((berkeley_indicator,
stanford_indicator,
mit_indicator))
# Show the matrix nicely with LaTeX markup
Matrix(X)
These indicator columns are dummy variables where the values code for the group membership.
Now the $B$ vector will be:
$$ \left[ \begin{array}{B} \mu_{Berkeley} \\ \mu_{Stanford} \\ \mu_{MIT} \\ \end{array} \right] $$When we estimate these using the least squares method, what estimates will we get for $B$?
B = npl.pinv(X).dot(psychopathy)
B
print(np.mean(psychopathy[:4]))
print(np.mean(psychopathy[4:8]))
print(np.mean(psychopathy[8:]))
It looks like the MIT students are a bit more psychopathic. Are they more psychopathic than Berkeley and Stanford?
We can test whether the mean for the MIT students is greater than the mean of (mean for Berkeley, mean for Stanford):
B, t, df, p = t_stat(psychopathy, X, [-0.5, -0.5, 1])
t, p
Ah - yes - just as we suspected.
The model above expresses the effect of group membersip. It is the expression of an analysis of variance (ANOVA) model using $y = XB + \epsilon$.
ANCOVA in the General Linear Model¶
Our formulation $y = XB + \epsilon$ makes it very easy to add extra regressors to models with group membership. For example, here is a simple ANCOVA model (analysis of covariance).
ANCOVA is a specific term for the case where we have a model with both group membership (ANOVA model) and one or more continuous covariate.
For example, we can add back our clamminess score to the mix. Does it explain anything once we know which school the student is at?
X = np.column_stack((berkeley_indicator,
stanford_indicator,
mit_indicator,
clammy))
X
No - it looks like the MIT students also have clammy hands, and once we know that the student is from MIT, the clammy score is not as useful.
Displaying the design matrix as an image¶
We can show the design as an image, by scaling the values with columns.
We scale within columns because we care more about the seeing variation within the regressor than between regressors. For example, if we have a regressor varying between 0 and 1, and another between 0 and 1000, without scaling, the column with the larger numbers will swamp the variation in the column with the smaller numbers.
def scale_design_mtx(X):
"""utility to scale the design matrix for display
This scales the columns to their own range so we can see the variations
across the column for all the columns, regardless of the scaling of the
column.
"""
mi, ma = X.min(axis=0), X.max(axis=0)
col_neq = (ma - mi) > 1.e-8
Xs = np.ones_like(X)
mi = mi[col_neq]
ma = ma[col_neq]
Xs[:,col_neq] = (X[:,col_neq] - mi)/(ma - mi)
return Xs
Then we can display this scaled design with a title and some default image display parameters:
def show_design(X, design_title):
""" Show the design matrix nicely """
plt.imshow(scale_design_mtx(X),
interpolation='nearest',
cmap='gray') # Gray colormap
plt.title(design_title)
We can then see our ANCOVA design above at a glance:
show_design(X, 'ANCOVA')
