机器算法验证 - 为什么这个 lme4 置信区间的覆盖率小于 95%？ - 吾爱随笔录

为什么这个 lme4 置信区间的覆盖率小于 95%？

机器算法验证混合模式置信区间 lme4-nlme 模拟覆盖概率

2022-03-26 07:10:03

我有可以使用模型描述的数据 $y_{i j} = \alpha_i + \epsilon_{i j}$ ，在哪里 $\alpha_i \sim \text{N}(\mu_{\alpha}, \sigma_{\alpha}^2)$ 和 $\epsilon_{i j} \sim \text{N}(0, \sigma_{\epsilon}^2)$ . 有 48 组 768 个观测值，每组 16 个观测值，所以 $i = 1, \ldots, 48$ 和 $j = 1, \ldots, 16$ . 我想计算 95% 的置信区间 $\mu_{\alpha}$ . 我决定尝试拟合模型 $\tt{y \sim (1 \ | \ group)}$ 使用 lme4 包的数据并使用包的计算置信区间 $\tt{confint()}$ 功能。我创建了一个模拟来检查置信区间的覆盖范围；代码如下。覆盖率似乎在 94.5% 左右，而不是 95%。我的模拟中是否存在缺陷，或者这个结果是意料之中的？

# Load packages -----------------------------------------------------------

library(lme4)
library(tidyverse)

# Define functions --------------------------------------------------------

make_dataset <- function(groups, group_size,
                         mu_alpha, sigma_alpha, sigma_epsilon) {
  replicate(
    groups,
    {
      group_mean <- rnorm(1, mean = mu_alpha, sd = sigma_alpha)
      errors <- rnorm(group_size, sd = sigma_epsilon)
      tibble(y = group_mean + errors)
    },
    simplify = FALSE
  ) %>%
    bind_rows(.id = "group") %>%
    mutate(group = as_factor(group))
}

check_conf_int <- function(dataset, mu_alpha) {
  mod <- lmer(y ~ (1 | group), dataset)
  conf_int <- confint(mod, "(Intercept)")
  between(mu_alpha, conf_int[, "2.5 %"], conf_int[, "97.5 %"])
}

# Set simulation parameters -----------------------------------------------

groups <- 48
group_size <- 16
mu_alpha <- 0
sigma_alpha <- 1
sigma_epsilon <- 0.25
runs <- 10000

# Make datasets -----------------------------------------------------------

set.seed(1)

# Takes 5 minutes
datasets <- replicate(
  runs,
  make_dataset(groups, group_size, mu_alpha, sigma_alpha, sigma_epsilon),
  simplify = FALSE
)

# Run simulation ---------------------------------------------------------

# Takes 30 minutes
results <- map_lgl(datasets, check_conf_int, mu_alpha = mu_alpha)

tibble(
  runs = seq_len(runs),
  coverage = cummean(results)
) %>%
  ggplot(aes(runs, coverage)) +
  geom_line() +
  geom_hline(linetype = "dashed", yintercept = 0.95) +
  coord_cartesian(ylim = c(0.93, 0.97))

下面是最后一个代码块产生的图：

下面是输出 $\tt{sessionInfo()}$ ：

R version 4.0.0 (2020-04-24)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] forcats_0.5.0   stringr_1.4.0   dplyr_0.8.5     purrr_0.3.4     readr_1.3.1    
 [6] tidyr_1.0.3     tibble_3.0.1    ggplot2_3.3.0   tidyverse_1.3.0 lme4_1.1-23    
[11] Matrix_1.2-18  

loaded via a namespace (and not attached):
 [1] statmod_1.4.34   tidyselect_1.1.0 xfun_0.13        splines_4.0.0   
 [5] haven_2.2.0      lattice_0.20-41  colorspace_1.4-1 vctrs_0.3.0     
 [9] generics_0.0.2   htmltools_0.4.0  yaml_2.2.1       rlang_0.4.6     
[13] nloptr_1.2.2.1   pillar_1.4.4     withr_2.2.0      glue_1.4.0      
[17] DBI_1.1.0        dbplyr_1.4.3     modelr_0.1.7     readxl_1.3.1    
[21] lifecycle_0.2.0  cellranger_1.1.0 munsell_0.5.0    gtable_0.3.0    
[25] rvest_0.3.5      evaluate_0.14    labeling_0.3     knitr_1.28      
[29] fansi_0.4.1      broom_0.5.6      Rcpp_1.0.4.6     backports_1.1.6 
[33] scales_1.1.1     jsonlite_1.6.1   farver_2.0.3     fs_1.4.1        
[37] hms_0.5.3        packrat_0.5.0    digest_0.6.25    stringi_1.4.6   
[41] grid_4.0.0       cli_2.0.2        tools_4.0.0      magrittr_1.5    
[45] crayon_1.3.4     pkgconfig_2.0.3  MASS_7.3-51.6    ellipsis_0.3.0  
[49] xml2_1.3.2       reprex_0.3.0     lubridate_1.7.8  assertthat_0.2.1
[53] minqa_1.2.4      rmarkdown_2.1    httr_1.4.1       rstudioapi_0.11 
[57] R6_2.4.1         boot_1.3-25      nlme_3.1-147     compiler_4.0.0

更新

如果我将组增加到 480，下面是我得到的图：

这是我将 group_size 更改为 160 时得到的图：

1个回答

这里的问题是，随着组数变大，平均覆盖率只会收敛到 0.95。当组数“低”时，我相信您会发现平均覆盖率低于 95%。48既不是特别高也不是特别低。但是，由于您的参数的典型覆盖率为 94.5%，我认为在这种情况下这不是一个大问题。

此外，您在情节中看到的模式仅与您设置的特定种子相关。

为了看到这一点，我们可以运行更多的模拟。请注意，我已修改您的代码，使其在您可用的尽可能多的本地内核上并行运行。这种类型的模拟特别适合并行化。这样做将我机器上的运行时间从超过 28 分钟减少到不到 6 分钟，有 48 组，每组 16 名受试者:)

我还为主模拟添加了一个进度条，这样你就不用等着想知道需要多长时间了：

library(lme4)
library(tidyverse)
library(parallel)
library(purrr)
library(furrr)

check_conf_int <- function(dataset, mu_alpha) {
  mod <- lmer(y ~ (1 | group), dataset)
  conf_int <- confint(mod, "(Intercept)", quiet = TRUE)
  between(mu_alpha, conf_int[, "2.5 %"], conf_int[, "97.5 %"])
}

runs <- 10000
mu_alpha <- 0  # note : also needs to be set in the cluster

for (seed in 1:10) {

  set.seed(seed)

  cl <- makeCluster(detectCores()-1)  

  # send data, packages and function to the cores
  clusterEvalQ(cl,
             {
              library(magrittr);library(dplyr)
              groups <- 96
              group_size <- 16
              mu_alpha <- 0
              sigma_alpha <- 1
              sigma_epsilon <- 0.25
 
              make_dataset <- function(groups, group_size,
                                      mu_alpha, sigma_alpha, sigma_epsilon) {
                  replicate(
                  groups,
                  {
                   group_mean <- rnorm(1, mean = mu_alpha, sd = sigma_alpha)
                   errors <- rnorm(group_size, sd = sigma_epsilon)
                   tibble(y = group_mean + errors)
                  },
                  simplify = FALSE
                  ) %>%
                   bind_rows(.id = "group") %>%
                  mutate(group = as.factor(group))
               }
             }
  ) 

  # create a parallel version of replicate
  parReplicate <- function(cl, n, expr, simplify = TRUE, USE.NAMES = TRUE) {
  parSapply(cl, integer(n), function(i, ex) eval(ex, envir = .GlobalEnv),
            substitute(expr), simplify = simplify, USE.NAMES = USE.NAMES)
  }

  datasets <- parReplicate(
    cl,
    runs,
    make_dataset(groups, group_size, mu_alpha, sigma_alpha, sigma_epsilon),
    simplify = FALSE
  )

  stopCluster(cl)

  future::plan(multiprocess)
  results <- future_map(datasets, check_conf_int, mu_alpha = mu_alpha, .progress = TRUE)

  p <- (tibble(
    runs = seq_len(runs),
    coverage = cummean(results)
  ) %>%
    ggplot(aes(runs, coverage)) +
    geom_line() +
    geom_hline(linetype = "dashed", yintercept = 0.95) +
    coord_cartesian(ylim = c(0.93, 0.97)) 
 )

  print(p)
}

以下两个图适用于 96 个组：

如您所见，这些显示了更好的覆盖范围。最后，一个有 200 个组。

其它你可能感兴趣的问题

上一篇明显变量和观察变量之间的差异下一篇混杂因素、共线性和交互项之间有什么区别？