我有数据，其中“飞行速度”是响应变量，group（实验/控制），test（第一/秒），FL（燃料负荷，瘦体重百分比：从0到〜 25％），wing（以毫米为单位的机翼长度）。由于我们已经对同一只鸟进行了两次测试（第一次和第二次测试，实验组被感染），因此我想执行混合模型（添加随机项~1|ring）。由于异方差性，我还为weight变量添加了test参数。
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
最佳答案
事实证明，这非常棘手。我认为问题在于，由于构造第二个公式的方式，R不会自动从模型矩阵中删除共线变量。
tl; dr这有点意识流，但是我认为基本要点是lme不一定为您检查/处理模型规范中的别名（与lm不同，或在较小程度上为lmer）
如果您违反了边际性，那么您可能会遇到R公式的麻烦，您可以通过在不包含test:group:wing和group:wing交互的情况下包含test:wing交互来完成此操作。 R让您可以做到这一点，但是该模型不一定有意义……让您惊讶的是您最终获得了此模型规范-通常是stepAIC和drop1，以及R的其他内置函数模型简化工具，请尝试尊重边缘性，因此不会让您到这里...
如果您真的想适合这类模型，请使用lmer（尽管处理异方差会更困难），或者使用model.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变量交互），但是也得到了wingcolnames(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
关于r - 在LME中向后选择，在反向求解中出现奇点，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/26449969/