當報告的相關性不接近 +1/-1 時如何簡化奇異隨機結構
我在這個網站上閱讀了幾個問題的答案，為混合效應模型選擇隨機結構的最佳方法是使用理論知識。另一方面，我還閱讀了 Barr 等人 (2013) 中關於“保持最大”的建議，即通過將包括交互在內的所有固定效應作為隨機斜率來擬合最大隨機結構。然而，這似乎經常導致模型要么不會收斂，要么會收斂，但會發出“奇異擬合”的警告。在 [this] 的公認答案中（是在 +/- 1 附近沒有相關性或方差為零的奇異擬合，是誤報嗎？）問題中指出應該簡化奇異模型。
但是，當報告的相關性不接近 +/-1 並且沒有理論知識可以幫助選擇時，如何做到這一點。
一個例子將不勝感激。
Bates 等人 (2015) 概述了解決此類問題的一個好方法。
但首先有一點背景。Bates 等人 (2015) 重新分析了幾組採用最大隨機結構的實驗數據。特別是，他們重新分析了 Barr 等人 (2013) 使用的數據集，該數據集被用作“保持最大”的示例，並發現該模型嚴重過度擬合。在 Barr 等人 (2013) 中，作者擬合了一個模型，該模型具有交叉隨機效應和隨機斜率，用於兩個分組因素的 8 個固定效應。這意味著 8 個方差分量和它們之間的 28 個相關性，對於 /each/ 分組因子，總共有 72 個參數。考慮到數據只有 56 名受試者對 32 個項目做出了回應，常識應該表明這樣的模型會嚴重過度擬合。貝茨，相當外交地評估了數據將支持如此復雜的隨機結構的想法，如“
lme4在 R 中，儘管正如 Bates 所指出的那樣，這相當“不幸”，因為他們繼續表明它確實是過度擬合的，並且他們使用主成分分析來識別這一點。更新版本的 lme4 實際上使用與下面解釋的完全相同的 PCA 程序來確定模型是否已以“奇異擬合”收斂並產生警告。很多時候，這也伴隨著 +1 或 -1 的隨機效應之間的估計相關性，和/或估計為零的方差分量，但是當隨機結構很複雜（通常為 3 維或更高維）時，這些“症狀”可以缺席。在 lme4 中，在估計期間使用方差協方差 (VCV) 矩陣的 Cholesky 分解。如果 Cholesky 因子（下三角矩陣）包含一列或多列零值，則它是秩虧的，這意味著一個或多個隨機效應沒有可變性。這等效於具有沒有可變性的方差分量。PCA 是一種降維過程，當應用於估計的隨機效應 VCV 矩陣時，將立即指示該矩陣是否為滿秩。如果我們可以降低VCV矩陣的維數，即如果占方差100%的主成分個數小於VCV矩陣的列數，
因此，Bates 建議採用以下迭代過程：
- 將 PCA 應用於 VCV 矩陣以確定模型是否過度擬合（奇異）。
- 擬合一個“零相關參數”（ZCP），它將識別具有零或非常小的方差的隨機效應
- 從模型中刪除這些隨機效應並擬合新簡化的模型並檢查任何其他接近零的隨機效應。根據需要重複。
- 重新引入剩餘隨機效應之間的相關性，如果獲得非奇異擬合，則使用似然比檢驗將此模型與前一個模型進行比較。如果仍然存在奇異擬合，則返回 2。此時值得注意的是，lme4 現在在擬合過程中合併了上面的步驟 1，並將產生擬合是奇異的警告。在隨機結構簡單的模型中，例如具有單個隨機斜率的隨機截距，通常很明顯問題出在哪裡，移除隨機斜率通常可以解決問題。需要注意的是，這並不意味著總體中沒有隨機斜率，只是當前數據不支持它。
但是，當 lme4 報告擬合是奇異的，但沒有 +/- 1 的相關性或零方差分量時，事情可能會有些混亂。但是應用上述過程通常會產生一個更簡潔的模型，而不是奇異的。一個工作示例可以證明這一點：
該數據集有 3 個變量被視為固定效應：
A、B和，以及一個具有 10 個級別的C分組因子。group響應變量是Y，每組有 15 個觀察值。正如 Barr 等人 (2013) 所建議的，我們首先擬合最大模型。
> library(lme4)數據可從以下網址下載： https ://github.com/WRobertLong/Stackexchange/blob/master/data/singular.csv
在這裡，它們被加載到 R 中的 dataframe
dt中。> m0 <- lmer(y ~ A * B * C + (A * B * C | group), data = dt) boundary (singular) fit: see ?isSingular請注意，這是一個單一的擬合。但是，如果我們檢查 VCV 矩陣，我們會發現在 1 或 -1 附近沒有相關性，在 zeroL 附近也沒有任何方差分量
> VarCorr(m0) Groups Name Variance Std.Dev. Corr group (Intercept) 3.710561 1.9263 A 4.054078 2.0135 0.01 B 7.092127 2.6631 -0.01 -0.03 C 4.867372 2.2062 -0.05 -0.02 -0.22 A:B 0.047535 0.2180 -0.05 -0.47 -0.83 -0.03 A:C 0.049629 0.2228 -0.24 -0.51 0.47 -0.74 0.01 B:C 0.048732 0.2208 -0.17 0.08 -0.40 -0.77 0.50 0.44 A:B:C 0.000569 0.0239 0.24 0.43 0.37 0.65 -0.72 -0.63 -0.86 Residual 3.905752 1.9763 Number of obs: 150, groups: group, 10現在我們使用以下
rePCA函數應用 PCAlme4：> summary(rePCA(m0)) $`group` Importance of components: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] Standard deviation 1.406 1.069 1.014 0.968 0.02364 0.000853 0.00000322 0 Proportion of Variance 0.389 0.225 0.202 0.184 0.00011 0.000000 0.00000000 0 Cumulative Proportion 0.389 0.613 0.816 1.000 1.00000 1.000000 1.00000000 1這表明 VCV 矩陣有 8 列，但是秩不足，因為前 4 個主成分解釋了 100% 的方差。因此是奇異擬合，這意味著它是過擬合的，我們可以移除部分隨機結構。
所以接下來我們擬合一個“零相關參數”模型：
> m1 <- lmer(y ~ A * B * C + (A * B * C || group), data = dt) boundary (singular) fit: see ?isSingular正如我們所看到的，這也是奇異的，但是我們可以立即看到幾個方差分量現在非常接近於零：
> VarCorr(m1) Groups Name Variance Std.Dev. group (Intercept) 3.2349037958 1.7985838 group.1 A 0.9148149412 0.9564596 group.2 B 0.4766785339 0.6904191 group.3 C 1.0714133159 1.0350910 group.4 A:B 0.0000000032 0.0000565 group.5 A:C 0.0000000229 0.0001513 group.6 B:C 0.0013923672 0.0373144 group.7 A:B:C 0.0000000000 0.0000000 Residual 4.4741626418 2.1152217這些恰好是所有交互項。此外再次運行 PCA，我們再次發現 4 個組件是多餘的：
> summary(rePCA(m1)) $`group` Importance of components: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] Standard deviation 0.8503 0.4894 0.4522 0.32641 0.01764 7.152e-05 2.672e-05 0 Proportion of Variance 0.5676 0.1880 0.1605 0.08364 0.00024 0.000e+00 0.000e+00 0 Cumulative Proportion 0.5676 0.7556 0.9161 0.99976 1.00000 1.000e+00 1.000e+00 1所以現在我們從隨機結構中移除相互作用：
> m2 <- lmer(y ~ A * B * C + (A + B + C || group), data = dt)該模型現在在沒有警告的情況下收斂，並且 PCA 顯示 VCV 是滿秩的：
> summary(rePCA(m2)) $`group` Importance of components: [,1] [,2] [,3] [,4] Standard deviation 1.5436 0.50663 0.45275 0.35898 Proportion of Variance 0.8014 0.08633 0.06894 0.04334 Cumulative Proportion 0.8014 0.88772 0.95666 1.00000所以我們現在重新引入相關性：
m3 <- lmer(y ~ A * B * C + (A + B + C | group), data = dt) boundary (singular) fit: see ?isSingular…現在擬合又是奇異的，這意味著至少不需要其中一個相關性。然後我們可以繼續使用相關性較少的模型，但之前的 PCA 表明不需要 4 個組件，因此在這種情況下，我們將選擇沒有交互的模型：
Random effects: Groups Name Variance Std.Dev. group (Intercept) 10.697 3.271 group.1 A 0.920 0.959 group.2 B 0.579 0.761 group.3 C 1.152 1.073 Residual 4.489 2.119 Fixed effects: Estimate Std. Error t value (Intercept) -44.2911 30.3388 -1.46 A 12.9875 2.9378 4.42 B 13.6100 3.0910 4.40 C 13.3305 3.1316 4.26 A:B -0.3998 0.2999 -1.33 A:C -0.2964 0.2957 -1.00 B:C -0.3023 0.3143 -0.96 A:B:C 0.0349 0.0302 1.16我們還可以從固定效應估計中觀察到交互項具有相當大的標準誤差，因此在這種情況下，我們還將刪除它們，生成最終模型：
> m4 <- lmer(y ~ A + B + C + (A + B + C || group), data = dt) > summary(m4) Random effects: Groups Name Variance Std.Dev. group (Intercept) 4.794 2.189 group.1 A 0.794 0.891 group.2 B 0.553 0.744 group.3 C 1.131 1.064 Residual 4.599 2.145 Number of obs: 150, groups: group, 10 Fixed effects: Estimate Std. Error t value (Intercept) -14.000 1.868 -7.5 A 9.512 0.301 31.6 B 10.082 0.255 39.5 C 10.815 0.351 30.8我還要指出，我模擬了這個數據集，殘差和隨機截距的標準差為 2，所有隨機斜率為 1，斜率之間沒有相關性，固定截距為 -10，每個固定效應為 10 ，並且沒有交互。
因此，在這種情況下，我們選擇了一個能夠充分估計所有參數的模型。
參考：
Bates, D.、Kliegl, R.、Vasishth, S. 和 Baayen, H.，2015 年。簡約混合模型。arXiv 預印本 arXiv:1506.04967。https://arxiv.org/pdf/1506.04967.pdf
Barr, DJ, Levy, R., Scheepers, C. 和 Tily, HJ, 2013。驗證性假設檢驗的隨機效應結構：保持最大。記憶和語言雜誌，68（3），pp.255-278。