如何快速为列拟合lm并执行ncvTest？批量lm拟合性能优化

Nice move switching to the matrix-based lm() syntax—this vectorized approach cuts down on all the redundant overhead of fitting thousands of separate linear models, which is why it’s so much faster! Let’s get those ncvTest() results for every column without losing that speed boost.

Background: How Matrix-Based lm() Works #

When you fit lm(as.matrix(y_cols) ~ x), the resulting model object stores residuals and fitted values as matrices (one column per original y variable). We can leverage these precomputed values to run ncvTest() (from the car package) efficiently, instead of refitting individual models.

Method 1: Batch Compute NCV Test Stats Directly #

The ncvTest() is essentially an F-test for heteroscedasticity, which tests if squared residuals are related to fitted values. We can compute this manually for each column using the precomputed residuals and fitted values:

library(car)

# Extract precomputed residuals and fitted values from your matrix lm model
resid_mat <- residuals(lmFitGA)
fitted_mat <- fitted(lmFitGA)
resid_sq_mat <- resid_mat ^ 2

# Define a helper function to calculate NCV stats for a single column
calc_ncv <- function(resid_sq_col, fitted_col) {
  # Fit the auxiliary regression for heteroscedasticity test
  aux_model <- lm(resid_sq_col ~ fitted_col)
  # Extract F-statistic and p-value from ANOVA
  anova_res <- anova(aux_model)
  list(
    F_statistic = anova_res$F[1],
    p_value = anova_res$`Pr(>F)`[1]
  )
}

# Apply the helper function to every column (fast, vectorized under the hood)
ncv_results <- mapply(
  calc_ncv,
  split(resid_sq_mat, col(resid_sq_mat)),
  split(fitted_mat, col(fitted_mat)),
  SIMPLIFY = FALSE
)

# Name the results to match your original columns
names(ncv_results) <- colnames(MyValues[,1:626374])

Method 2: Get Full ncvTest() Output #

If you want the exact same output format as calling ncvTest() on a single lm object, you can create lightweight "fake" lm objects for each column (using the precomputed residuals/fitted values) and pass them to ncvTest():

# Helper function to create a minimal lm object for ncvTest
make_fake_lm <- function(col_idx) {
  fake_lm <- list(
    residuals = resid_mat[, col_idx],
    fitted.values = fitted_mat[, col_idx],
    # Include the original data to match lm object structure
    model = data.frame(
      y = MyValues[, col_idx],
      x = MyValues$Gestational.Age
    )
  )
  class(fake_lm) <- "lm"
  fake_lm
}

# Run ncvTest on each fake lm object
full_ncv_results <- lapply(
  seq(ncol(resid_mat)),
  function(i) ncvTest(make_fake_lm(i))
)

# Name the results
names(full_ncv_results) <- colnames(MyValues[,1:626374])

Speed Up with Parallel Processing #

If you still want to use parallelization for the final step (though the matrix-based prep makes this less necessary), you can use foreach with doParallel:

library(doParallel)
registerDoParallel(cores = 10)

parallel_ncv_results <- foreach(
  i = seq(ncol(resid_mat)),
  .packages = "car",
  .export = c("lmFitGA", "MyValues")
) %dopar% {
  ncvTest(make_fake_lm(i))
}

names(parallel_ncv_results) <- colnames(MyValues[,1:626374])

Why This Works So Much Faster #

By using the matrix-based lm() first, we avoid fitting 600k+ separate linear models. All the heavy lifting (calculating coefficients, residuals, fitted values) is done in a single vectorized operation. The subsequent ncvTest() steps are just lightweight auxiliary regressions, which are far faster than refitting the original models.

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