lapply: подгонка тысяч смешанных моделей и возможность извлечения lsmeans

Question

У меня есть список формул (> 10000) для линейных смешанных моделей (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

Ошибка в solve.default (L% % V0% % t (L ), L): подпрограмма Лапака dgesv: система в точности единственная: U [1,1] = 0

Причина, по которой я не думаю, что это проблема модели, заключается в том, что если я подгоняю модели по отдельности Я могу извлечь lsmeans просто отлично. Есть ли какие-либо комментарии по поводу 1) почему я не могу извлечь lsmeans, 2) как эффективно извлечь средства или 3) альтернативный, эффективный метод.

Спасибо!

__ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ ______ данные с повторяющимися выборками предметов со временем, с которыми я играю, поэтому> 10000 моделей имеют те же фиксированные и случайные эффекты, которые описывают экспериментальный дизайн. Ответ (ген) является единственной переменной, которая варьируется. Я попытался сделать это более явным в приведенном ниже коде. Признавая, что смешанная модель с идентификационной связью не может быть идеальной для данных, у меня есть новая оболочка ниже. У меня все еще есть проблемы с извлечением средств. Также приветствуется любой комментарий по поводу более подходящих и эффективных по времени методов для вычисления значений 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))

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 и перемаркирует результаты.