程序员经验分享-hwhale

不同 R/lme4 版本的单一拟合结果不匹配

2024-05-17

我试图将 R 版本 3.5.3 (lme4 1.1-18-1) 的随机效应估计与 R 版本 4.1.1 (lme4 1.1-27.1) 相匹配。然而，当存在奇异拟合时，这两个版本之间的随机效应存在微小差异。我对奇点警告很满意，但令人费解的是不同版本的 R/lme4 产生的结果略有不同。

以下脚本来自 R 版本 3.5.3 (lme4 1.1-18-1) 和 R 版本 4.1.1 (lme4 1.1-27.1)，数据集来自 lme4 的拟南芥。

> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

loaded via a namespace (and not attached):
 [1] minqa_1.2.4     MASS_7.3-51.1   compiler_3.5.3  Matrix_1.2-15  
 [5] tools_3.5.3     Rcpp_1.0.1      splines_3.5.3   nlme_3.1-137   
 [9] grid_3.5.3      nloptr_1.2.1    lme4_1.1-18-1   lattice_0.20-38
> library(lme4)
Loading required package: Matrix
> options(digits = 15)
> ##########
> #Example1#
> ##########
> fit1 <- lmer(total.fruits~(1|reg)+(1|reg:popu),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> VarCorr(fit1)
 Groups   Name        Std.Dev.       
 reg:popu (Intercept)  7.744768797534
 reg      (Intercept) 10.629179104291
 Residual             39.028818969641
> ##########
> #Example2#
> ##########
> fit2 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> fit2@theta
[1] 0.150979711638631 0.000000000000000 0.189968995915902
[4] 0.260818869156072
> VarCorr(fit2)
 Groups              Name        Std.Dev.       
 reg:popu:amd:status (Intercept)  5.841181759473
 reg:popu:amd        (Intercept)  0.000000000000
 reg:popu            (Intercept)  7.349619506926
 reg                 (Intercept) 10.090696322743
 Residual                        38.688521100461
> ##########
> #Example3#
> ##########
> devfun353 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"),devFunOnly = T)
> save.image('myEnvironment353.Rdata')

> sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

loaded via a namespace (and not attached):
 [1] minqa_1.2.4        MASS_7.3-54        compiler_4.1.1     minque_2.0.0       Matrix_1.3-4      
 [6] tools_4.1.1        Rcpp_1.0.7         tinytex_0.34       splines_4.1.1      nlme_3.1-152      
[11] grid_4.1.1         xfun_0.27          nloptr_1.2.2.2     boot_1.3-28        lme4_1.1-27.1     
[16] ADDutil_2.2.1.9005 lattice_0.20-44   
> library(lme4)
Loading required package: Matrix
Warning message:
package ‘lme4’ was built under R version 4.1.2 
> options(digits = 15)
> ##########
> #Example1#
> ##########
> fit1 <- lmer(total.fruits~(1|reg)+(1|reg:popu),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
> VarCorr(fit1)
 Groups   Name        Std.Dev.       
 reg:popu (Intercept)  7.744768797534
 reg      (Intercept) 10.629179104291
 Residual             39.028818969641
> ##########
> #Example2#
> ##########
> fit2 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
boundary (singular) fit: see ?isSingular
> fit2@theta
[1] 0.150979743348540 0.000000000000000 0.189969036985684 0.260818797487214
> VarCorr(fit2)
 Groups              Name        Std.Dev.       
 reg:popu:amd:status (Intercept)  5.841182965248
 reg:popu:amd        (Intercept)  0.000000000000
 reg:popu            (Intercept)  7.349621069388
 reg                 (Intercept) 10.090693513643
 Residual                        38.688520961140
> ##########
> #Example3#
> ##########
> devfun411 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"),devFunOnly = T)
> load('myEnvironment353.Rdata')
> devfun353 <- lme4:::mkdevfun(environment(devfun353))
> minqa::bobyqa(c(1,1,1,1),devfun353,0,control = list(iprint=2))
npt = 6 , n =  4 
rhobeg =  0.2 , rhoend =  2e-07 
start par. =  1 1 1 1 fn =  6443.44054431489 
rho:    0.020 eval:  11 fn:      6393.61 par: 0.00000 0.621363 0.744867 0.823498 
rho:   0.0020 eval:  38 fn:      6361.97 par:0.156855  0.00000 0.190090 0.234676 
rho:  0.00020 eval:  49 fn:      6361.94 par:0.150719  0.00000 0.190593 0.249106 
rho:  2.0e-05 eval:  67 fn:      6361.94 par:0.150988  0.00000 0.189943 0.260821 
rho:  2.0e-06 eval:  74 fn:      6361.94 par:0.150980  0.00000 0.189965 0.260811 
rho:  2.0e-07 eval:  82 fn:      6361.94 par:0.150980  0.00000 0.189969 0.260819 
At return
eval:  90 fn:      6361.9381 par: 0.150980  0.00000 0.189969 0.260819
parameter estimates: 0.150979722854965, 0, 0.189968942342717, 0.260818725554898 
objective: 6361.93810274656 
number of function evaluations: 90 
> minqa::bobyqa(c(1,1,1,1),devfun411,0,control = list(iprint=2))
npt = 6 , n =  4 
rhobeg =  0.2 , rhoend =  2e-07 
start par. =  1 1 1 1 fn =  6443.44054431489 
rho:    0.020 eval:  11 fn:      6393.61 par: 0.00000 0.621363 0.744867 0.823498 
rho:   0.0020 eval:  38 fn:      6361.97 par:0.156855  0.00000 0.190090 0.234676 
rho:  0.00020 eval:  49 fn:      6361.94 par:0.150719  0.00000 0.190593 0.249106 
rho:  2.0e-05 eval:  67 fn:      6361.94 par:0.150988  0.00000 0.189943 0.260821 
rho:  2.0e-06 eval:  74 fn:      6361.94 par:0.150980  0.00000 0.189965 0.260811 
rho:  2.0e-07 eval:  82 fn:      6361.94 par:0.150980  0.00000 0.189969 0.260819 
At return
eval:  90 fn:      6361.9381 par: 0.150980  0.00000 0.189969 0.260819
parameter estimates: 0.150979722854965, 0, 0.189968942342717, 0.260818725554898 
objective: 6361.93810274656 
number of function evaluations: 90

当模型更简单时，没有奇点警告并且结果匹配。 （参见两个脚本中的示例 1）当模型相对复杂时，会出现奇点警告并且结果略有偏差（参见两个脚本中的示例 2）。在这种情况下，差异为

不确定这是否有用，但我也从 3.5.3 中提取 devfun 并在 4.1.1 中运行它。结果匹配。 （参见示例 3）此外，当我从 BOBYQA 读取迭代历史记录时，导致奇点警告的术语的 $\theta$ 在 0 和小数字（大约 1e-7 到 1e-9）之间振荡。

这个帖子 https://stats.stackexchange.com/questions/392302/singular-fit-in-lmer-despite-no-high-correlations-of-random-effects讨论类似的话题。它还表明奇点警告导致估计略有不同。没有明显的变化LME4 新闻 https://cran.r-project.org/web/packages/lme4/news.html造成差异的原因。This FAQ https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#singular-models-random-effect-variances-estimated-as-zero-or-correlations-estimated-as---1 and ?是单数 https://rdrr.io/cran/lme4/man/isSingular.html对奇点警告给出了很好的解释，但没有直接解决不匹配的问题。

TL;DR：有时当出现奇点警告时（我可以接受），不同 R/lme4 版本下的随机效果略有不同。为什么会发生这种情况以及如何解决？


这在一般情况下是一个很难解决的问题，在特定情况下甚至是一个相当难解决的问题。

I think版本 1.1.27.1 和 1.1.28 之间出现的差异可能来自此新闻：

交互因子的构建（例如，当 f1:f2 或 f1/f2 以随机效应项出现时）现在对于部分交叉设计更有效（不尝试创建 f1 和 f2 的所有组合）（GH #635 和 #636）

我的猜测是，这改变了 Z 矩阵中分量的顺序，这又意味着各种线性代数运算的结果不相同（例如浮点算术是不具有关联性，所以虽然二进制加法是可交换的（a + b == b + a)，从左到右计算总和可能与从右到左计算不同((a+b) + c != a + (b+c)) ...)

我尝试使用相同版本的 R（“正在开发 2022-02-25 r81818”）来重现该问题，并且仅进行比较lme4软件包版本 1.18.1 和 1.1.28.9000（开发）；任何上游包，例如Rcpp, RcppEigen, Matrix使用相同的版本。 （我不得不从开发版本中向后移植一些更改lme4到 1.1.18.1 以使其安装在最新版本的 R 下，但我认为这些修改不会影响数值结果。）

我通过安装不同版本的进行了比较lme4在新的 R 会话中运行代码之前打包。我的结果在版本 1.1.18.1 和 1.1.28 之间的差异比你的要小（两种拟合都是单一的，并且theta估计约为 2e-7 - 仍然大于您期望的 1e-8 容差，但远小于 1e-4 ...）

1.1.18.1 和 1.1.27.1 的结果相同。

  • Q1: Why are your results more different between versions than mine?
    • 一般来说/有趣的是，Windows 上的数值结果稍微不稳定/与其他平台差异更大
    • 你的两个测试平台之间的差异比我的更多：R版本，上游包（Matrix/Rcpp/RcppEigen/minqa），可能是用于构建所有内容的编译器版本和设置[所有这些都可能产生影响]
  • Q2: how should one deal with this kind of problem?
    • 作为一个次要的框架挑战，为什么（除了不理解正在发生的事情，这是一个完全合理的担忧理由）这会让您担心？结果的差异是way小于统计不确定性的大小，并且即使对于其他相同的环境（R版本，lme4，以及其他包）。
    • 您现在可以恢复到版本 1.1.27.1 ...
    • 我确实将 1.1.27.1 之间的差异视为一个 bug，至少它是包中未记录的更改。如果它的优先级足够高，我可以调查上面描述的代码更改，看看是否有一种方法可以解决它们所解决的问题而不破坏向后兼容性（理论上这应该是可能的，但可能会非常困难......）
## R CMD INSTALL ~/R/misc/lme4
library(lme4)
packageVersion("lme4")
## 1.1.18.1
fit2 <- lmer(total.fruits~(1|reg)+(1|reg:popu)+(1|reg:popu:amd)+(1|reg:popu:amd:status),data=Arabidopsis,control=lmerControl(optimizer="bobyqa"))
dput(getME(fit2, "theta"))
t1 <- c(`reg:popu:amd:status.(Intercept)` = 0.150979711638631, `reg:popu:amd.(Intercept)` = 0,
`reg:popu.(Intercept)` = 0.189968995915902, `reg.(Intercept)` = 0.260818869156072
)

在 1.1.28.9000 下运行（新的 R 会话，重新运行 package-loading/lmer上面的代码）

## R CMD INSTALL ~/R/pkgs/lme4git/lme4
packageVersion("lme4")
## [1] ‘1.1.28.9000’
dput(getME(fit2, "theta"))
t2 <- c(`reg:popu:amd:status.(Intercept)` = 0.15097974334854, `reg:popu:amd.(Intercept)` = 0,
`reg:popu.(Intercept)` = 0.189969036985684, `reg.(Intercept)` = 0.260818797487214
)

(t1-t2)/((t1+t2)/2)
## reg:popu:amd:status.(Intercept)        reg:popu:amd.(Intercept)
##                   -2.100276e-07                             NaN
##            reg:popu.(Intercept)                 reg.(Intercept)
##                  -2.161920e-07                    2.747841e-07

第二个元素是NaN因为两个版本都给出奇异拟合（0/0 == NaN）。

在 1.1.27.1 下运行（新的 R 会话，重新运行 package-loading/lmer上面的代码）

## remotes::install_version("lme4", "1.1-27.1")

t3 <- c(`reg:popu:amd:status.(Intercept)` = 0.150979711638631, `reg:popu:amd.(Intercept)` = 0,
`reg:popu.(Intercept)` = 0.189968995915902, `reg.(Intercept)` = 0.260818869156072)

identical(t1, t3) ## TRUE
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

r

lme4

mixedmodels

不同 R/lme4 版本的单一拟合结果不匹配 的相关文章

  • dplyr +“meta”-columns：当列包含要使用的其他列的名称而不是数据时

    我想知道以下问题在 dplyr 中是否有一个优雅的解决方案 要提供一个简单的可重现示例 请考虑以下 data frame df lt data frame a 1 5 b 2 6 c 3 7 ref c a a b b c stringsA

  • 为什么在 data.frame 中预先指定类型会比较慢？

    我预先分配了一个大 data frame 以便稍后填写 我通常这样做NA是这样的 n lt 1e6 a lt data frame c1 1 n c2 NA c3 NA 我想知道如果我预先指定数据类型是否会让事情变得更快 所以我测试了 f1

  • 将函数应用于矩阵列表

    我有一个矩阵列表 注意 它们的维度与此示例不同 x lt matrix 1 10 ncol 2 y lt x 300 mylist lt list x y 我想运行一个函数networklevel在矩阵列表中的每个矩阵上 该函数有各种可以计

  • ess-rdired：我收到此错误“现在没有 ESS 进程与此缓冲区关联”

    To use ess rdired为了浏览对象 我按照 ESS 手册并将以下内容添加到我的 emacs autoload ess rdired ess rdired View R objects in a dired like buffer

  • tidyverse 干扰 ggplot2 吗？无法访问map_data

    在控制台中运行这些命令 输出为 gt cty0 ggplot2 map data county gt library tidyverse Loading tidyverse ggplot2 Loading tidyverse tibble

  • 距数据帧中最近的非 NA 值的距离

    我有以下数据帧 df 我想添加一列 其中包含与每行最接近的非 NA 值的距离 df lt data frame x 1 20 df c 1 3 4 5 11 14 15 16 x lt NA 换句话说 我正在寻找以下值 df distanc

  • mlogit：需要 TRUE/FALSE 时缺少值

    我有来自离散选择实验 DCE 的数据 该实验研究了来自不同行业的个人的招聘偏好 我已经格式化为长格式 我想使用 mlogit 进行建模 我已导出数据 并且可以使用 asclogit 命令在 Stata 中成功运行模型 但在 R 中运行时遇到

  • 从受密码保护的站点读取信息

    我一直在 R 教程中使用 readLines 从网站上抓取信息 我现在希望从我自己的网站提取数据 特别是 awstats 数据 但是该域受密码保护 有没有一种方法可以通过用户名和密码传递我需要的特定 awstats 数据的 url url

  • R中使用余弦距离的层次聚类

    我想通过使用余弦相似度与 R 编程语言对文档语料库进行层次聚类 但出现以下错误 if is na n n gt 65536L stop 大小不能为 NA 或 超过 65536 需要 TRUE FALSE 时缺少值 我应该怎么办 为了重现它

  • 使用faceting()时如何连接geom_point()和geom_line？

    我有一个问题 但我在互联网上没有找到任何相关信息 我很高兴得到一些提示 我有一个数据集 其中 x 轴是离散的 但我想将这些点相互连接 我可以做到 我的问题是当我添加分面选项时 我无法再将这些点相互链接起来 我找到了一个替代方案 但看起来不太

  • RMySQL fetch - 找不到继承的方法

    使用 RMySQL 我想将数据从数据库加载到 R 中的数据帧中 为此 我使用以下代码 R连接数据库 con lt dbConnect MySQL user root password password dbname prediction h

  • 在 Windows / Linux 中创建 Mac 包

    我自己努力制作一个 r 包 我按照 stackoverflow 中上一个问题的说明进行操作如何为外行开发软件包 http cran r project org bin windows Rtools 以下是我根据上一个问题采取的步骤 在新的

  • 如何将数据从长格式重塑为宽格式

    我在重新排列以下数据框时遇到问题 set seed 45 dat1 lt data frame name rep c firstName secondName each 4 numbers rep 1 4 2 value rnorm 8 d

  • Shiny ：针对所有错误显示一条消息

    我在 R 的 Shiny 中有一个应用程序 我想处理消息 以便用户看不到发生了什么错误 我知道通过 tags style type text css shiny output error visibility hidden shiny ou

  • 基本 dyplr 函数给出错误：“check_dots_used”

    试图找出为什么我会收到此错误 以前从未见过 谷歌没有帮助 check dots used action warn 中的错误 未使用参数 action warn 我在下面的非常基本的试验中收到错误 而且在 group by count 中也收

  • 使用 R 将日期格式的字符串列表/向量转换为 posix 日期类

    我有一个日期格式的字符串列表 我想将其转换为可以使用 R 操作的 posix 日期列表 我该怎么做 这就是我所拥有的 但我最终得到了一个列表 a lt c 2009 01 01 00 00 00 2009 01 01 00 00 00 z

  • 使用矢量相应地更改传单线条的颜色

    无论如何 是否可以根据某些变量的值更改传单线条的颜色 我用谷歌搜索 发现了这个link http hgoebl github io Leaflet MultiOptionsPolyline demo 然而 我想知道是否有一种简单的方法可以在

  • R ggplot2 分面保持比率但覆盖/定义输出图大小

    我目前正在使用 ggplot2 来比较不同组的统计数据 每个组属于不同的区域 这是通过运行 R 脚本的 Web 应用程序 tikiwiki CMS 插件 R 完成的 每个区域我可以有 2 到 30 个或更多组 相同的 R 脚本针对唯一网页中

  • 使用 data.table 左连接

    假设我有两个数据表 s dataA A B 1 1 12 2 2 13 3 3 14 4 4 15 dataB A B 1 2 13 2 3 14 我有以下代码 merge test merge dataA dataB by A all d

  • SparkR 和 Sparklyr 之间导入 parquet 文件所需的时间差异

    我正在使用 databricks 导入镶木地板文件SparkR and sparklyr data1 SparkR read df dbfs data202007 source parquet header TRUE inferSchema

随机推荐

热门标签