lapply: 拟合数千个混合模型并能够提取lsmeans

Question

我有一个包含超过10,000个线性混合模型(lme4)公式的列表，已经成功将它们拟合到一组数据中。我使用了lapply()函数以及一个包含tryCatch()函数的自定义函数来拟合这些模型。现在，我想提取所有这些模型的P值和lsmeans值。我成功地提取了P值，但是lsmeans函数遇到了错误。

library(lme4)
library(lmerTest)
library(pbkrtest)
library(lsmeans)

formulaS <- list() #Not going to detail generation of list, generically: 'Yvar~X1*X2+(1|subject)'
dataSET <- X #dataframe with first 3 columns containing fixed and random factors, 
             # as well as >10,000 columns of variables of interest

modelSeq <- function (x, dat) {
  return(tryCatch(lmer(x, data = dat), error=function(e) NULL))
}

modelsOutput <- lapply(formulaS, function(x) modelSeq(x, dat = dataSET))

lsmeans(modelsOutput[[1]], pairwise ~ X1:X2) #recieves error

解决默认值中的错误(L %% V0 %% t(L), L) : Lapack例程dgesv: 系统恰好是奇异的：U[1,1] = 0

我认为这不是模型问题的原因是，如果我单独拟合这些模型，我可以很好地提取lsmeans。是否有任何关于1)为什么我不能提取lsmeans、2)如何高效地提取均值、或3)替代的高效方法的评论。

谢谢！

__ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __

更新与编辑：这是我正在处理的RNAseq数据，具有随时间重复的受试者的多个样本，因此>10,000个模型具有描述实验设计的相同固定和随机效应。响应（基因）是唯一变化的变量。我已经在下面的代码中尝试了更明确地表达这一点。由于混合模型具有身份链接可能不适合这些数据，因此我有以下新包装器。我仍然无法提取平均值。此外，任何关于计算P值的更合适且时间有效的方法的评论都将不胜感激。

library(lme4)
library(blmeco)
library(ggeffects)

formulaS <- list() #Not going to detail generation of list, generically: 'GeneI~TRT*TIME+(1|subject)'
dataSET <- X #dataframe with first 3 columns containing fixed and random factors, 
             # as well as >10,000 columns of variables of interest (gene TPM)

wrap.glmer.nb <- function (modelForm, dat) {
  m <- tryCatch(glmer.nb(formula = modelForm, data = dat), error = function(e) NULL)
  if (!is.null(m)) {
    m.disp <- tryCatch(dispersion_glmer(m), error = function(e) NULL)
    m.wald <- tryCatch(anova(m), error = function(e) NULL)
    m.means.c <- tryCatch(ggemmeans(model = m, terms = c('TRT')), error = function(e) NULL)
    m.means.e <- tryCatch(ggemmeans(model = m, terms = c('TIME')), error = function(e) NULL)
    m.means.cxe <- tryCatch(ggemmeans(model = m, terms = c('TRT', 'TIME')), error = function(e) NULL)
    x <- list(m.disp, m.wald, m.means.c, m.means.e, m.means.cxe)
    print(paste0('Done with a model at ', Sys.time()))
    return(x)
  } else{
    x <- m
    return(x)
  }
}

startTime <- Sys.time()
modelOUTPUTS <- lapply(formulaS, function(modelForm) wrap.glmer.nb(modelForm, dat = dataSET))
endTime <- Sys.time()
print(paste('Victory! The analysis took:', endTime - startTime))

- rsp

---

没有看到模型和数据，很难说清楚，但从我的经验来看，这种错误通常与模型拟合有关。是否可能一个二元变量或交互项完全是单侧的？在数据中只有一个值？ - svenhalvorson

从统计学家的角度来看，我无法想象有任何好的理由去拟合超过10k个模型到同一数据集中，特别是如果你关心p值的话（这些p值将无效）。很可能，你应该用完全不同的方法来实现你的实际目标。 - Roland

@Roland 每个人一个模型，共计10k个人？ - Conor Neilson

这可能是一个作用域问题。Lsmeans试图使用模型调用来重现数据，但由于公式有自己的环境，它可能正在错误的位置查找。我建议在lsmeans()调用中明确指定数据。 - Russ Lenth

如果 @rsp 解决了自己的问题，请发布答案；如果您认为您的问题不再有用，请删除它... - Ben Bolker

Answer 1

如果您在modelSeq()中添加一行代码，您原来的设置将会起作用：

modelSeq <- function (x, dat) {
  environment(x) <- environment()
  return(tryCatch(lmer(x, data = dat), error=function(e) NULL))
}

这将公式的环境设置为函数体的环境，使得可以找到名为dat的数据集。

类似的例子：

fitm <- function(formula, data, ...) {
    environment(formula) <- environment()
    lm(formula, data = data, ...)
}

fl <- list(breaks ~ tension, breaks ~ wool + tension, breaks ~ wool*tension)

md <- lapply(fl, fitm, data = warpbreaks[c(1,2,3,5,8,13,21,34,54), ])

lapply(md, function(m) emmeans(m, "tension"))

这将产生：

NOTE: Results may be misleading due to involvement in interactions

[[1]]
 tension emmean    SE df lower.CL upper.CL
 L         41.2  6.64  6    24.91     57.4
 M         17.0 16.27  6   -22.82     56.8
 H         26.0 11.51  6    -2.16     54.2

Confidence level used: 0.95 

[[2]]
 tension emmean    SE df lower.CL upper.CL
 L         41.6  8.91  5    18.73     64.5
 M         17.7 19.41  5   -32.21     67.6
 H         26.0 12.59  5    -6.38     58.4

Results are averaged over the levels of: wool 
Confidence level used: 0.95 

[[3]]
 tension emmean   SE df lower.CL upper.CL
 L         41.1 10.9  4     10.9     71.3
 M       nonEst   NA NA       NA       NA
 H         26.0 14.1  4    -13.0     65.0

Results are averaged over the levels of: wool 
Confidence level used: 0.95 

顺便提一下，你不需要安装lsmeans包；它只是emmeans的前端。实际上，lsmeans函数本身就在emmeans中；它只是运行emmeans并重新标记结果。

- Russ Lenth

---

嗨，我偶然发现了这个问题。正如您在我的帖子中所看到的，我也遇到了相同的问题，因为我还需要提取与比较相关的统计数据。我不太确定该怎么做。所以如果您不介意的话，您的帮助将不胜感激。 - 12666727b9