r - LME 向后选择，反向求解出现奇点

Question

我有数据，其中“飞行速度”是一个响应变量，group(实验/控制)，test(第一/第二)，FL (燃料负荷，来自瘦体重的百分比:从 0 到 ~25%)，wing(机翼长度，以毫米为单位)。由于我们对同一只鸟进行了两次测试(第一次和第二次测试，实验组被感染)，我想执行混合模型(添加一个随机项 ~1|ring)。由于异方差性，我还为 test 变量添加了 weight 参数。

mod<-lme(speed~test* group * FL * wing,weight=~1|test,random=~1|ring,data=data,method="ML")

这是完整模型的样子(我使用 nlme 包)。之后我开始向后选择。我手动执行(根据最低 AIC)，然后使用函数 stepAIC(MASS 包)检查结果。在这种情况下，前两个选择步骤很好，但是当我从模型开始时:

mod3<-lme(speed~test+group + FL + wing+ test:group + group:FL + FL:wing + test:group:wing, weight=~1|test,random=~1|ring,data=data,method="ML")

我遇到了一个错误:

 Error in MEEM(object, conLin, control$niterEM) : 
    Singularity in backsolve at level 0, block 1

据我了解，这意味着并非所有因素的相互作用都存在。但是那时我应该已经在完整模型上遇到了同样的错误。与其他响应变量一起使用效果很好。如果你们有任何想法，我会很高兴!

原始数据

ring    group   wing    speed_aver  FL  test
1   XZ13125 E   75  0.62    16.2950000  first
2   XZ13125 E   75  0.22    12.5470149  second
3   XZ13126 E   68  0.39    7.7214876   first
4   XZ13127 C   75  0.52    9.1573643   first
5   XZ13127 C   75  0.17    -1.9017391  second
6   XZ13129 C   73  0.46    10.3821705  first
7   XZ13129 C   73  0.33    -0.5278261  second
8   XZ13140 C   73  0.48    13.0774436  first
9   XZ13140 C   73  0.27    18.0092199  second
10  XZ13144 C   73  0.36    7.5144000   first
11  XZ13144 C   73  0.36    9.6820312   second
12  XZ13146 E   73  0.32    14.3651852  first
13  XZ13146 E   73  0.28    20.8171233  second
14  XZ13159 C   74  0.55    20.2760274  first
15  XZ13159 C   74  0.37    19.1687500  second
16  XZ13209 E   72  0.35    8.1464000   first
17  XZ13209 E   72  0.43    10.9945736  second
18  XZ13213 E   74  0.57    5.3682927   first
19  XZ13213 E   74  0.26    1.3584746   second
20  XZ13220 C   73  0.30    6.0105691   first
21  XZ13220 C   73  0.36    -8.0439252  second
22  XZ13230 E   74  0.44    5.3682927   first
23  XZ13230 E   74  0.31    3.0025000   second
24  XZ13231 C   75  0.28    6.2504000   first
25  XZ13231 C   75  0.37    7.7267717   second
26  XZ13232 C   74  0.34    16.8592857  first
27  XZ13232 C   74  0.33    13.7800000  second
28  XZ13271 C   73  0.32    16.2268116  first
29  XZ13271 C   73  0.28    14.3651852  second
30  XZ13278 E   72  0.45    15.5757353  first
31  XZ13278 E   72  0.37    14.9503704  second
32  XZ13280 C   74  0.33    15.0386861  first
33  XZ13280 C   74  0.36    7.6214286   second
34  XZ13340 E   73  0.62    16.8294964  first
35  XZ13340 E   73  0.26    13.7261194  second
36  XZ13367 E   75  0.42    23.4071895  first
37  XZ13370 E   71  0.25    13.6159091  first

Answer 1

事实证明这是非常棘手的。我认为问题在于，由于您构建第二个公式的方式，R 不会自动从模型矩阵中删除共线变量。

tl;dr 这有点意识流，但我认为基本的要点是

• lme 不一定会为您检查/处理模型规范中的别名(与 lm 不同，或者在较小程度上 lmer)
• 如果您违反边际，您可能会遇到 R 公式的麻烦，您在此处通过包含 test:group:wing 交互而不包含 group:wing 和 test:wing 交互。 R 允许你这样做，但模型不一定有意义......我有点惊讶你最终得到了这个模型规范——通常是 stepAIC 和 drop1 和 R 的其他内置模型简化工具，尽量尊重边际性，因此不会让你在这里结束......
• 如果您真的想要拟合这些类型的模型，请使用lmer(尽管处理异方差性更难)，或者使用构建您自己的数字虚拟变量模型.matrix() ...
• 最好使用 model.matrix() 检查这些类型的混叠问题，这超出了模型拟合的范围 (lm/lme /lmer) 函数本身 ...

为简单起见，我将省略方差模型 (weights=varIdent(form=~1|test))，因为它似乎与这个特定问题无关(我不知道先验，但使用和不使用它的测试没有区别)。

library("nlme")
form1 <- speed_aver~test* group * FL * wing
form2 <- speed_aver~test+group + FL + wing+
                      test:group + group:FL + FL:wing +
                      test:group:wing
mod <- lme(form1,random=~1|ring,data=dd,method="ML") ## OK
update(mod,form2)
## fails with "Singularity in backsolve" error

如果我们用 lme4 试试呢？

## ugh, I wish I knew a better way to append to a formula
form1L <- formula(paste(deparse(form1),"(1|ring)",sep="+"))
form2L <- formula(paste(deparse(form2,width=100),"(1|ring)",sep="+"))
library("lme4")
mod2 <- lmer(form1L, data=dd)
mod3 <- lmer(form2L, data=dd)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

啊哈! lmer 自动检测模型矩阵是否秩亏。 lm 自动执行此操作并用NA 值替换别名项。目前 lmer 只是删除了它们，尽管在 lme4 的最新版本中(已记录但未公开)选项 add.dropped=TRUE 到 fixef() 会将 NA 值放回适当的位置。

那么让我们研究一下模型矩阵:

X0 <- model.matrix(form1,data=dd)
c(rankMatrix(X0)==ncol(X0))  ## TRUE; both are 16

X <- model.matrix(form2,data=dd)
c(rankMatrix(X))==ncol(X)  ## FALSE; 11<12

尝试识别别名列:svd(X)$d 的第 12 个元素很小 (1e-15)

ss <- svd(X)
(zz <- zapsmall(ss$v[,12]))  ## elements of collinear grouping
##  [1]  0.0000000  0.0000000  0.0000000  0.0000000 -0.4472136  0.0000000
##  [7]  0.0000000  0.0000000  0.4472136  0.4472136  0.4472136  0.4472136

因此第 9-12 列的总和与第 5 列完全相同(相同的值，相反的符号)。这是怎么回事？

colnames(X)[zz!=0]
## [1] "wing" "testfirst:groupC:wing"  "testsecond:groupC:wing"
## [4] "testfirst:groupE:wing"  "testsecond:groupE:wing"

看起来我们以某种方式获得了与 wing 交互的组测试交互的所有级别，以及 wing 变量本身......

mm <- X[,zz!=0]
colnames(mm) <- gsub("(test|group|:wing)","",colnames(mm))
head(mm)
##   wing first:C second:C first:E second:E
## 1   75       0        0      75        0
## 2   75       0        0       0       75
## 3   68       0        0      68        0
## 4   75      75        0       0        0
## 5   75       0       75       0        0
## 6   73      73        0       0        0

我仍然不能 100% 确定为什么会发生这种情况，但您可以看到 R 扩展了三向交互，包括 所有 双向交互的四个级别(这又与连续的 wing 变量)，但它也有 wing

colnames(X)
##  [1] "(Intercept)"  "testsecond"    "groupE"                
##  [4] "FL"           "wing"          "testsecond:groupE"     
##  [7] "groupE:FL"    "FL:wing"       "testfirst:groupC:wing" 
## [10] "testsecond:groupC:wing" "testfirst:groupE:wing"
##      "testsecond:groupE:wing"
colnames(X0)
##  [1] "(Intercept)"               "testsecond"               
##  [3] "groupE"                    "FL"                       
##  [5] "wing"                      "testsecond:groupE"        
##  [7] "testsecond:FL"             "groupE:FL"                
##  [9] "testsecond:wing"           "groupE:wing"              
## [11] "FL:wing"                   "testsecond:groupE:FL"     
## [13] "testsecond:groupE:wing"    "testsecond:FL:wing"       
## [15] "groupE:FL:wing"            "testsecond:groupE:FL:wing"

如果我们定义一个尊重边际的模型，那么我们又没问题了......

form3 <- speed_aver~test*group*wing+FL*(group+wing)
X1 <- model.matrix(form3,dd)
c(rankMatrix(X1)== ncol(X1))  ## TRUE

我们可以通过这种方式更简单地重现问题:

form4 <- speed_aver~wing+test:group:wing
X2 <- model.matrix(form4,dd)
c(rankMatrix(X2)== ncol(X2))  ## FALSE

此模型具有三向交互(明确地)，但缺少双向交互。如果我们使用~wing*test*group，甚至~wing+wing*test*group，我们就可以了......

form5 <- speed_aver~wing+test*group*wing
X3 <- model.matrix(form5,dd)
c(rankMatrix(X3)== ncol(X3))  ## TRUE