Aviso "O modelo falhou ao convergir" no lmer ()

Com o seguinte conjunto de dados, eu queria ver se a resposta (efeito) muda em relação a sites, temporada, duração e suas interações. Alguns fóruns on-line sobre estatísticas sugeriram que eu continuasse com os Modelos de efeitos mistos lineares, mas o problema é que, como as réplicas são randomizadas em cada estação, tenho pouca chance de coletar a amostra exatamente do mesmo local em estações sucessivas (por exemplo, repl-1 de s1 de pós-monção pode não ser o mesmo que o de monção). É diferente dos ensaios clínicos (com design dentro do assunto) em que você mede o mesmo assunto repetidamente ao longo das estações. No entanto, considerando sites e temporada como um fator aleatório, executei os seguintes comandos e recebi uma mensagem de aviso:

Warning messages:
1: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: Model failed to converge: degenerate Hessian with 1 negative eigenvalues

Alguém pode me ajudar a resolver o problema? Os códigos são dados abaixo:

library(lme4)
read.table(textConnection("duration season  sites   effect
                          4d    mon s1  7305.91
                          4d    mon s2  856.297
                          4d    mon s3  649.93
                          4d    mon s1  10121.62
                          4d    mon s2  5137.85
                          4d    mon s3  3059.89
                          4d    mon s1  5384.3
                          4d    mon s2  5014.66
                          4d    mon s3  3378.15
                          4d    post    s1  6475.53
                          4d    post    s2  2923.15
                          4d    post    s3  554.05
                          4d    post    s1  7590.8
                          4d    post    s2  3888.01
                          4d    post    s3  600.07
                          4d    post    s1  6717.63
                          4d    post    s2  1542.93
                          4d    post    s3  1001.4
                          4d    pre s1  9290.84
                          4d    pre s2  2199.05
                          4d    pre s3  1149.99
                          4d    pre s1  5864.29
                          4d    pre s2  4847.92
                          4d    pre s3  4172.71
                          4d    pre s1  8419.88
                          4d    pre s2  685.18
                          4d    pre s3  4133.15
                          7d    mon s1  11129.86
                          7d    mon s2  1492.36
                          7d    mon s3  1375
                          7d    mon s1  10927.16
                          7d    mon s2  8131.14
                          7d    mon s3  9610.08
                          7d    mon s1  13732.55
                          7d    mon s2  13314.01
                          7d    mon s3  4075.65
                          7d    post    s1  11770.79
                          7d    post    s2  4254.88
                          7d    post    s3  753.2
                          7d    post    s1  11324.95
                          7d    post    s2  5133.76
                          7d    post    s3  2156.2
                          7d    post    s1  12103.76
                          7d    post    s2  3143.72
                          7d    post    s3  2603.23
                          7d    pre s1  13928.88
                          7d    pre s2  3208.28
                          7d    pre s3  8015.04
                          7d    pre s1  11851.47
                          7d    pre s2  6815.31
                          7d    pre s3  8478.77
                          7d    pre s1  13600.48
                          7d    pre s2  1219.46
                          7d    pre s3  6987.5
                          "),header=T)->dat1


m1 = lmer(effect ~ duration + (1+duration|sites) +(1+duration|season),
          data=dat1, REML=FALSE)

"Resolver" o problema que você enfrenta no sentido de não receber avisos sobre convergência com falha é bastante direto: você não usa o otimizador BOBYQA padrão, mas prefere usar a rotina de otimização do Nelder-Mead usada por padrão nas 1.0.xversões anteriores. Ou você instala o pacote optimxpara poder diretamente uma rotina L-BFGS-B ou nlminb(igual às lme4versões anteriores à versão 1). Por exemplo:

m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(optimizer ="Nelder_Mead")
library(optimx)
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='nlminb')))

tudo funciona bem (sem avisos). As perguntas interessantes são:

1. por que você recebeu esses avisos e
2. porque quando você usou REML = TRUEvocê não recebeu nenhum aviso.

Sucintamente, 1. você recebeu esses avisos porque definiu durationtanto um efeito fixo quanto uma inclinação aleatória para o fator sitese também season. O modelo esgotou efetivamente os graus de liberdade para estimar as correlações entre as inclinações e as interceptações que você definiu. Se você usou um modelo marginalmente mais simples como:

m1 = lmer(effect~duration+ (1+duration|sites) + (0+duration|season) + (1|season),
          data=dat1, REML = FALSE)

você não enfrentaria problemas de convergência. Este modelo estimaria efetivamente interceptações aleatórias não correlacionadas e inclinações aleatórias para cada uma season.

REML = FALSE$X$$y = X\beta + Z\gamma + \epsilon$$K$$KX = 0$$y$$Ky$$Z$$KZ$$Z$

Uma observação final é que não tenho certeza se faz sentido usar seasoncomo efeito aleatório para começar. Afinal, existem poucas temporadas, então você pode tratá-las como efeitos fixos.

A questão é estatística e não técnica. Na verdade, eu usei um modelo de efeito aleatório em vez de um modelo de efeito fixo. Acho que nenhum dos fatores deve ser tratado como fator aleatório, pois precisamos de pelo menos 5 ou 6 níveis ou réplicas para tratar um fator como efeito aleatório (consulte aqui Qual é o número mínimo recomendado de grupos para um fator de efeitos aleatórios? ).

O conjunto de dados acima contém apenas amostras / local / estação em triplicado, insuficiente para um modelo de efeito aleatório. No conjunto de dados, a duração de 4 e 7 dias pertence a duas experiências paralelas separadas, realizadas ao mesmo tempo. Portanto, citar os dados definidos por duração (4 e 7 dias) e realizar uma anova bidirecional para cada duração com a estação e os locais, pois os fatores seriam suficientes para modelar o efeito (variável de resposta) aqui. O modelo deve ser o seguinte:

lm(day_4_effect~sites*season, data=dat1)

lm(day_7_effect~sites*season, data=dat1)

Agradeço a Bodo Winter ( http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf ) e a @ usεr11852 que me ajudaram a resolver o problema.