不同 R/lme4 版本的奇异拟合结果不匹配
Mismatching results for singular fit with different R/lme4 versions
我正在尝试匹配 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）。在这种情况下，差异是 <1e-5，但我之前已经观察到 <1e-4。任何人都可以阐明为什么结果略有不同吗？甚至有可能将结果匹配到至少 1e-8 吗？
不确定这是否有用，但我还从 3.5.3 中提取了 devfun，并在 4.1.1 中 运行 提取了它。结果一致。 （参见示例 3）此外，当我从 BOBYQA 读取迭代历史记录时，导致奇点警告的项的 $\theta$ 在 0 和小数（大约 1e-7 到 1e-9）之间振荡。
This post discusses similar topics. It also shows the singularity warning leads to slightly different estimate. There is no obvious change in LME4 NEWS that cause the difference. This FAQ and ?isSingular对奇点警告给出了很好的解释，但没有直接解决不匹配的问题。
TL;DR: 有时出现奇点警告时（我没意见），不同 R/lme4 版本下的随机效果略有不同。为什么会发生这种情况以及如何解决？
这是一个总体上很难解决的问题，在特定情况下甚至是比较难解决的问题。
我认为 版本 1.1.27.1 和 1.1.28 之间存在差异，可能来自此新闻条目：
construction of interacting factors (e.g. when f1:f2 or f1/f2 occur in random effects terms) is now more efficient for partially crossed designs (doesn't try to create all combinations of f1 and f2) (GH #635 and #636)
我的猜测是，这会改变 Z 矩阵中分量的顺序，这反过来意味着各种线性代数运算的结果并不相同（例如，浮点运算 不相关，因此虽然二进制加法是可交换的 (a + b == b + a)，但求和的 left-to-right 可能与 right-to-left 求值 ((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 下，但我认为这些修改中的任何一个都不会影响数值结果。 )
我通过在 运行 新 R 会话中的代码之前安装不同版本的 lme4 包来进行比较。我的结果在版本 1.1.18.1 和 1.1.28 之间的差异比你的小（两个拟合都是单一的，theta 估计的相对差异是 2e-7 的数量级——仍然大于你想要的 1e- 8 公差但比 1e-4 小得多 ...)
1.1.18.1 和 1.1.27.1 的结果相同。
- Q1：为什么你的版本之间的结果比我的差异更大？
- 在 general/anecdotally 中，Windows 上的数值结果略多 unstable/differ 来自其他平台
- 你的两个测试平台之间的差异比我的更多：R 版本、上游包 (
Matrix/Rcpp/RcppEigen/minqa)，可能用于构建所有内容的编译器版本和设置[所有这些都可能有所作为]
- Q2：遇到这种问题应该怎么处理？
- 作为一个小的框架挑战，为什么（除了不了解发生了什么，这是一个完全合理的关注理由）你担心吗？结果的差异 小于统计不确定性的大小，而且这么大的差异也可能发生在不同的平台上（OS/compiler version/etc。）即使对于其他相同的环境（R 版本、
lme4 和其他包）。
- 您现在可以恢复到版本 1.1.27.1 ...
- 我确实将 1.1.27.1 之间的差异视为某种错误 — 至少它是程序包中未记录的更改。如果足够 high-priority 我可以调查上面描述的代码更改，看看是否有办法在不破坏向后兼容性的情况下解决他们解决的问题（理论上这应该是可能的，但它可能非常困难.. .)
## 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
)
运行 under 1.1.28.9000 (fresh R session, re-run 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)。
运行 under 1.1.27.1（新鲜的R会话，re-run 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
我正在尝试匹配 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）。在这种情况下，差异是 <1e-5，但我之前已经观察到 <1e-4。任何人都可以阐明为什么结果略有不同吗？甚至有可能将结果匹配到至少 1e-8 吗？
不确定这是否有用，但我还从 3.5.3 中提取了 devfun，并在 4.1.1 中 运行 提取了它。结果一致。 （参见示例 3）此外，当我从 BOBYQA 读取迭代历史记录时，导致奇点警告的项的 $\theta$ 在 0 和小数（大约 1e-7 到 1e-9）之间振荡。
This post discusses similar topics. It also shows the singularity warning leads to slightly different estimate. There is no obvious change in LME4 NEWS that cause the difference. This FAQ and ?isSingular对奇点警告给出了很好的解释，但没有直接解决不匹配的问题。
TL;DR: 有时出现奇点警告时（我没意见），不同 R/lme4 版本下的随机效果略有不同。为什么会发生这种情况以及如何解决？
这是一个总体上很难解决的问题，在特定情况下甚至是比较难解决的问题。
我认为 版本 1.1.27.1 和 1.1.28 之间存在差异，可能来自此新闻条目：
construction of interacting factors (e.g. when f1:f2 or f1/f2 occur in random effects terms) is now more efficient for partially crossed designs (doesn't try to create all combinations of f1 and f2) (GH #635 and #636)
我的猜测是，这会改变 Z 矩阵中分量的顺序，这反过来意味着各种线性代数运算的结果并不相同（例如，浮点运算 不相关，因此虽然二进制加法是可交换的 (a + b == b + a)，但求和的 left-to-right 可能与 right-to-left 求值 ((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 下，但我认为这些修改中的任何一个都不会影响数值结果。 )
我通过在 运行 新 R 会话中的代码之前安装不同版本的 lme4 包来进行比较。我的结果在版本 1.1.18.1 和 1.1.28 之间的差异比你的小（两个拟合都是单一的，theta 估计的相对差异是 2e-7 的数量级——仍然大于你想要的 1e- 8 公差但比 1e-4 小得多 ...)
1.1.18.1 和 1.1.27.1 的结果相同。
- Q1：为什么你的版本之间的结果比我的差异更大？
- 在 general/anecdotally 中，Windows 上的数值结果略多 unstable/differ 来自其他平台
- 你的两个测试平台之间的差异比我的更多：R 版本、上游包 (
Matrix/Rcpp/RcppEigen/minqa)，可能用于构建所有内容的编译器版本和设置[所有这些都可能有所作为]
- Q2：遇到这种问题应该怎么处理？
- 作为一个小的框架挑战，为什么（除了不了解发生了什么，这是一个完全合理的关注理由）你担心吗？结果的差异 小于统计不确定性的大小，而且这么大的差异也可能发生在不同的平台上（OS/compiler version/etc。）即使对于其他相同的环境（R 版本、
lme4和其他包）。 - 您现在可以恢复到版本 1.1.27.1 ...
- 我确实将 1.1.27.1 之间的差异视为某种错误 — 至少它是程序包中未记录的更改。如果足够 high-priority 我可以调查上面描述的代码更改，看看是否有办法在不破坏向后兼容性的情况下解决他们解决的问题（理论上这应该是可能的，但它可能非常困难.. .)
- 作为一个小的框架挑战，为什么（除了不了解发生了什么，这是一个完全合理的关注理由）你担心吗？结果的差异 小于统计不确定性的大小，而且这么大的差异也可能发生在不同的平台上（OS/compiler version/etc。）即使对于其他相同的环境（R 版本、
## 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
)
运行 under 1.1.28.9000 (fresh R session, re-run 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)。
运行 under 1.1.27.1（新鲜的R会话，re-run 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