r - 使用带有 varIdent 的 lme 的平面残差问题

Question

我在运行带有 nlme 包的分层（混合）模型时遇到了一个奇怪的问题。当沿一个因子添加异质方差时，残差在其中一个变得平坦。

这是代码：

创建数据框

x <- data.frame("time"    = factor(c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3)), 
                "rep"     =        c(1,1,1,1,2,2,2,2,3,3,3,3,1,1,1,1,2,2,2,2,3,3,3,3,1,1,1,1,2,2,2,2,3,3,3,3),
                "varB"    = factor(c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2)), 
                "site"    = factor(c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2)),
                "varA"    = factor(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2)),
                "subsite" = factor(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2)),
                "response"= c(0.657, 0.634, 0.723, 0.410, 0.649, 0.820, 0.382, 0.394, 0.586, 0.572, 0.603,
                              0.441, 0.563, 0.514, 1.120, 1.011, 0.449, 0.822, 0.209, 0.559, 0.440, 0.374, 
                              0.511, 0.640, 0.599, 0.245, 0.582, 0.186, 0.407, 0.271, 0.351, 0.202, 0.154, 
                              0.328, 0.220, 0.233))

运行分层模型（同方差）

library(nlme)
m1 <- lme(response ~ (varA+time+varB)^2, 
            random = ~1|rep/site/subsite, 
            data = x, 
            na.action = na.omit, 
            method = "REML")

检查同方差模型的残差

plot(m1, resid(., type = "p") ~ fitted(.) | time, 
     abline = c(-2,0,2), main = "Residuals by time", id=0.05)

从图中可以看出，残差中有一个随时间变化的异方差结构，所以我决定对其进行建模（时间 2 的离散度大于其他两个）

创建异方差模型

m2 <- update(m1, weights = varIdent(form = ~1|time), control = lmeControl(niterEM = 1000))

比较模型

anova(m1,m2)

   Model df      AIC      BIC   logLik   Test  L.Ratio p-value
m1     1 14 25.06879 42.68214 1.465604                        
m2     2 16 20.96663 41.09617 5.516685 1 vs 2 8.102163  0.0174

模型 2（异方差）比 1 更适合

当我检查新模型的残差时，问题就出现了：

plot(m2, resid(., type = "p") ~ fitted(.) | time, 
     abline = c(-2,0,2), main = "Residuals by time", id=0.05)

从图中可以看出，时间 1 的残差变平了......可能是什么问题？

我不确定它是否有帮助，但也许检查模型 2 的方差结构可能是找出问题的线索：

summary(m2$modelStruct$varStruct)

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | time 
 Parameter estimates:
       1        2        3 
   1.000 4557.519 3274.876 

标签： rhierarchical-datamixed-modelsnlmeheterogeneous