# load libraries

library(data.table)
library(ggplot2)

# clean work space

rm(list = ls())

# init colorscheme

COL <- c("#2271B2", "#E69F00", "#D55E00")
names(COL) <- c("blue", "orange", "red")
theme_set(
  theme_minimal(base_size = 13) +
    theme(
      panel.grid.minor = element_blank(),
      strip.text = element_text(face = "bold"),
      legend.position = "bottom"
    )
)
update_geom_defaults("point", list(size = 2))
update_geom_defaults("line", list(linewidth = 0.8))

Overview

So far, you have worked with simple linear regression, where one predictor is used to explain one outcome. Multiple regression extends that framework by including two or more predictors at the same time.

In this tutorial, you will use a small V1 tuning dataset. The example is built around a familiar receptive-field idea: a V1 neuron tends to respond more strongly to stimuli near its preferred orientation and preferred spatial frequency.

You will model spike rate using:

  • one continuous predictor: ori_distance
  • a second continuous predictor: sf_distance

The raw stimulus variables are stim_ori and stim_sf. For the regression, we use distance from the neuron’s preferred value. A value of 0 means the stimulus matches the neuron’s preference. Larger values mean the stimulus is farther away from what this neuron prefers.

Part 1 - Introducing the data

Load the V1 tuning file.

d <- fread("data/v1_tuning_summary.csv")
str(d)
## Classes 'data.table' and 'data.frame':   160 obs. of  6 variables:
##  $ stim_ori           : num  125.8 100.2 25.2 51.4 100 ...
##  $ stim_sf            : num  0.768 2.882 12.744 0.739 0.768 ...
##  $ attention_condition: chr  "ignore" "ignore" "ignore" "ignore" ...
##  $ ori_distance       : num  35.76 10.18 64.77 38.57 9.97 ...
##  $ sf_distance        : num  2.38 0.473 1.672 2.436 2.38 ...
##  $ spike_rate         : num  32.3 41.2 20.6 28.7 31.2 ...
##  - attr(*, ".internal.selfref")=<externalptr>
head(d)
##    stim_ori stim_sf attention_condition ori_distance sf_distance spike_rate
##       <num>   <num>              <char>        <num>       <num>      <num>
## 1:   125.76   0.768              ignore        35.76       2.380      32.28
## 2:   100.18   2.882              ignore        10.18       0.473      41.16
## 3:    25.23  12.744              ignore        64.77       1.672      20.59
## 4:    51.43   0.739              ignore        38.57       2.436      28.74
## 5:    99.97   0.768              ignore         9.97       2.380      31.24
## 6:     4.52   1.857              ignore        85.48       1.107      21.76

Each row is one stimulus condition in one attention condition for one V1 neuron. The variables are:

  • stim_ori - stimulus orientation in degrees
  • stim_sf - stimulus spatial frequency in cycles per degree
  • attention_condition - whether the stimulus was attended or ignored
  • ori_distance - distance from the neuron’s preferred orientation
  • sf_distance - distance from the neuron’s preferred spatial frequency
  • spike_rate - mean firing rate in spikes per second

For this teaching dataset, the preferred orientation is 90 degrees and the preferred spatial frequency is 4 cycles per degree.

Part 2 - Visualise the tuning relationships

Before fitting a model, plot the data.

ggplot(d, aes(x = ori_distance, y = spike_rate)) +
  geom_point(colour = COL[1], alpha = 0.8) +
  labs(
    x = "Distance from preferred orientation (degrees)",
    y = "Spike rate",
    title = "Spike rate and orientation tuning"
  )
Spike rate as a function of orientation distance.

Spike rate as a function of orientation distance.

This plot asks whether the neuron fires less as the stimulus moves farther from its preferred orientation. That is the same question as a simple linear regression.

Now make the same plot for spatial frequency.

ggplot(d, aes(x = sf_distance, y = spike_rate)) +
  geom_point(colour = COL[2], alpha = 0.8) +
  labs(
    x = "Distance from preferred spatial frequency",
    y = "Spike rate",
    title = "Spike rate and spatial-frequency tuning"
  )
Spike rate as a function of spatial-frequency distance.

Spike rate as a function of spatial-frequency distance.

Multiple regression lets you test whether both stimulus properties explain spike rate when they are included in the model at the same time.

Part 3 - Simple regression baseline

Start with the simplest model: predict spike rate from orientation distance alone.

fit_simple <- lm(spike_rate ~ ori_distance, data = d)
summary(fit_simple)
## 
## Call:
## lm(formula = spike_rate ~ ori_distance, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.7994  -5.5674  -0.3523   5.8066  17.7584 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  42.26639    1.15085  36.726  < 2e-16 ***
## ori_distance -0.17574    0.02191  -8.019 2.23e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.313 on 158 degrees of freedom
## Multiple R-squared:  0.2893, Adjusted R-squared:  0.2848 
## F-statistic: 64.31 on 1 and 158 DF,  p-value: 2.233e-13

This model asks whether spike rate changes as the stimulus moves farther away from the neuron’s preferred orientation.

Plot the simple regression line.

ggplot(d, aes(x = ori_distance, y = spike_rate)) +
  geom_point(colour = COL[1], alpha = 0.55) +
  geom_smooth(method = "lm", se = FALSE, colour = COL[3]) +
  labs(
    x = "Distance from preferred orientation (degrees)",
    y = "Spike rate",
    title = "Simple regression model"
  )
Simple regression of spike rate on orientation distance.

Simple regression of spike rate on orientation distance.

Part 4 - Add a second continuous predictor

Now add sf_distance as a second continuous predictor.

fit_multiple <- lm(spike_rate ~ ori_distance + sf_distance, data = d)
summary(fit_multiple)
## 
## Call:
## lm(formula = spike_rate ~ ori_distance + sf_distance, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.0392  -2.9391   0.3859   3.0573   9.5280 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  52.39230    0.95745   54.72   <2e-16 ***
## ori_distance -0.20310    0.01372  -14.80   <2e-16 ***
## sf_distance  -7.03029    0.44231  -15.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.542 on 157 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7241 
## F-statistic: 209.7 on 2 and 157 DF,  p-value: < 2.2e-16

This is a multiple regression model because it has two predictors. The model is:

\[ \widehat{spike\_rate} = \beta_0 + \beta_1 ori\_distance + \beta_2 sf\_distance \]

The coefficient for ori_distance estimates the relationship between orientation distance and spike rate while holding sf_distance constant.

The coefficient for sf_distance estimates the relationship between spatial-frequency distance and spike rate while holding ori_distance constant.

Compare the R-squared values from the two models.

summary(fit_simple)$r.squared
## [1] 0.2892814
summary(fit_multiple)$r.squared
## [1] 0.727604

If R-squared increases, the second predictor explains additional variance above and beyond orientation distance alone.

Now plot the predictions from the multiple regression model.

d[, predicted_multiple := predict(fit_multiple)]

ggplot(d, aes(x = predicted_multiple, y = spike_rate)) +
  geom_point(colour = COL[1], alpha = 0.55) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = COL[3]) +
  labs(
    x = "Predicted spike rate",
    y = "Observed spike rate",
    title = "Multiple regression predictions"
  )
Observed spike rates and predictions from the two-predictor model.

Observed spike rates and predictions from the two-predictor model.

Points close to the dashed line are observations that the model predicts well.

A second way to visualise the model is to draw predicted lines for a few different values of sf_distance.

new_dat <- CJ(
  ori_distance = seq(min(d$ori_distance), max(d$ori_distance), length.out = 100),
  sf_distance = c(0.5, 1.5, 2.5)
)

new_dat[, predicted_spike_rate := predict(fit_multiple, newdata = new_dat)]
new_dat[, sf_label := paste("sf_distance =", sf_distance)]

ggplot(d, aes(x = ori_distance, y = spike_rate)) +
  geom_point(colour = "grey60", alpha = 0.35) +
  geom_line(
    data = new_dat,
    aes(y = predicted_spike_rate, colour = sf_label),
    linewidth = 1
  ) +
  labs(
    x = "Distance from preferred orientation (degrees)",
    y = "Spike rate",
    colour = "Spatial frequency",
    title = "Two-predictor regression model"
  ) +
  scale_colour_manual(values = COL)
Predicted orientation-distance lines at three spatial-frequency distances.

Predicted orientation-distance lines at three spatial-frequency distances.

The lines are separated because the model also uses sf_distance to predict spike rate.

Part 5 - Model comparison

You can compare the simple model with the two-predictor model using anova(). This asks whether the model improves after adding sf_distance.

anova(fit_simple, fit_multiple)
## Analysis of Variance Table
## 
## Model 1: spike_rate ~ ori_distance
## Model 2: spike_rate ~ ori_distance + sf_distance
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    158 8449.1                                  
## 2    157 3238.3  1    5210.8 252.63 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This table asks whether adding sf_distance improves fit beyond the simple regression model with ori_distance alone.

Adjusted R-squared is also useful in multiple regression because it penalises models for adding predictors.

summary(fit_simple)$adj.r.squared
## [1] 0.2847832
summary(fit_multiple)$adj.r.squared
## [1] 0.724134

Part 6 - Multiple regression assumptions

Because each row is one stimulus condition, the independence assumption is reasonable for this simplified teaching example. The remaining assumptions are the usual regression checks:

  1. linearity
  2. homoscedasticity
  3. approximately normal residuals

Multiple regression also adds a practical concern: the predictors should not be so strongly correlated with one another that it becomes difficult to separate their individual contributions.

cor(d$ori_distance, d$sf_distance)
## [1] -0.1254647

A high correlation between predictors is a warning sign for multicollinearity. It does not automatically make a model invalid, but it does make individual coefficient estimates harder to interpret.

d_diag <- data.table(fitted = fitted(fit_multiple), resid = resid(fit_multiple))

ggplot(d_diag, aes(x = fitted, y = resid)) +
  geom_point(colour = COL[1], alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    x = "Fitted values",
    y = "Residuals",
    title = "Residuals vs fitted"
  )

ggplot(d_diag, aes(sample = resid)) +
  stat_qq(colour = COL[1]) +
  stat_qq_line(colour = COL[2]) +
  labs(
    x = "Theoretical quantiles",
    y = "Sample quantiles",
    title = "Q-Q plot of residuals"
  )

Look for residuals scattered around zero and Q-Q points that stay reasonably close to the diagonal.

Part 7 - Add a categorical predictor

Now add attention_condition as a categorical predictor.

In R, categorical predictors need to be represented with numbers before they can go into a regression model. For a two-level factor, R does this by creating a 0/1 comparison. Here, ignore is the reference level:

  • ignore is coded as 0
  • attend is coded as 1
d[, attention_fac := factor(attention_condition, levels = c("ignore", "attend"))]

contrasts(d$attention_fac)
##        attend
## ignore      0
## attend      1

fit_group <- lm(spike_rate ~ ori_distance + sf_distance + attention_fac, data = d)
summary(fit_group)
## 
## Call:
## lm(formula = spike_rate ~ ori_distance + sf_distance + attention_fac, 
##     data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.2978  -2.5522  -0.4035   2.6445   8.4309 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         48.71253    0.89906  54.182  < 2e-16 ***
## ori_distance        -0.20579    0.01134 -18.153  < 2e-16 ***
## sf_distance         -6.12729    0.38014 -16.119  < 2e-16 ***
## attention_facattend  5.32197    0.61832   8.607 7.62e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.752 on 156 degrees of freedom
## Multiple R-squared:  0.8153, Adjusted R-squared:  0.8118 
## F-statistic: 229.6 on 3 and 156 DF,  p-value: < 2.2e-16

The coefficient for attention_facattend is the estimated difference between the attended and ignored conditions, after accounting for orientation distance and spatial-frequency distance.

If the coefficient is positive, the model predicts higher spike rates in the attended condition than in the ignored condition. If it is negative, the model predicts lower spike rates in the attended condition.

Compare this model with the two-continuous-predictor model:

anova(fit_multiple, fit_group)
## Analysis of Variance Table
## 
## Model 1: spike_rate ~ ori_distance + sf_distance
## Model 2: spike_rate ~ ori_distance + sf_distance + attention_fac
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    157 3238.3                                  
## 2    156 2195.6  1    1042.7 74.084 7.617e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Part 8 - Test an interaction

The additive model above assumes that the relationship between orientation distance and spike rate is the same in both attention conditions. To test whether that slope differs by attention condition, add an interaction.

fit_interact <- lm(spike_rate ~ sf_distance + ori_distance * attention_fac,
                   data = d)
summary(fit_interact)
## 
## Call:
## lm(formula = spike_rate ~ sf_distance + ori_distance * attention_fac, 
##     data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.9914  -2.4585  -0.2718   2.7918   8.4052 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      49.19252    1.00778  48.813  < 2e-16 ***
## sf_distance                      -6.12473    0.38001 -16.117  < 2e-16 ***
## ori_distance                     -0.21684    0.01545 -14.038  < 2e-16 ***
## attention_facattend               4.24124    1.19828   3.539  0.00053 ***
## ori_distance:attention_facattend  0.02377    0.02258   1.053  0.29408    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.75 on 155 degrees of freedom
## Multiple R-squared:  0.8166, Adjusted R-squared:  0.8119 
## F-statistic: 172.6 on 4 and 155 DF,  p-value: < 2.2e-16

Interpret the interaction term as the difference between attention conditions in the ori_distance slope, after accounting for sf_distance.

ggplot(d, aes(x = ori_distance, y = spike_rate, colour = attention_fac)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    x = "Distance from preferred orientation (degrees)",
    y = "Spike rate",
    colour = "Attention",
    title = "Attention-specific regression lines"
  ) +
  scale_colour_manual(values = COL)
Regression lines from the interaction model.

Regression lines from the interaction model.

Apply this to your assigned dataset

Your assigned dataset is determined by your student ID number. Take the last digit and compute last_digit %% 3 in R.

Load your assigned dataset and fit a multiple regression model. Specifically:

  1. Decide on an independent unit of analysis before you fit the model.
  2. Create one summary row per participant or unit if needed.
  3. Identify a continuous outcome and at least two continuous predictors.
  4. Fit a simple regression, then add a second predictor.
  5. Add a categorical predictor.
  6. Compare models with anova().
  7. Check the residual plots.
  8. Test an interaction.
  9. Write a short paragraph interpreting the final model.