lme4 allFit () дает непонятные результаты при переносе в функцию - PullRequest
Новая
       20

lme4 allFit () дает непонятные результаты при переносе в функцию

1 голос
/ 04 июля 2019

Я использую allFit() в lme4 для автоматического сканирования возможных оптимизаторов, поскольку оптимизатор по умолчанию обычно не сходится в этой ситуации.Мой код работает нормально, когда я запускаю it line-by-line, но когда я запускаю его, завернутый в простую функцию, он дает разные результаты.

Я посмотрел на вывод вызова allFit, и кажется, что когда онНЕ внутри функции, она возвращает список lmerModLmerTest объектов по желанию.

Однако внутри функции возвращается список со значениями simpleError, error и condition.Почему он это делает?

Я использую RStudio, R 3.6, lme4 1.1-21, lmerTest 3.1-0.

ОБНОВЛЕНИЕ: Проблема в том, что обновление ()Метод, используемый allFit, не может найти фрейм данных 'tt' при повторной подгонке моделей.Я установил в коде точки останова, и кажется, что в функциональной среде существуют «тестовые» данные, поэтому я не понимаю, почему они не могут их найти ...

ОБНОВЛЕНИЕ 2: Похоже, что если я изменю назначение тестовых данных на <<-, это сработает.Это опасно, однако, нарушая функциональное программирование, и я думаю, что оно может потерпеть неудачу, когда я пытаюсь распараллелить.Я тестирую дальше ... все еще открыт для предложений!

Вот код, который работает, а не внутри функции: 

library(lme4)

multi_arm_var_sim <- function(nsub = 20, nclust = 100, narm = 2, iccs = c(.01, .04), betas = c(0,.3)){
  sig_b2 <- -1*iccs / (iccs - 1)
  n <- nsub * nclust * narm
  y <- rep_len(NA, n)
  arm <- as.factor(rep(0:(narm-1), each = nsub*nclust))
  clustid <- rep(1:(nclust*narm), each = nsub)
  clustRElist <- rnorm(narm*nclust, mean = 0, sd = rep(sqrt(sig_b2), each = nclust))
  clustRE <- rep(clustRElist, each = nsub)
  sig_b2 <- rep(sig_b2, each = nclust*nsub)
  error <- rnorm(n, mean = 0, sd = 1)
  beta <- rep(betas, each = nclust*nsub)
  linpred <- beta + clustRE + error
  output <- cbind.data.frame(arm, clustid, sig_b2, clustRE, linpred)
  return(output)
}

set.seed(2)
test_1 <- multi_arm_var_sim()
model_flex_1 <- lmer(linpred ~ arm + (1 + arm | clustid),
                             data = test_1)

diff_optims_1 <- allFit(model_flex_1, verbose = TRUE)
print(class(diff_optims_1[[1]]))
is.OK_1 <- sapply(diff_optims_1, is, "lmerMod")
print(is.OK_1)

А вот код, который не 'т, та же настройка, завернутая в функцию. 

library(lme4)

multi_arm_var_sim <- function(nsub = 20, nclust = 100, narm = 2, iccs = c(.01, .04), betas = c(0,.3)){
  sig_b2 <- -1*iccs / (iccs - 1)
  n <- nsub * nclust * narm
  y <- rep_len(NA, n)
  arm <- as.factor(rep(0:(narm-1), each = nsub*nclust))
  clustid <- rep(1:(nclust*narm), each = nsub)
  clustRElist <- rnorm(narm*nclust, mean = 0, sd = rep(sqrt(sig_b2), each = nclust))
  clustRE <- rep(clustRElist, each = nsub)
  sig_b2 <- rep(sig_b2, each = nclust*nsub)
  error <- rnorm(n, mean = 0, sd = 1)
  beta <- rep(betas, each = nclust*nsub)
  linpred <- beta + clustRE + error
  output <- cbind.data.frame(arm, clustid, sig_b2, clustRE, linpred)
  return(output)
}

get_pval <- function(){
  tt <- multi_arm_var_sim()
  model_flex <- lme4::lmer(linpred ~ arm + (1 + arm | clustid),
                               data = tt)
  diff_optims <- lme4::allFit(model_flex, data = tt, verbose = TRUE)
  print(class(diff_optims[[1]]))
  is.OK <- sapply(diff_optims, is, "merMod")
  print(is.OK)
}
set.seed(2)
get_pval()

Спасибо !!

...