为什么这个简单的混合模型不能收敛？

机器算法验证 r 混合模式 lme4-nlme 重复测量 统计能力
2022-04-07 20:56:23

我有数据，我想比较第value28 天和第 83 天之间变量的平均值。因为经验涉及伪复制 ( culture)，我正在考虑使用具有随机效应的混合模型， culture我有两个问题。第一个是关于这个数据集的：

  library(lme4)
#> Loading required package: Matrix
  library(lmerTest)
#> 
#> Attaching package: 'lmerTest'
#> The following object is masked from 'package:lme4':
#> 
#>     lmer
#> The following object is masked from 'package:stats':
#> 
#>     step

df <- structure(list(
  day_true = c(28, 83, 28, 83, 28, 83), value = c(
    758453.333333333,
    575133.333333333, 684160, 656933.333333333, 816840, 734700
  ),
  culture = c(1L, 1L, 2L, 2L, 3L, 3L)
), row.names = c(NA, -6L), class = c("data.frame"))

如果我按如下方式安装它，我没有警告。

lme4::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 102.7974
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 47535   
#>  Residual             55990   
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)83  
#>             753151              -97562

但是，如果我适合它，lmerTest我会收到警告。我应该关心它吗？

lmerTest::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
#> eigenvalue close to zero: 2.6e-09
#> Linear mixed model fit by REML ['lmerModLmerTest']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 102.7974
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 47535   
#>  Residual             55990   
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)83  
#>             753151              -97562

根据得到的模型，没有显着差异。但在视觉上它似乎有一个。我只是想确保了解收敛问题。

boxplot(value ~ factor(day_true), data = df)

我的第二个问题是关于这些数据的，对此我有一个独特的适合 信息。我找不到原因。是因为我的分数很少（每组 n = 3）吗？或者，是否有更好的分析来比较这两组之间重复测量的平均值？

df <- structure(list(day_true = c(0, 28, 0, 28, 0, 28), value = c(
  34.6732447526395,
  31.5635584852635, 34.2763264775584, 32.1719125021771, 35.0747566289866,
  31.7318622838194
), culture = c(1L, 1L, 2L, 2L, 3L, 3L)), row.names = c(
  NA,
  -6L
), class = c("data.frame"))

  lme4::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> singular fit
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 5.3578
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 0.0000  
#>  Residual             0.3592  
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)28  
#>             34.675              -2.852  
#> convergence code 0; 1 optimizer warnings; 0 lme4 warnings

由reprex 包（v0.2.1）于 2019 年 2 月 6 日创建

会话信息

devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.5.2 (2018-12-20)
#>  os       Linux Mint 19.1             
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language en_CA:en                    
#>  collate  en_CA.UTF-8                 
#>  ctype    en_CA.UTF-8                 
#>  tz       America/Toronto             
#>  date     2019-02-06                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package     * version    date       lib source                          
#>  assertthat    0.2.0      2017-04-11 [1] CRAN (R 3.5.1)                  
#>  backports     1.1.3      2018-12-14 [1] CRAN (R 3.5.1)                  
#>  callr         3.1.1      2018-12-21 [1] CRAN (R 3.5.2)                  
#>  cli           1.0.1      2018-09-25 [1] CRAN (R 3.5.1)                  
#>  colorspace    1.4-0      2019-01-13 [1] CRAN (R 3.5.2)                  
#>  crayon        1.3.4      2017-09-16 [1] CRAN (R 3.5.1)                  
#>  curl          3.3        2019-01-10 [1] CRAN (R 3.5.2)                  
#>  desc          1.2.0      2018-05-01 [1] CRAN (R 3.5.1)                  
#>  devtools      2.0.1      2018-10-26 [1] CRAN (R 3.5.1)                  
#>  digest        0.6.18     2018-10-10 [1] CRAN (R 3.5.1)                  
#>  dplyr         0.8.0      2019-02-01 [1] Github (tidyverse/dplyr@e9902f1)
#>  evaluate      0.12       2018-10-09 [1] CRAN (R 3.5.1)                  
#>  fs            1.2.6      2018-08-23 [1] CRAN (R 3.5.1)                  
#>  ggplot2       3.1.0      2018-10-25 [1] CRAN (R 3.5.1)                  
#>  glue          1.3.0      2018-07-17 [1] CRAN (R 3.5.1)                  
#>  gtable        0.2.0      2016-02-26 [1] CRAN (R 3.5.1)                  
#>  highr         0.7        2018-06-09 [1] CRAN (R 3.5.1)                  
#>  htmltools     0.3.6      2017-04-28 [1] CRAN (R 3.5.1)                  
#>  httr          1.4.0      2018-12-11 [1] CRAN (R 3.5.1)                  
#>  knitr         1.21       2018-12-10 [1] CRAN (R 3.5.1)                  
#>  lattice       0.20-38    2018-11-04 [1] CRAN (R 3.5.1)                  
#>  lazyeval      0.2.1      2017-10-29 [1] CRAN (R 3.5.1)                  
#>  lme4        * 1.1-20     2019-02-04 [1] CRAN (R 3.5.2)                  
#>  lmerTest    * 3.0-1      2018-04-23 [1] CRAN (R 3.5.2)                  
#>  magrittr      1.5        2014-11-22 [1] CRAN (R 3.5.1)                  
#>  MASS          7.3-51.1   2018-11-01 [1] CRAN (R 3.5.1)                  
#>  Matrix      * 1.2-15     2018-11-01 [1] CRAN (R 3.5.1)                  
#>  memoise       1.1.0      2017-04-21 [1] CRAN (R 3.5.1)                  
#>  mime          0.6        2018-10-05 [1] CRAN (R 3.5.1)                  
#>  minqa         1.2.4      2014-10-09 [1] CRAN (R 3.5.1)                  
#>  munsell       0.5.0      2018-06-12 [1] CRAN (R 3.5.1)                  
#>  nlme          3.1-137    2018-04-07 [1] CRAN (R 3.5.2)                  
#>  nloptr        1.2.1      2018-10-03 [1] CRAN (R 3.5.1)                  
#>  numDeriv      2016.8-1   2016-08-27 [1] CRAN (R 3.5.1)                  
#>  pillar        1.3.1.9000 2019-02-01 [1] Github (r-lib/pillar@3a54b8d)   
#>  pkgbuild      1.0.2      2018-10-16 [1] CRAN (R 3.5.1)                  
#>  pkgconfig     2.0.2      2018-08-16 [1] CRAN (R 3.5.1)                  
#>  pkgload       1.0.2      2018-10-29 [1] CRAN (R 3.5.1)                  
#>  plyr          1.8.4      2016-06-08 [1] CRAN (R 3.5.1)                  
#>  prettyunits   1.0.2      2015-07-13 [1] CRAN (R 3.5.1)                  
#>  processx      3.2.1      2018-12-05 [1] CRAN (R 3.5.1)                  
#>  ps            1.3.0      2018-12-21 [1] CRAN (R 3.5.2)                  
#>  purrr         0.3.0      2019-01-27 [1] CRAN (R 3.5.2)                  
#>  R6            2.3.0      2018-10-04 [1] CRAN (R 3.5.1)                  
#>  Rcpp          1.0.0      2018-11-07 [1] CRAN (R 3.5.1)                  
#>  remotes       2.0.2      2018-10-30 [1] CRAN (R 3.5.1)                  
#>  rlang         0.3.1.9000 2019-02-01 [1] Github (tidyverse/rlang@7243c6d)
#>  rmarkdown     1.11       2018-12-08 [1] CRAN (R 3.5.1)                  
#>  rprojroot     1.3-2      2018-01-03 [1] CRAN (R 3.5.1)                  
#>  scales        1.0.0      2018-08-09 [1] CRAN (R 3.5.1)                  
#>  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 3.5.1)                  
#>  stringi       1.2.4      2018-07-20 [1] CRAN (R 3.5.1)                  
#>  stringr       1.3.1      2018-05-10 [1] CRAN (R 3.5.1)                  
#>  testthat      2.0.1      2018-10-13 [1] CRAN (R 3.5.1)                  
#>  tibble        2.0.1      2019-01-12 [1] CRAN (R 3.5.2)                  
#>  tidyselect    0.2.5      2018-10-11 [1] CRAN (R 3.5.1)                  
#>  usethis       1.4.0      2018-08-14 [1] CRAN (R 3.5.1)                  
#>  withr         2.1.2      2018-03-15 [1] CRAN (R 3.5.1)                  
#>  xfun          0.4        2018-10-23 [1] CRAN (R 3.5.1)                  
#>  xml2          1.2.0      2018-01-24 [1] CRAN (R 3.5.1)                  
#>  yaml          2.2.0      2018-07-25 [1] CRAN (R 3.5.1)                  
#> 
#> [1] /home/pmassicotte/R/x86_64-pc-linux-gnu-library/3.5
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library

2个回答

这很可能不是您需要担心的警告。如您所见，两种情况下的参数估计值是相同的。lmerin的版本lmertest显然比当前版本具有更保守的收敛性检查lme4。

中的问题lmertest::lmer是由于变量的尺度差异很大，这可能使一些优化例程不稳定，从而产生警告。如果重新缩放，问题就会消失：

> df$value <- df$value / 10000

> lmerTest::lmer(value ~ factor(day_true) + (1 | culture), data = df)

Linear mixed model fit by REML ['lmerModLmerTest']
Formula: value ~ factor(day_true) + (1 | culture)
   Data: df
REML criterion at convergence: 29.1147
Random effects:
 Groups   Name        Std.Dev.
 culture  (Intercept) 4.753   
 Residual             5.599   
Number of obs: 6, groups:  culture, 3
Fixed Effects:
       (Intercept)  factor(day_true)83  
            75.315              -9.756

请注意，由于您只有 3 个组，因此最好建模culture为固定：

> lm(value ~ factor(day_true) + culture, data = df)
Coefficients:
   (Intercept)  factor(day_true)83             culture  
        64.417              -9.756               5.449

day_true显然，在两种分析中对 的固定效应的估计是相同的。

之所以没有找到统计上显着的估计值，是因为样本量太小。在收集数据和拟合模型之前运行“功率分析”是非常可取的。

我不能谈论计算问题，但是六个数据点根本不是很多，如果你想拟合随机效应就更不用说了——每个数据点只有 2 个，每个只有 3 个例子。

无统​​计学意义的结果在这里很有意义：您有少量数据暗示某些事情，但不足以计算出足够好的任何东西以得出任何一般性结论。

您的箱形图误导了您：它是计算均值、四分位数等，每个框只有 3 个数字。在这种情况下，可视化特征多于数据。

其它你可能感兴趣的问题
上一篇有没有不给出概率的分类算法？ 下一篇使用 R 的二次模型。使用 poly(..) 和 I(..) 函数（R 语言）