GLM代码运行卡顿崩溃求助：逻辑回归执行过久致电脑冻结

我在运行以下R代码时遇到了严重的性能问题：电脑持续冻结，只能查看相关矩阵，无法获取GLM结果和Datarray数据。已经尝试缩小n_people的范围，但问题依旧，电脑甚至出现异响，最后只能强制退出程序。恳请大家提供解决建议：

# running the assay
#which_p_value = "x1"
which_p_value = "groupcategory"
#which_p_value = "x1:groupcategory"
run_anova = FALSE
simulate_mixed_effect = TRUE
mixed_effect_sd = 3.094069
mixed_effect_sd_slope = 3.098661
library(tidyverse)
n_people <- c(2,5,10,15,20)
coef1 <- 1.61
coef2 <- -0.01
#coef3 <- 5
#coef4 <- 0
g1 = 0
g2 = 1
g3 = 2
distances <- c(60,90,135,202.5,303.75,455.625)/100
n_trials <- 35
oneto1000 <- 25
n_track_lengths <- length(distances)
groupcategory = c(rep(g1, n_track_lengths), rep(g2, n_track_lengths),rep(g3,n_track_lengths))
z = c(n_people)
emptydataframeforpowerplots = NULL
coef3s <- c(-5, -4, -3, -2,-1, 0, 1, 2, 3, 4, 5)
coef4s <- c(-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1)
Datarray <- array(dim=c(length(coef3s), length(coef4s),length(n_people)))
coef3_counter =1
for (coef3 in coef3s) {
 coef4_counter =1
 for (coef4 in coef4s) {
 z1_g2 <- coef1 + coef2*distances + coef3*g2 + coef4*g2*distances
 z1_g3 <- coef1 + coef2*distances + coef3*g3 + coef4*g3*distances
 d = NULL
 pr1 = 1/(1+exp(-z1_g2))
 pr2 = 1/(1+exp(-z1_g3))
 counter=1
 for (i in n_people) {
 for (j in 1:oneto1000){
 df <- c()
 for (k in 1:i){
 # random effect from drawing a random intercept with sd = x
 if (simulate_mixed_effect){
 coef1_r = rnorm(1, mean=coef1, sd=mixed_effect_sd)
 coef2_r = rnorm(1, mean=coef1, sd=mixed_effect_sd_slope)
 } else {
 coef1_r = coef1
 coef2_r = coef2
 }
 z_g1 <- coef1_r + coef2*distances + coef3*g1 + coef4*g1*distances
 pr = 1/(1+exp(-z_g1))
 z1_g2 <- coef1_r + coef2*distances + coef3*g2 + coef4*g2*distances
 pr1 = 1/(1+exp(-z1_g2))
 if (run_anova) {
 df <- rbind(df, data.frame(x1 = c(rep(distances, 3)), y = c(rbinom(n_track_lengths,n_trials,pr), rbinom(n_track_lengths,n_trials,pr1),rbinom(n_track_lengths,n_trials,pr2)), groupcategory = groupcategory, id = c(rep(k,18))))
 } else {
 # this is for glmer data organisation
 for (m in 1:n_trials) {
 df <- rbind(df, data.frame(x1 = c(rep(distances, 3)), y = c(rbinom(n_track_lengths,1,pr),rbinom(n_track_lengths,1,pr1),rbinom(n_track_lengths,1,pr2)),groupcategory = groupcategory,id = c(rep(k,18))))
 }
 }
 }
 if (run_anova) {
 #df_aov <- aov(y~x1*groupcategory+Error(id/(x1*groupcategory)),data=df)
 #df_aov_sum <- summary(df_aov)
 #pvalue <- df_aov_sum[[5]][[1]][which_p_value,"Pr(>F)"]
 df_aov <- aov(y~x1*groupcategory+Error(id),data=df)
 df_aov_sum <- summary(df_aov)
 pvalue <- df_aov_sum[[2]][[1]][which_p_value, "Pr(>F)"]
 }
 checkme <- df %>% group_by(groupcategory,id) %>% summarise(miny=min(y),maxy=max(y)) %>% mutate(expectfail = miny==maxy)
 else {
 mod_group_glmer <- glmer(y ~ x1 + groupcategory + (1+x1|id), data = df, family = "binomial")
 sum <- summary(mod_group_glmer)
 pvalue <- sum$coefficients[which_p_value, "Pr(>|z|)"]
 }
 d = rbind(d,data.frame(pvalue))
 }
 count <- plyr::ldply(d,function(c) sum(c<=0.05))
 Datarray[coef3_counter,coef4_counter,counter] <- count$V1/oneto1000
 counter = counter +1
 d = NULL
 }
 coef4_counter = coef4_counter + 1
 }
 coef3_counter = coef3_counter + 1
 }

哇，这种代码跑着跑着直接卡死还异响的情况真的太闹心了！我帮你拆解下问题根源，再给你几个可落地的优化方向：

1. 先修复语法错误（最紧急） #

你代码里有个明显的语法bug：在if (run_anova)分支结束后，直接写了checkme <- ...然后接else，这会导致逻辑完全混乱，R根本不知道这个else对应哪个if！这不仅会让代码报错，还可能导致无意义的计算，加剧资源消耗。

把这段代码改成这样（给checkme加个判断逻辑，或者如果是冗余代码直接删掉）：

if (run_anova) {
  # 原anova代码...
} else {
  # 先做checkme的判断（如果需要的话）
  checkme <- df %>% 
    group_by(groupcategory,id) %>% 
    summarise(miny=min(y),maxy=max(y)) %>% 
    mutate(expectfail = miny==maxy)
  # 再拟合glmer
  mod_group_glmer <- glmer(y ~ x1 + groupcategory + (1+x1|id), data = df, family = "binomial")
  sum <- summary(mod_group_glmer)
  pvalue <- sum$coefficients[which_p_value, "Pr(>|z|)"]
}

2. 干掉嵌套循环里的rbind（内存杀手） #

你现在每次循环都用rbind(df, ...)拼接数据框，R的rbind每次都会复制整个数据框，循环次数多了会导致内存碎片化严重，越跑越慢。预分配数据框或者用列表存储再合并才是正确姿势：

比如在k循环前，先创建一个空列表：

df_list <- list()
for (k in 1:i){
  # 生成单个人的data.frame，存到列表里
  person_df <- data.frame(...) # 这里写原来生成单条数据的代码
  df_list[[k]] <- person_df
}
# 最后一次性合并
df <- bind_rows(df_list)

同样，d这个存储pvalue的数据框，也可以改成列表存储，最后用bind_rows合并，不要每次rbind。

3. 优化glmer的拟合效率 #

混合效应模型拟合本身就是计算密集型任务，你还要循环多次，这对CPU和内存都是极大的考验：

• 先缩小测试规模：把oneto1000先改成10或者20，确认代码逻辑没问题、能正常输出结果后，再逐步增大数值。
• 换更高效的优化器：给glmer加控制参数，比如control = glmerControl(optimizer = "bobyqa")，这个优化器比默认的Nelder-Mead更适合二项式模型，收敛更快。
• 试试glmmTMB包：glmmTMB是lme4的替代包，处理二项式混合模型时速度更快，而且对边界情况的处理更稳定，语法和glmer几乎一致，直接替换就行：

  library(glmmTMB)
  mod_group_glmer <- glmmTMB(y ~ x1 + groupcategory + (1+x1|id), data = df, family = "binomial")
  

4. 减少重复计算 #

代码里很多重复计算的部分可以提前预处理，比如distances、groupcategory这些固定向量，还有z1_g2、z1_g3在coef4循环里其实可以提前算好，不用每次都重复计算，减少循环内的计算量。

5. 手动清理内存 #

每次循环结束后，手动删除没用的对象并强制垃圾回收，释放内存：

# 在j循环结束后
rm(df, mod_group_glmer, sum, checkme)
gc() # 强制R回收内存

6. 用并行计算分担压力 #

如果你的电脑是多核CPU，可以把oneto1000的循环改成并行的，把任务分到多个核心上，减少总运行时间，也避免单核心满载导致的卡顿。比如用foreach结合doParallel：

library(doParallel)
# 留一个核心给系统，避免卡死
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

# 把原来的j循环改成foreach
d <- foreach(j = 1:oneto1000, .combine = rbind, .packages = c("glmmTMB", "dplyr")) %dopar% {
  # 这里放原来j循环里的所有代码：生成df、拟合模型、计算pvalue
  # 最后返回data.frame(pvalue = pvalue)
}

stopCluster(cl)

最后一步：分步测试 #

先把所有参数调到最小（比如coef3s只取一个值，coef4s只取一个值，n_people=2，oneto1000=1），跑通整个流程，确认能正常输出Datarray，再逐步扩大参数范围，这样能快速定位是不是某个参数组合导致的问题。

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