Проблема с плоскими остатками при использовании lme с varIdent

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))

Запуск иерархической модели (homoscedasti c)

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

Проверка остатков homoscedasti c модель

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

Как видно на рисунке, в остатках имеется структура гетероскедаста c по времени, поэтому я решил смоделировать ее (большая дисперсия в время 2, чем в двух других)

Создание гетероскедастий c модель

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 (heteroscedasti c) подходит лучше, чем 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