後向選擇，奇異在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