[解決済み] LMEにおける後方選択、バックソルブにおける特異点の発生

質問

飛行速度を応答変数とするデータがあり group (実験/対照)です。 test (1回目/2回目)となります。 FL (燃料負荷、除脂肪体重からの％：0～25％）。 wing (翼長（mm））。同じ鳥を2回試験しているので（1回目、2回目、実験群は感染）、混合モデルを実行したい（ランダム項 ~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パッケージ)。この場合、最初の2つのステップの選択はうまくいっていますが、モデルから始めると、このようになります。

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

解決方法は？

これは、結果的にかなり厄介なことです。 I 思う 問題は、2つ目の式を構築する方法によって、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

エイリアスされたカラムを識別しようとする。の12番目の要素。 svd(X)$d が小さい (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 と相互作用する test-by-group 交互作用の水準。 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が三者間相互作用を拡大するのは、次のような場合であることがわかるでしょう。 すべて 双方向交互作用の4つのレベル（これらは順番に連続的な 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"

マージンを尊重したモデルを定義すれば、またOKなのですが．

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 ということであれば、OKでしょう。

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