1

R でシミュレーションベースの電力分析を行っています。R を RStudio (0.98.932) で実行し、関数plyr::rdplyとを使用してlme4::glmer、それぞれデータを生成し、モデルを適合させます (R 環境とパッケージ バージョンについては、以下の再現可能な例の最後を参照してください)。 ）。

手順は、特定のパラメーター化のデータセットをランダムに生成し、それにモデルを適合させることです。ただし、時々、モデルは収束に失敗します。これが発生すると、次の警告が表示されます

[1]「スケーリングされた勾配を評価できません」
[2]「モデルは収束に失敗しました: 1 つの負の固有値を持つ縮退ヘッセ行列」

R はブラウザ モードに入りc、シミュレーション ループに戻るには手動で介入する必要があります (たとえば、 を押します)。数日間にわたって何千もの反復を実行する必要があるため、これは本当に苦痛ですが、この特定の収束エラーが発生するたびに、キーを押すまで停止します。

Rがブラウザモードに入らないようにする方法はありますか? 各シミュレーションで発生するすべての警告を保存するので、唯一の問題は、この特定の収束エラーが発生したときに手動で介入しなければならないことです。purrr::quietly関数と関数を使用してみましpurrr::safelyたが、成功しませんでした (以下のコードの例を参照)。

これは私のコンピューターで動作するMWEです（set.seed再現性のために使用しているので、パッケージのバージョンなどに関係なく同じ結果になることを願っています）。この例では、同じロジックを適用しますが、実際のシミュレーションで使用するのとは異なり、より単純なパラメーター化を適用します。

library(lme4)
library(plyr)
library(purrr)

# function to generate data that will lead to convergence failure
mini_simulator <- function() {
  nb_items <- 10  # observations per subject
  nb_subj <- 10  # subjects per group
  generate_data <- function() {
    A <- rbinom(nb_items * nb_subj, 1, .99)
    B <- rbinom(nb_items * nb_subj, 1, .8)
    simdata <- data.frame(
      Group = rep(c("A", "B"), each = nb_items * nb_subj),
      Subj = rep(1 : (nb_subj * 2), each = nb_items),
      Items = 1:nb_items,
      Response = c(A, B)
      )
  }
}

# Sanity check that the function is generating data appropriately.
# d should be a dataframe with 200 obs. of 4 variables
d <- mini_simulator()()
head(d, 3)
# Group Subj Items Response
# 1     A    1     1        1
# 2     A    1     2        1
# 3     A    1     3        1
rm(d)

## Functions to fit model

# basic function to fit model on simulated data
fit_model <- function(data_sim) {
  fm <- glmer(
    formula = Response ~ Group + (1|Subj) + (1|Items),
    data = data_sim, family = "binomial")
  out <- data.frame(summary(fm)$coef)
  out
}

# similar but using purrr::quietly (also tried purrr::safely with no success)
# see http://r4ds.had.co.nz/lists.html section "Dealing with failure"
fit_model_quietly <- function(data_sim) {
  purrr_out <- purrr::quietly(glmer)(
    formula = Response ~ Group + (1|Subj) + (1|Items),
    data = data_sim, family = "binomial")
  fm <- purrr_out$result
  out <- data.frame(summary(fm)$coef)
  # keeps track of convergence failures and other warnings
  out$Warnings <- paste(unlist(purrr_out$warnings), collapse = "; ")
  out
}

# this seed creates the problematic convergence failure on the first evaluation
# of rdply
set.seed(2)
# When I run the next line R goes into Browse mode and I need to enter "c"
# in order to continue
simulations <- plyr::rdply(.n = 3, fit_model(mini_simulator()()))
simulations

# problem persists using the quietly adverb from purrr
set.seed(2)
simulations <- plyr::rdply(.n = 3, fit_model_quietly(mini_simulator()()))
simulations

# sessionInfo()

# R version 3.1.2 (2014-10-31)
# Platform: i386-w64-mingw32/i386 (32-bit)
# 
# locale:
#   [1] LC_COLLATE=Swedish_Sweden.1252  LC_CTYPE=Swedish_Sweden.1252    LC_MONETARY=Swedish_Sweden.1252
# [4] LC_NUMERIC=C                    LC_TIME=Swedish_Sweden.1252    
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] purrr_0.2.1  plyr_1.8.1   lme4_1.1-8   Matrix_1.1-4
# 
# loaded via a namespace (and not attached):
#   [1] grid_3.1.2      lattice_0.20-29 magrittr_1.5    MASS_7.3-35     minqa_1.2.4     nlme_3.1-118    nloptr_1.0.4   
# [8] Rcpp_0.11.3     splines_3.1.2   tools_3.1.2

更新 (r2evans のコメントに基づく)

両方のコンピューターoptions("error")で利回り

(関数 () { .rs.breakOnError(TRUE) })()

これはある種の RStudio のデフォルトのようであり、実際にstop()呼び出しが発生したときに R をブラウザー モードにするようです (これは、メニュー ツールバーの [デバッグ] > [エラー時] > ... を介してグラフィカル インターフェイスで変更できることがわかります)。とにかく、設定するoptions(error = NULL)と、問題はなくなります。これは、うまく機能する新しい (簡略化された) 例です (この最小限の例でも、実際のシミュレーションに適用された場合でも)。

library(lme4)
library(plyr)
library(purrr)

options(error=NULL)

## Function to generate data
# Generates data that will lead to convergence failure
mini_simulator <- function() {
  nb_items <- 10  # observations per subject
  nb_subj <- 10  # subjects per group
  generate_data <- function() {
    A <- rbinom(nb_items * nb_subj, 1, .99)
    B <- rbinom(nb_items * nb_subj, 1, .8)
    simdata <- data.frame(
      Group = rep(c("A", "B"), each = nb_items * nb_subj),
      Subj = rep(1 : (nb_subj * 2), each = nb_items),
      Items = 1:nb_items,
      Response = c(A, B)
    )
  }
}

## Function to fit model
# Fits model on simulated data with purrr::quietly to capture warnings
# (http://r4ds.had.co.nz/lists.html section "Dealing with failure")
fit_model_quietly <- function(data_sim) {
  purrr_out <- purrr::quietly(glmer)(
    formula = Response ~ Group + (1|Subj) + (1|Items),
    data = data_sim, family = "binomial")
  fm <- purrr_out$result
  out <- data.frame(summary(fm)$coef)
  # keeps track of convergence failures and other warnings
  out$Warnings <- paste(unlist(purrr_out$warnings), collapse = "; ")
  out
}

# this seed creates the problematic convergence failure on the first evaluation
# of rdply
set.seed(2)
simulations <- plyr::rdply(.n = 3, fit_model_quietly(mini_simulator()()))
simulations
4

1 に答える 1