Primer: the generalised linear model

On By @thinquanautIn Primers

This is the first of three short primers on methods that keep turning up in the lessons, written for someone who has done lesson 2 and lesson 3 and wants the idea behind the tools. Each one is introductory, runs on the same real data as the lessons, and takes about twenty minutes.

A generalised linear model, or GLM, is the thing you reach for when you want to know what a neuron responds to and there is more than one candidate. The tuning curve in lesson 2 answered “how does firing depend on contrast?” by averaging. A GLM answers “how does firing depend on contrast, and on what the mouse did, and on how long the experiment has been running, all at once, with each effect measured while the others are held fixed?” It is the workhorse of systems neuroscience, and the good news is that you have already fitted one: logistic regression in lesson 3 is a GLM.

One recipe, three dishes

Every GLM has the same three parts. A weighted sum of the predictors, exactly as in linear regression: so many points per unit of contrast, so many for turning the wheel, plus a baseline. A link, a fixed function that bends that sum into the range the data can actually occupy. And a noise model, a statement of how the data scatter around the prediction. Change the link and the noise model and you get a different member of the family.

The same recipe three times, on made-up data. Left: no link, Gaussian noise; the familiar straight line. Middle: a log-odds link squashes the line into a probability between 0 and 1, and the data are yes/no outcomes; this is logistic regression, lesson 3. Right: a log link makes the prediction an always-positive rate, and the data are counts that scatter more when the rate is higher; this is Poisson regression, the one neuroscience uses most.
The same recipe three times, on made-up data. Left: no link, Gaussian noise; the familiar straight line. Middle: a log-odds link squashes the line into a probability between 0 and 1, and the data are yes/no outcomes; this is logistic regression, lesson 3. Right: a log link makes the prediction an always-positive rate, and the data are counts that scatter more when the rate is higher; this is Poisson regression, the one neuroscience uses most.

Spike counts are the right-hand panel. They cannot be negative, so the prediction should not be either, which the log link guarantees: the model predicts log(rate), and the rate is exp of that. And their noise is not constant. Look at neuron 141 from lesson 2, counting spikes in the 200 ms after the stimulus:

import numpy as np
import matplotlib.pyplot as plt

dat = np.load("steinmetz_session11.zip", allow_pickle=True)
spks = dat["spks"]
contrast = dat["contrast_right"]
contrast_left = dat["contrast_left"]
response = dat["response"]                 # -1 wheel left, 0 no move, 1 wheel right
dt = float(dat["bin_size"])
t = np.arange(spks.shape[2]) * dt - 0.5
post = (t >= 0.05) & (t < 0.25)

y = spks[141][:, post].sum(axis=1)         # neuron 141's spike count on each trial, the thing we will model
levels = np.unique(contrast)
for c in levels:
    print(f"contrast {c}: mean count {y[contrast == c].mean():5.2f}, variance {y[contrast == c].var():5.2f}")
# contrast 0.0: mean count  1.53, variance  4.13
# contrast 0.25: mean count  7.06, variance 13.66
# contrast 0.5: mean count 10.10, variance 24.17
# contrast 1.0: mean count 13.26, variance 21.34
Neuron 141's spike-count variance against its mean, one point per contrast level. The dashed line is what a Poisson process would give: variance equal to the mean. Real neurons are a bit noisier than that, but the pattern is unmistakable: the more it fires, the more it varies.
Neuron 141’s spike-count variance against its mean, one point per contrast level. The dashed line is what a Poisson process would give: variance equal to the mean. Real neurons are a bit noisier than that, but the pattern is unmistakable: the more it fires, the more it varies.

A straight-line fit with constant noise would treat a two-spike miss on a blank trial and a two-spike miss on a full-contrast trial as equally surprising. They are not. The Poisson model knows that, and that is most of why it fits spike data better.

Fitting one in eight lines

Fitting means choosing the weights that make the observed counts most probable. For a Poisson GLM the gradient of that probability has a form so simple it is worth remembering: for each predictor, add up (observed count minus predicted count) times the predictor, over trials. Walk uphill along that gradient and you arrive at the best weights. This is the same loop as lesson 3’s decoder with exp in place of the sigmoid.

def poisson_glm(X, y, lr=0.05, steps=20000):
    """Fit log(rate) = b0 + X @ b by gradient ascent on the Poisson log-likelihood."""
    X1 = np.column_stack([np.ones(len(y)), X])        # a column of ones for the intercept
    w = np.zeros(X1.shape[1])
    for _ in range(steps):
        rate = np.exp(X1 @ w)                          # the inverse link: weights to a positive rate
        w += lr * X1.T @ (y - rate) / len(y)           # the gradient is just (observed - predicted) times X
    return w

def predict(X, w):
    return np.exp(np.column_stack([np.ones(len(X)), X]) @ w)

X_lin = contrast[:, None]                             # one predictor: contrast as a number
w_lin = poisson_glm(X_lin, y)
print("log rate =", w_lin.round(2), "-> predicted counts:", predict(levels[:, None], w_lin).round(2))

X_cat = np.column_stack([(contrast == c) for c in levels[1:]]).astype(float)   # one column per non-zero contrast
w_cat = poisson_glm(X_cat, y)
print("one-hot ->", predict(np.vstack([np.zeros(3), np.eye(3)]), w_cat).round(2))
# log rate = [1.   1.71] -> predicted counts: [ 2.73  4.18  6.4  15.02]
# one-hot -> [ 1.53  7.06 10.1  13.26]

Two fits of the same neuron. The first uses contrast as a single number and predicts 2.7, 4.2, 6.4 and 15 spikes at the four levels. The real means are 1.5, 7.1, 10.1 and 13.3. The fit is poor, and the reason is instructive: with one weight on contrast, the log link forces the rate to grow exponentially with contrast, and this neuron saturates instead. The second fit gives each contrast level its own weight, and recovers the four means exactly. That is not a coincidence: a tuning curve is a GLM with one weight per condition.

Neuron 141's contrast tuning (black, with standard errors) and the two GLM fits. One weight on contrast as a number forces an exponential shape that misses every point. One weight per contrast level reproduces the tuning curve from lesson 2 exactly.
Neuron 141’s contrast tuning (black, with standard errors) and the two GLM fits. One weight on contrast as a number forces an exponential shape that misses every point. One weight per contrast level reproduces the tuning curve from lesson 2 exactly.

The lesson from the orange curve generalises: the link function imposes a shape on how a numeric predictor acts, and if the real relationship has a different shape you must give the model room, with categories, with a square term, or with a transform of the predictor. A GLM is only as good as the predictors you hand it.

Now add everything else

Here is what averaging cannot do. The mouse in this experiment also saw a stimulus on the left on some trials, turned the wheel left or right or not at all, and got tired or practised as the session went on. Any of those might move neuron 141. Put them all in as columns and fit once.

def log_likelihood(y, rate):
    """How probable the observed counts are under predicted Poisson rates (up to a constant)."""
    return np.sum(y * np.log(rate) - rate)

X_full = np.column_stack([
    X_cat,                                   # contrast on the right, one-hot
    contrast_left > 0,                       # was there anything on the left?
    response == -1,                          # did the mouse turn the wheel left?
    response == 1,                           # or right?
    np.arange(len(y)) / len(y),              # how far into the session (0 to 1)
]).astype(float)
names = ["contrast 0.25", "contrast 0.5", "contrast 1.0", "left stimulus", "turned left", "turned right", "time in session"]
w_full = poisson_glm(X_full, y)
for name, b in zip(names, w_full[1:]):
    print(f"{name:16s} x{np.exp(b):.2f}")
print(f"baseline {np.exp(w_full[0]):.2f} spikes")
print("log-likelihood, contrast only:", log_likelihood(y, predict(X_cat, w_cat)).round(1))
print("log-likelihood, full model:    ", log_likelihood(y, predict(X_full, w_full)).round(1))
# contrast 0.25    x3.77
# contrast 0.5     x5.20
# contrast 1.0     x7.24
# left stimulus    x1.12
# turned left      x1.76
# turned right     x1.52
# time in session  x2.24
# baseline 0.69 spikes
# log-likelihood, contrast only: 2386.7
# log-likelihood, full model:     2449.1
Every effect in the full model, expressed as a multiplier on firing because the link is a log: a weight of 0.57 becomes x1.76. Contrast dominates, turning the wheel in either direction adds half as much again, a stimulus on the left does almost nothing, and the neuron fires more than twice as much at the end of the session as at the start.
Every effect in the full model, expressed as a multiplier on firing because the link is a log: a weight of 0.57 becomes x1.76. Contrast dominates, turning the wheel in either direction adds half as much again, a stimulus on the left does almost nothing, and the neuron fires more than twice as much at the end of the session as at the start.

Read the multipliers. A full-contrast stimulus multiplies this neuron’s firing by about seven. A stimulus on the left screen, in the other visual field, does almost nothing, which is what a visual cortex neuron with a receptive field on the right should do. Turning the wheel multiplies firing by 1.5 to 1.8 regardless of direction; this is the movement signal that lesson 3 warned about and lesson 2 saw as the second hump, now measured and separated from the stimulus. And the last line is the surprise: the predictor “time in session” multiplies firing by 2.2, meaning the neuron fires more than twice as much at the end of the hour as at the start. No tuning curve would have shown that, because tuning curves average over time. The log-likelihood, our measure of how probable the data are under the model, improves by 62 units, which is a large amount.

rate_full = predict(X_full, w_full)
fig, ax = plt.subplots(figsize=(10, 3.8))
ax.plot(y, lw=0.8, label="observed count")
ax.plot(rate_full, lw=1.2, label="predicted by the full model")
ax.set(xlabel="trial", ylabel="spikes in the window", xlim=(0, 120))
ax.legend()
plt.show()
The first 120 trials: observed counts, the contrast-only prediction, and the full model's prediction. The full model follows the trial-to-trial swings that contrast alone cannot, and climbs slowly as the session goes on.
The first 120 trials: observed counts, the contrast-only prediction, and the full model’s prediction. The full model follows the trial-to-trial swings that contrast alone cannot, and climbs slowly as the session goes on.

What a GLM is for, and what it is not

  • It separates effects that co-occur. Stimulus and movement are correlated in this task; the GLM gives each its own weight while accounting for the other. This is the main reason the field uses it.
  • Its weights have units and meaning. With a log link, exp(weight) is a multiplier on firing, which is a sentence you can put in a paper.
  • It is still linear in the weights. The log link does not let it discover that the contrast response saturates; you have to give it categories. Interactions (does movement matter more at high contrast?) also have to be added by hand as a product column.
  • Real neurons are over-dispersed. Their variance is a bit more than their mean, as the second figure showed. The weights are still fine; the confidence you put on them should be a little looser than textbook Poisson theory says.
  • Correlated predictors are its weakness. If two columns always move together, the model cannot tell which one matters, and the weights become unstable. Check your design before trusting the weights.

Exercises

  1. Fit the full model to a neuron in MOs (secondary motor cortex) or MD (thalamus) instead of 141. Which multiplier is largest there?
  2. Add an interaction: a column equal to (contrast == 1.0) * (response != 0). Does movement matter more on high-contrast trials?
  3. Replace the Poisson fit with ordinary least squares (np.linalg.lstsq on the same columns) and compare the predictions on blank trials. Where does the linear model go wrong?
  4. Harder: fit a separate GLM to every VISp neuron and plot the distribution of the “turned left” multiplier across the population. How many visual neurons carry a movement signal?

Next primer: independent component analysis. Data: Steinmetz et al., Nature 2019, CC-BY 4.0, via Neuromatch Academy.

Keep reading