Lab 3 Reference Guide

Overview

Subjects listen to podcast stories while an fMRI scanner records brain activity every 2 seconds at ~94,000 voxels. We build encoding models: text embeddings in, predicted BOLD signal out.

Neuroimaging Terms

fMRI and the BOLD Signal

Hemodynamic Response

Consequence for modeling: if a subject hears the word "dog" at time $t$, the brain's BOLD response to "dog" appears at roughly $t + 4$ to $t + 8$ seconds. This is why we need the make_delayed function: it gives the ridge regression access to features from 2, 4, 6, and 8 seconds before each BOLD measurement.

Text Embeddings

An embedding converts a word (or subword) into a numerical vector. This is what allows us to use text as input to a regression model.

Static Embeddings

BERT

Key properties:

• Bidirectional: attends to both preceding and following words (unlike GPT, which is left-to-right only)
• Max sequence length: 512 tokens. Long stories must be processed in overlapping windows
• Uses WordPiece tokenizer, which splits rare words into subword units (no OOV problems)
• Embeddings are typically extracted from the last hidden layer, with mean pooling over subwords

Fine-tuning vs transfer learning:

• Transfer learning: freeze BERT's weights, only train new layers on top (e.g., a classification head). Embeddings come from the frozen model
• Fine-tuning: allow BERT's weights to update during a second round of training on your specific data. More expressive but risks overfitting on small datasets

LoRA (Low-Rank Adaptation)

Low-rank intuition: a full $768 \times 768$ matrix has ~590K parameters per layer. A rank-4 factorization has only $768 \times 4 + 768 \times 4 = 6{,}144$ parameters. The assumption is that the direction BERT needs to shift for this task lives in a low-dimensional subspace, like PCA for weight updates.

Why LoRA for this lab? Full fine-tuning of BERT on ~24k words of story transcripts risks overfitting. LoRA's low parameter count acts as implicit regularization while still allowing the model to adapt its contextual representations to the narrative domain. Students fine-tune with LoRA on the MLM objective (same pretraining task, different data), then extract embeddings from the adapted model.

The Modeling Pipeline

Step 1: Lanczos Downsampling

Words occur at irregular intervals (~3 per second on average, but varying). The fMRI scanner captures snapshots at regular intervals (every TR $\approx$ 2 seconds). We need to convert from "one embedding per word, at irregular times" to "one embedding per timepoint, at regular times," which means interpolating onto a new time grid.

Example: adollshouse has $W_s = 1656$ words and tr_times has 256 entries, so downsampling gives $1656 \times 768 \to 256 \times 768$. After dropping the 5 + 10 padding rows (Step 3), you get $241 \times 768$, matching Y. Each row is a weighted blend of nearby word embeddings, with words closest in time weighted most heavily.

Step 2: Temporal Delays

Due to the HRF, the brain's response to a word heard at time $t$ appears in the fMRI data at approximately $t + 4$ to $t + 8$ seconds. The function make_delayed(stim, delays=[1,2,3,4]) accounts for this.

Why these specific delays? The HRF peaks at 4-6 seconds post-stimulus. With TR $\approx$ 2s, delays of 1-4 TR intervals correspond to 2, 4, 6, and 8 seconds. The ridge regression learns a weight for each delay at each voxel, effectively fitting the HRF shape per voxel. Some voxels respond faster (more weight on delay 1-2), others slower (more weight on delay 3-4).

Step 3: Aligning X to Y

The scanner started ~10 seconds before each story (giving 5 TRs with negative timestamps) and continued after the story ended. The raw data had 256 timepoints for adollshouse, but the Y data provided to you has only 241 rows: 15 timepoints were removed during preprocessing because fMRI responses at the edges of stories are noisier due to scanner transients and detrending artifacts.

When you downsample your embeddings using all 256 tr_times entries, you get 256 rows. You need to drop the same 15 rows that were removed from Y so the dimensions match. The question is which 15.

The Huth lab's own tutorial code uses [5+trim:-trim] with trim=5, giving [10:-5]: drop 10 from the start (5 pre-story padding + 5 noise buffer) and 5 from the end (noise buffer). The lab instructions say "trim the first 5 and last 10," which reverses this, but the [10:-5] split has stronger reference support and tends to give better correlations in practice. Either way, you need to remove exactly 15 rows total, and document your choice.

If you want to trim a few additional TRs from both X and Y beyond these 15 (e.g., because the BOLD signal still shows some edge effects in the kept region), that's a defensible judgment call. The edge artifacts are gradual rather than sharp, so the pre-trimming removed the worst of it but didn't eliminate residual effects entirely.

Stacking stories: after alignment, the per-story matrices are vertically concatenated across all training stories to form the full design matrix $X \in \mathbb{R}^{T \times 4p}$ and response matrix $Y \in \mathbb{R}^{T \times V}$, where $T = \sum_s T_s \approx 23{,}000$ is the total number of timepoints across all stories.

Step 4: Ridge Regression

The approach is sometimes called "mass-univariate" in the neuroimaging literature: we fit one independent regression model per voxel (~94k models), all sharing the same feature matrix $X$ but with different target columns $y_i$. Here's the notation.

Cross-validation: leave-one-story-out (LOSO). Train on all stories except one, predict on the held-out story. The regularization parameter $\alpha_{\text{ridge}}$ is selected via grid search over the LOSO folds.

Train/Test Split Strategies

The lab instructions say "split the stories into a training and test set" but leave the details to you. There are several reasonable approaches, and the choice matters for how you report and compare results.

Fixed holdout: reserve a small number of stories (e.g., 2-5) as a test set, train on the rest. Simple, but your CC estimates depend on which specific stories you held out. A long, complex story might be harder to predict than a short simple one, so the particular test set can shift your numbers noticeably.

Leave-one-story-out (LOSO): hold out each story in turn, train on the remaining ~100, predict the held-out story. This gives you a CC distribution across held-out stories rather than a single number, which is more informative. It's also what Huth et al. use for $\alpha_{\text{ridge}}$ selection. The downside is computational cost (101 ridge fits per embedding method), though the incremental trick mentioned above helps a lot.

Multiple random splits: randomly partition stories into train/test several times, fit and evaluate each split, report averaged CCs. This gives you error bars on your evaluation metrics, which strengthens any comparison between embedding methods. If BERT beats GloVe on one split but not another, that's important to know.

Whatever you choose, use the same splits across all embedding methods so comparisons are fair. And if you report a single number (e.g., "mean CC = 0.12"), be clear about what it's averaging over: voxels, stories, splits, or some combination.

Fit entirely separate models per subject. Subjects have different numbers of voxels and different brain geometry, so sharing data across subjects would introduce data leakage and wouldn't make sense anatomically.

Interpretation (Lab 3.3)

The ridge weights $\hat{\beta}_i \in \mathbb{R}^{4p}$ don't tell you "which words matter" directly, because each BERT embedding dimension doesn't correspond to anything human-readable. (With one-hot bag-of-words you'd get this for free, but we lose it with distributed representations.) To get word-level interpretations, we use post-hoc attribution methods.

SHAP

Theoretically grounded (additive, symmetric), but computationally expensive. In practice, tends to highlight function words and discourse markers.

LIME

Fast and model-agnostic, but deletion changes the sequence length and shifts positions, which can be problematic for transformers. Tends to highlight content words rather than structural ones.

The Wrapper Pattern (how to make SHAP/LIME work here)

SHAP and LIME expect a function f(text) → prediction, but our model is a ridge regression on delayed, downsampled embeddings. The trick is to wrap the entire pipeline in a class so SHAP/LIME can perturb the raw text and observe how predictions change.

class fMRI_Predictor:
    def __init__(self, ridge_model, tokenizer, bert_model,
                 data_times, tr_times, voxel_index):
        # Store all pipeline components
        ...

    def predict(self, texts):
        # texts: list of strings (perturbed versions from SHAP/LIME)
        predictions = []
        for text in texts:
            words = text.split()
            embeddings = get_bert_embeddings(words)  # re-embed
            downsampled = lanczos_downsample(embeddings)
            delayed = make_delayed(downsampled, [1,2,3,4])
            pred = self.ridge_model.predict(delayed)
            predictions.append(pred[:, self.voxel_index].mean())
        return np.array(predictions)

Then hand the wrapper to the explainers:

# SHAP: masks words with [MASK], preserving position
masker = shap.maskers.Text(r"\s+", mask_token="[MASK]")
explainer = shap.Explainer(predictor.predict, masker)
shap_values = explainer([text_input])

# LIME: deletes words, fits local linear surrogate
lime_explainer = LimeTextExplainer()
exp = lime_explainer.explain_instance(
    text_input, predictor.predict, num_features=10)

SHAP replaces words with [MASK], so BERT still sees a token at that position. LIME deletes words entirely, changing the sequence length. This is why they tend to flag different words: SHAP picks up positional/structural importance, LIME picks up semantic content.

Practical note: this is slow because every perturbation requires a full forward pass through BERT + the preprocessing pipeline. For LIME with num_samples=100, that's 100 BERT forward passes per voxel. Run this on GPU via sbatch, not the login node.

PCS Framework

Data Summary

Transcripts: raw_text.pkl contains 109 stories (superset of fMRI stories). Each story is a DataSequence object with: word tokens (data), word onset times (data_times), and TR times (tr_times). Loading requires ridge_utils on the Python path.

Memory math: full Subject 2 in float64 $\approx$ 18 GB; in float32 $\approx$ 9 GB; single story in float32 $\approx$ 87 MB.

Infrastructure

Bridges-2 (PSC)

Pittsburgh Supercomputing Center HPC cluster. Data at: /ocean/projects/mth250011p/shared/215a/final_project/data. Use /ocean for storage (not your home directory, which has limited quota). Shared class storage: 512 GB (provisioning an additional 512 GB).

SCF (Berkeley Stats)

Stats Computing Facility, local alternative. Data mirrored at: /scratch/users/s214/lab3. Use /scratch/users/your_username for your working files. Some students have reported permission issues accessing the shared data path.

HuggingFace

100GB free storage, useful for sharing large files (embeddings, trained models) within your group. Can be accessed programmatically from either cluster.

Compute Rules

Do not run jobs on login nodes. This slows down the entire cluster for everyone. Submit jobs via sbatch or use an interactive JupyterHub session. Request memory explicitly (e.g., --mem=64G).

Only use GPUs for BERT fine-tuning and embedding extraction. Ridge regression and scikit-learn operations run on CPU only and won't benefit from a GPU. Don't hold a GPU allocation while running ridge fits.

Common Pain Points

• Kernel crashes: usually from running on the login node or insufficient memory allocation
• Memory: cast to float32 (halves memory), use memory-mapped loading (np.load(..., mmap_mode='r')), process voxels in batches. Bag-of-words embeddings are sparse (scipy.sparse), but Word2Vec, GloVe, and BERT embeddings are dense
• Storage: don't store unnecessarily large files. Compress where possible
• Permissions: ensure GitHub org and classroom invitations are accepted

Timeline & Deadlines

Quick Reference