基于lmer提取fMRI研究中随机效应的方差矩阵

Hi there! I’ve worked with mixed models for neuroimaging data before, so I can help you sort out how to extract and reshape those random effect variance components properly. Let’s break this down step by step:

First, Fix Your Dataset (Critical Missing Info!) #

The biggest issue with your current setup is that you don’t track which element of your stacked vector corresponds to which (row, column) position in the original 20x20 node matrices. Without this index info, you can’t map variance values back to the brain node structure. Let’s fix that first:

library(lme4)
set.seed(49830)

# Helper function to convert matrix to data frame with node indices
mat_to_df <- function(mat, subj_id) {
  # Get all (node_i, node_j) pairs
  node_pairs <- expand.grid(node_i = 1:nrow(mat), node_j = 1:ncol(mat))
  # Stack matrix values and combine with indices
  data.frame(
    subj = factor(subj_id),
    node_i = node_pairs$node_i,
    node_j = node_pairs$node_j,
    value = as.vector(mat)
  )
}

# Create your subject matrices and convert to indexed data frames
S.1 <- matrix(rexp(400, rate = 10), 20, 20)
S.2 <- matrix(rexp(400, rate = 10), 20, 20)
S.3 <- matrix(rexp(400, rate = 10), 20, 20)

df_list <- lapply(1:3, function(x) mat_to_df(get(paste0("S.", x)), x))
data.comp <- do.call(rbind, df_list)

# Optional: Center your predictor if needed (matches your original step)
data.comp$constant <- scale(data.comp$value, center = TRUE, scale = FALSE)[, 1]

Adjust Your Model to Match Node-Level Variability #

Your original model (constant|Group.all) is fitting a single random intercept and slope per subject, which doesn’t capture node-specific variability (the whole point of your fMRI analysis!). Instead, we want to model node-specific random effects across subjects—so each (node_i, node_j) pair has its own random deviation per subject.

Here’s the appropriate model:

# First, create a unique identifier for each node pair
data.comp$node_pair <- interaction(data.comp$node_i, data.comp$node_j, sep = "_")

# Fit model with fixed effects for each node pair, and random effects per subject per node pair
mod <- lmer(
  value ~ 0 + node_pair + (0 + node_pair | subj),
  data = data.comp,
  REML = TRUE,
  control = lmerControl(calc.derivs = FALSE)
)

• 0 + node_pair removes the global intercept and fits a fixed mean for each node pair.
• (0 + node_pair | subj) fits a random deviation for each node pair, per subject—this is where your node-specific variance lives.

Extract Random Effect Variance Components #

Now we can pull out the variance for each node pair. The easiest way is using VarCorr():

# Extract variance-covariance of random effects
vc <- VarCorr(mod)
subj_var_df <- as.data.frame(vc$subj)

# The result has row names matching your node pairs (e.g., "1_1", "1_2") and a `vcov` column with variance
head(subj_var_df)

Convert Variance Back to 20x20 Matrix #

Finally, map these variance values back to the original node matrix structure:

# Initialize empty 20x20 variance matrix
var_matrix <- matrix(NA, nrow = 20, ncol = 20)

# Fill the matrix using node indices from the variance data frame
for (pair in rownames(subj_var_df)) {
  # Split the pair identifier into row and column indices
  idx <- as.integer(strsplit(pair, "_")[[1]])
  i <- idx[1]
  j <- idx[2]
  # Assign variance to the correct matrix position
  var_matrix[i, j] <- subj_var_df[pair, "vcov"]
}

# Now var_matrix is your 20x20 variance matrix matching the original node structure!

Why Your Original Pseudocode Didn’t Work #

The code you shared calculates the variance-covariance matrix of the entire response vector y, which is a huge 1200x1200 matrix (3 subjects × 400 nodes). That’s not the node-specific random effect variance you need—we’ve focused on isolating the variance of the subject-specific deviations for each node pair instead.

内容的提问来源于stack exchange，提问作者Hevioss