后向选择，奇异在backsolve发生

Question

我有数据，其中“的飞行速度”是从稀体响应变量和group（实验/对照），test（第一/第二），FL（燃料负载，％质量：从0到〜25％），wing（机翼长度，单位为mm）。由于我们已经对同一只鸟进行了两次测试（第一次和第二次测试，实验组被感染），我想执行混合模型（添加一个随机词~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;博士这是一个比特流意识的-，但我认为根本拿回家点

• lme不一定检查/在你的模型规范处理走样（不像lm，或在较小程度上lmer）
• 您可以在麻烦，如果你违反边缘化，您已通过包括test:group:wing互动在这里完成，而不包括得到有R的公式和test:wing的相互作用。 R可以让你做到这一点，但是模型并不一定是有意义的......我有点惊讶，你最终得到了这个模型规范 - 通常是stepAIC和drop1，以及R的其他内置模型简化工具，尽量尊重边缘化，因而不会让你在这里结束...
• 如果你真的想适合这样的模式，使用lmer（虽然处理异方差性是很难），或者构建自己的数字虚变量与model.matrix() ...
• 检查出这些种类的混叠问题最好用model.matrix()完成，超出了模型拟合的范围（lm/lme/lmer）功能本身...

为简单起见，我要离开了方差模型（weights=varIdent(form=~1|test)），因为它似乎并没有被相关与这个特定问题（我不知道先验，但测试与和没有它没有不同）。

​​3210

如果我们用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变量本身沿水平所有 ...

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