Предупреждение lme4: модели не удалось сойтись с max | grad |

Question

Мне нужно запустить lmer с измененной логарифмической переменной ответа, непрерывной переменной как фиксированный эффект и вложенным случайным эффектом:

first<-lmer(logterrisize~spm + (1|studyarea/teriid),
            data = Data_table_for_analysis_Character_studyarea,
            control=lmerControl(optimizer="Nelder_Mead",
                                 optCtrl=list(maxfun=1e4)))

Я получил это сообщение об ошибке: Ошибка в длине (значение <- as.numeric (value)) == 1L: Downdated VtV не является положительно определенным </p>

Я пробовал это с bobyqa () в качестве аргумента оптимизации и получил это предупреждение:

1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : 
Model failed to converge with max|grad| = 0.753065 (tol = 0.002, component
1) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,:
Model is nearly unidentifiable: very large eigenvalue-Rescale variables?

myСводка выглядит следующим образом:

Linear mixed model fit by REML ['lmerMod'] 
Formula: logterrisize ~ spm + (1 studyarea/teriid) Data: Data_table_for_analysis_Character_studyareaControl: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000)) REML criterion at convergence: -6079.6Scaled residuals: 
   Min         1Q     Median         3Q        Max 
-3.639e-07 -4.962e-08  3.310e-09  5.293e-08  9.725e-07 
Random effects:
 Groups           Name        Variance  Std.Dev.
 teriid:studyarea (Intercept) 1.291e-01 3.593e-01
 studyarea        (Intercept) 1.944e-02 1.394e-01
 Residual                     4.506e-15 6.712e-08
Number of obs: 273, groups:  teriid:studyarea, 66; studyarea, 22
Fixed effects:
 Estimate Std. Error t value
(Intercept)  1.480e+00  5.631e-02   26.28
spm         -5.785e-16  8.507e-10    0.00
Correlation of Fixed Effects:
 (Intr) spm 0.000  convergence code: 0
Model failed to converge with max|grad| = 0.753065 (tol = 0.002, component1)
Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

мои данные выглядят следующим образом:

show(logterrisize) [1] 1.3317643 1.3317643 1.3317643 0.1295798 0.1295798 1.5051368 1.5051368 1.5051368 1.5051368 [10] 1.5051368 1.5051368 1.5051368 1.5051368 1.5051368 1.5051368 1.5051368 1.5051368 1.5051368 [19] 1.5051368 1.5051368 1.5051368 1.4665993 1.4665993 1.4665993 1.8282328 1.8282328 1.9252934 [28] 1.9252934 1.9252934 2.3006582 2.3006582 2.5160920 2.7774040 2.7774040 3.3398623 3.3398623 [37] 3.4759297 1.2563594 1.6061204 1.6061204 1.7835139 1.7835139 2.1669498 2.1669498 2.1669498 [46] 2.1669498 0.7264997 0.7458155 0.8380524 0.8380524 0.8380524 0.8380524 0.8380524 0.8380524

 show(spm)  [1] 18.461538 22.641509 35.172414 10.418006 15.611285  3.482143  3.692308  4.483986  4.821429 [10]  6.000000  6.122449  6.176471  6.220736  6.260870  6.593407  7.010309  9.200000  9.473684 [19]  9.600000 12.600000 14.200000 16.146179 28.125000 30.099010 13.731343 14.432990 11.089109 [28] 17.960526 32.903226  8.955224 33.311688  8.800000 11.578947 20.000000 14.455446 18.181818 [37] 28.064516 25.684211 17.866667 23.142857 18.208955 20.536913 11.419355 11.593220 12.703583 [46] 20.000000  3.600000 11.320755  6.200000  6.575342 12.800000 19.109589 20.124224 22.941176 [55]  4.600000  6.600000  6.771160  8.000000 19.200000 19.400000 22.773723  3.333333  4.214047

Studyarea - это имена персонажей, а teriID представляет непрерывные номера сайтов исследования.

Myфрейм данных выглядит следующим образом:

Я забыл что-нибудь включить в уравнение при использовании преобразованной лог-переменной?Спасибо!

РЕДАКТИРОВАТЬ: я использовал? Схождение, чтобы проверить на ошибки сходимости.Я попробовал это:

## 3. пересчитать градиент и Гессиана с экстраполяцией Ричардсона

devfun <- update(first, devFunOnly=TRUE)
if (isLMM(first)) {
pars <- getME(first,"theta")
} else {## GLMM: requires both random and fixed parameters
pars <- getME(first, c("theta","fixef"))
}
if (require("numDeriv")) {
 cat("hess:\n"); print(hess <- hessian(devfun, unlist(pars)))
cat("grad:\n"); print(grad <- grad(devfun, unlist(pars)))
cat("scaled gradient:\n")
print(scgrad <- solve(chol(hess), grad))}

и получил ответ:

hess:
      [,1]      [,2]
[1,] 147.59157 -14.37956
[2,] -14.37956 120.85329
grad:
[1] -222.1020 -108.1038
scaled gradient:
[1] -19.245584  -9.891077

К сожалению, я не знаючто ответ должен сказать мне.

2-е РЕДАКТИРОВАНИЕ:

Я пробовал многочисленные оптимизаторы и при использовании этого:

first<-lmer(logterrisize~spm + (1|studyarea/teriid),REML=FALSE,
            data = Data_table_for_analysis_Character_studyarea,
            control=lmerControl(optimizer="optimx",
                                 optCtrl=list(method='nlminb')))

Я получил только одно предупреждение: In optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp), : convergence code 1 from optimx

теперь мое резюме выглядит так:

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: logterrisize ~ spm + (1 | studyarea/teriid)
Data: Data_table_for_analysis_Character_studyarea
Control: lmerControl(optimizer = "optimx", optCtrl = list(method ="nlminb"))
AIC      BIC   logLik deviance df.resid 
 -3772.4  -3754.3   1891.2  -3782.4      268 
Scaled residuals: 
   Min         1Q     Median         3Q        Max 
-1.523e-04 -1.693e-05  1.480e-06  1.436e-05  3.332e-04 
Random effects:
Groups           Name        Variance  Std.Dev. 
teriid:studyarea (Intercept) 8.219e-02 0.2866882
studyarea        (Intercept) 7.478e-02 0.2734675
Residual                     3.843e-10 0.0000196
Number of obs: 273, groups:  teriid:studyarea, 66; studyarea, 22
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.551e+00  7.189e-02   21.58
spm         3.210e-11  2.485e-07    0.00
Correlation of Fixed Effects:
(Intr)spm 0.000 
convergence code: 1

Так могу ли я закрыть глаза на это предупреждающее сообщение или это будет огромная ошибка?

Answer 1

tl; dr каждое наблюдение на территории разделяет один и тот же размер территории, поэтому ваш случайный эффект идентификатора территории по сути объясняет все и не оставляет никаких изменений ни для фиксированного эффекта log(terrsize), ни дляостаточная дисперсия.Исключение случайного эффекта идентификатора территории из модели дает разумные ответы;смоделированный набор данных довольно хорошо копирует эту патологию, но предполагает, что в результате вы получите заниженную оценку эффекта spm ...

чтения данных и построения графика

library(readxl)
library(dplyr)

dd <- (read_excel("lme4_terr_dataset.xlsx")
    %>% rename(spm="scans per min",
               studyarea="Study areaID",
               teriid="TerritoryID",
               terrsize="Territory_Size")
)

library(ggplot2); theme_set(theme_bw())
library(ggalt)
(ggplot(dd, aes(spm,terrsize,colour=studyarea))
    +geom_point()
    +geom_encircle(aes(group=teriid))
    +theme(legend.position="none")
    + scale_y_log10()
)

Этот график, с его горизонтальными линиями значений из того же идентификатора территории, помог мне диагностировать проблему.Подтверждение того, что у каждого идентификатора территории есть единый размер территории для всех наблюдений:

tt <- with(dd,table(terrsize,teriid))
all(rowSums(tt>0)==1)  ## TRUE

подбор модели

library(lme4)
m1 <- lmer(log(terrsize) ~ spm + (1|studyarea/teriid), dd)
## replicate warnings    
m2 <- lmer(log(terrsize) ~ spm + (1|studyarea), dd)
## no warnings

теперь моделирует похожие данные

set.seed(101)
## experimental design: rep within f2 (terr_id) within f1 (study area)
ddx <- expand.grid(studyarea=factor(letters[1:10]),
                   teriid=factor(1:4),rep=1:5)
## study-area, terr_id effects, and spm
b_studyarea <- rnorm(10)
b_teriid <- rnorm(40)
ddx <- within(ddx, {
    int <- interaction(studyarea,teriid)
    spm <- rlnorm(nrow(ddx), meanlog=1,sdlog=1)
})
## compute average spm per terr/id
## (because response will be identical across id)
spm_terr <- aggregate(spm~int, data=ddx, FUN=mean)[,"spm"]
ddx <- within(ddx, {
    mu <- 1+0.2*spm_terr[int]+b_studyarea[studyarea] + b_teriid[int]
    tsize <- rlnorm(length(levels(int)), meanlog=mu, sdlog=1)
    terrsize <- tsize[int]
})
gg1 %+% ddx

соответствует моделируемым данным

Это дает поведение, подобное реальным данным:

lmer(log(terrsize) ~ spm + (1|studyarea/teriid), ddx)

Мы можем избежать предупреждений, сбросив teriid:

m1 <- lmer(log(terrsize) ~ spm + (1|studyarea), ddx)

Но истинный эффект spm (0,2) будет недооценен (из-за игнорируемого шума от teriid ...)

round(confint(m1, parm="beta_"),3)
##             2.5 % 97.5 %
## (Intercept) 1.045  2.026
## spm         0.000  0.070

агрегирование

На основе этого единственного моделирования это выглядит как агрегирование на уровне территории (как рекомендовано, например, в Murtaugh 2007, "Простота и сложность в анализе экологических данных" Экология ) и взвешиваниепо количеству образцов на территорию дает разумную оценку истинного spm эффекта ...

ddx_agg <- (ddx
    %>% group_by(studyarea,terrsize,teriid)
    %>% summarise(spm=mean(spm),
                  n=n())
)
library(nlme)
m3x <- lme(log(terrsize) ~ spm, random=~1|studyarea, data=ddx_agg,
          weights=varFixed(~I(1/n)))
round(summary(m3x)$tTab,3)
            Value Std.Error DF t-value p-value
(Intercept) 0.934     0.465 29   2.010   0.054
spm         0.177     0.095 29   1.863   0.073