O que fazer com a correlação de efeitos aleatórios igual a 1 ou -1?

Uma ocorrência não tão incomum ao lidar com modelos mistos máximos complexos (estimar todos os efeitos aleatórios possíveis para dados e modelo) é perfeita (+1 ou -1) ou correlação quase perfeita entre alguns efeitos aleatórios. Para o propósito da discussão, vamos observar o seguinte modelo e resumo do modelo

Model: Y ~ X*Cond + (X*Cond|subj)

# Y = logit variable  
# X = continuous variable  
# Condition = values A and B, dummy coded; the design is repeated 
#             so all participants go through both Conditions  
# subject = random effects for different subjects  

Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         CondB            0.54367  0.7373   -0.37  0.37      
         X:CondB          0.14812  0.3849    0.26 -0.26 -0.56
Number of obs: 39401, groups:  subject, 219

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.49686    0.06909   36.14  < 2e-16 ***
X                -1.03854    0.03812  -27.24  < 2e-16 ***
CondB            -0.19707    0.06382   -3.09  0.00202 ** 
X:CondB           0.22809    0.05356    4.26 2.06e-05 ***

A suposta razão por trás dessas correlações perfeitas é que criamos um modelo complexo demais para os dados que possuímos. O conselho comum que é dado nessas situações é (por exemplo, Matuschek et al., 2017; paper ) fixar os coeficientes superparameterizados em 0, porque esses modelos degenerados tendem a diminuir a potência. Se observarmos uma mudança acentuada nos efeitos fixos em um modelo reduzido, devemos aceitá-lo; se não houver alteração, não haverá problema em aceitar a original.

No entanto, vamos supor que não estamos interessados ​​apenas em efeitos fixos controlados por ER (efeitos aleatórios), mas também na estrutura de ER. Nesse caso, seria teoricamente correto supor que Intercepte a inclinação Xtenha correlação negativa diferente de zero. Várias perguntas a seguir:

1. O que fazer em tais situações? Devemos relatar a correlação perfeita e dizer que nossos dados não são "bons o suficiente" para estimar a correlação "real"? Ou devemos relatar o modelo de correlação 0? Ou deveríamos tentar definir outra correlação como 0, na esperança de que a "importante" não seja mais perfeita? Não acho que haja respostas 100% corretas aqui. Gostaria principalmente de ouvir suas opiniões.

2. Como escrever o código que fixa a correlação de 2 efeitos aleatórios específicos para 0, sem influenciar as correlações entre outros parâmetros?

Matrizes de covariância de efeito aleatório singulares

Obter uma estimativa de correlação de efeito aleatório de +1 ou -1 significa que o algoritmo de otimização atingiu "um limite": as correlações não podem ser maiores que +1 ou menores que -1. Mesmo se não houver erros ou avisos explícitos de convergência, isso potencialmente indica alguns problemas com a convergência, porque não esperamos que correlações verdadeiras fiquem no limite. Como você disse, isso geralmente significa que não há dados suficientes para estimar todos os parâmetros de maneira confiável. Matuschek et al. 2017 dizem que nessa situação o poder pode ser comprometido.

Outra maneira de atingir um limite é obter uma estimativa de variação de 0: Por que obtenho uma variação zero de um efeito aleatório no meu modelo misto, apesar de alguma variação nos dados?

Ambas as situações podem ser vistas como a obtenção de uma matriz de covariância degenerada de efeitos aleatórios (no seu exemplo, a matriz de covariância de saída é ); uma variância zero ou uma correlação perfeita significa que a matriz de covariância não é de classificação completa e [pelo menos] um de seus valores próprios é zero. Essa observação sugere imediatamente que existem outras maneiras mais complexas de obter uma matriz de covariância degenerada: pode-se ter uma matriz de covariância sem zeros ou correlações perfeitas, mas, no entanto, deficiente na classificação (singular). Bates et al. 2015 modelos mistos parcimoniosos4 × $4\times 4$$4\times 4$(pré-impressão não publicada) recomenda o uso da análise de componentes principais (PCA) para verificar se a matriz de covariância obtida é singular. Se for, eles sugerem tratar essa situação da mesma maneira que as situações singulares acima.

Então o que fazer?

Se não houver dados suficientes para estimar todos os parâmetros de um modelo de maneira confiável, devemos considerar a simplificação do modelo. Tomando o seu modelo de exemplo X*Cond + (X*Cond|subj), existem várias maneiras possíveis de simplificá-lo:

1. Remova um dos efeitos aleatórios, geralmente a correlação de maior ordem:

   X*Cond + (X+Cond|subj)

2. Livre-se de todos os parâmetros de correlação:

   X*Cond + (X*Cond||subj)

   Atualização: como o @Henrik observa, a ||sintaxe somente removerá correlações se todas as variáveis ​​à esquerda forem numéricas. Se variáveis ​​categóricas (como Cond) estiverem envolvidas, é melhor usar seu afexpacote conveniente (ou soluções manuais complicadas). Veja a resposta dele para mais detalhes.

3. Livre-se de alguns dos parâmetros de correlação dividindo o termo em vários, por exemplo:

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)

4. Restrinja a matriz de covariância de alguma maneira específica, por exemplo, definindo uma correlação específica (a que atingiu o limite) como zero, como você sugere. Não há uma maneira integrada de lme4conseguir isso. Veja a resposta do @ BenBolker no SO para uma demonstração de como conseguir isso através de alguns hackers inteligentes.

Ao contrário do que você disse, não acho que Matuschek et al. 2017 recomendam especificamente # 4. A essência de Matuschek et al. 2017 e Bates et al. Parece que 2015 começa com o modelo máximo a la Barr et al. 2013 e, em seguida, diminui a complexidade até que a matriz de covariância esteja no nível máximo. (Além disso, eles recomendariam frequentemente reduzir ainda mais a complexidade, a fim de aumentar a potência.) Atualização: Em contraste, Barr et al. recomende reduzir APENAS a complexidade se o modelo não convergir; eles estão dispostos a tolerar matrizes de covariância singulares. Veja a resposta de @ Henrik.

Se alguém concorda com Bates / Matuschek, acho que é bom tentar maneiras diferentes de diminuir a complexidade para encontrar a pessoa que faz o trabalho enquanto faz "o menor dano". Olhando para a minha lista acima, a matriz de covariância original possui 10 parâmetros; # 1 tem 6 parâmetros, # 2 tem 4 parâmetros, # 3 tem 7 parâmetros. Qual modelo se livrará das correlações perfeitas é impossível dizer sem ajustá-las.

Mas e se você estiver interessado neste parâmetro?

A discussão acima trata a matriz de covariância de efeito aleatório como um parâmetro incômodo. Você levanta uma questão interessante sobre o que fazer se estiver especificamente interessado em um parâmetro de correlação que precisa "desistir" para obter uma solução de classificação completa significativa.

Observe que fixar o parâmetro de correlação em zero não produzirá necessariamente BLUPs ( ranef) não correlacionados; na verdade, eles podem nem ser afetados tanto assim (veja a resposta de @ Placidia para uma demonstração ). Portanto, uma opção seria examinar as correlações dos BLUPs e reportar isso.

Outra opção, talvez menos atraente, seria usar o tratamento subjectcomo um efeito fixo Y~X*cond*subj, obter as estimativas para cada sujeito e calcular a correlação entre eles. Isso é equivalente a executar Y~X*condregressões separadas para cada sujeito separadamente e obter as estimativas de correlação deles.

Veja também a seção sobre modelos singulares nas perguntas frequentes sobre modelos mistos de Ben Bolker:

É muito comum que modelos mistos com excesso de ajuste resultem em ajustes singulares. Tecnicamente, singularidade significa que alguns dos parâmetros (decomposição de Cholesky da variância-covariância) correspondentes a elementos diagonais do fator de Cholesky são exatamente zero, que é a borda do espaço viável, ou equivalente que a matriz de variância-covariância tem zero autovalores (ou seja, é positivo semidefinido em vez de positivo positivo) ou (quase equivalente) que algumas das variações são estimadas como zero ou que algumas das correlações são estimadas como +/- 1.$\theta$

Concordo com tudo o que foi dito na resposta da ameba, que fornece um ótimo resumo da discussão atual sobre esse assunto. Tentarei acrescentar alguns pontos adicionais e, de outra forma, refiro-me à apostila do meu curso recente de modelos mistos, que também resume esses pontos.

Suprimir os parâmetros de correlação (opções 2 e 3 na resposta da ameba) via ||funciona apenas para covariáveis ​​numéricas lmere não para fatores. Isso é discutido em alguns detalhes com o código de Reinhold Kliegl .

No entanto, meu afexpacote fornece a funcionalidade para suprimir a correlação também entre fatores se argumento expand_re = TRUEna chamada para mixed()(consulte também a função lmer_alt()). Essencialmente, ele faz isso implementando a abordagem discutida por Reinhold Kliegl (isto é, transfomindo os fatores em covariáveis ​​numéricas e especificando a estrutura de efeitos aleatórios sobre eles).

Um exemplo simples:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

Para quem não sabe afex, a principal funcionalidade dos modelos mistos é fornecer valores de p para os efeitos fixos, por exemplo:

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

Dale Barr, de Barr et al. (2013) é mais cauteloso ao recomendar a redução da estrutura de efeitos aleatórios do que a apresentada na resposta da ameba. Em uma recente troca no Twitter, ele escreveu:

• "reduzir o modelo apresenta um risco desconhecido de anticonservatividade e deve ser feito com cautela, se for o caso". e
• "Minha principal preocupação é que as pessoas entendam os riscos associados à redução do modelo e que a minimização desse risco exija uma abordagem mais conservadora do que é comumente adotada (por exemplo, cada inclinação testada em 0,05)."

Portanto, recomenda-se cautela.

Como um dos revisores, também posso fornecer algumas dicas sobre por que nós, Bates et al. (2015) permaneceu inédito. Eu e os outros dois revisores (que assinamos, mas permaneceremos sem nome aqui), fizemos algumas críticas à abordagem PCA (parece sem princípios e não há evidências de que seja superior em termos de poder). Além disso, acredito que os três criticaram o fato de o artigo não se concentrar na questão de como especificar a estrutura de efeitos aleatórios, mas também tenta introduzir GAMMs. Assim, o artigo de Bates et al (2015) se transformou no Matuschek et al. (2017) artigo que aborda a questão da estrutura de efeitos aleatórios com simulações e o Baayen et al. (2017) artigo que apresenta GAMMs.

Minha revisão completa do Bates et al. rascunho pode ser encontrado aqui . IIRC, as outras análises tiveram pontos principais semelhantes.

Eu também tive esse problema ao usar a estimativa de máxima verossimilhança - apenas eu uso o algoritmo Goldstein IGLS conforme implementado pelo software MLwiN e não o LME4 em R. No entanto, em cada caso, o problema foi resolvido quando mudei para a estimativa MCMC usando o mesmo Programas. Eu até tive uma correlação acima de 3, que foi resolvida quando mudei de estimativa. Usando o IGLS, a correlação é calculada após a estimativa como a covariância dividida pelo produto da raiz quadrada do produto das variações associadas - e isso não leva em conta a incerteza em cada uma das estimativas constituintes.

O software IGLS não 'sabe' que a covariância implica uma correlação e apenas calcula estimativas de uma função de variação constante, linear, quadrática etc. Em contraste, a abordagem MCMC baseia-se na suposição de amostras de uma distribuição normal multivariada que corresponde a variações e covariâncias com boas propriedades e propagação total de erros, de modo que a incerteza na estimativa das covariâncias é levada em consideração na estimativa das variações. e vice versa.

O MLwin é a cadeia de estimativa do MCMC com estimativas do IGLS e a matriz de covariância de variância definida não-negativa pode precisar ser alterada alterando-se a covariância para zero no início antes de iniciar a amostragem.

Para um exemplo trabalhado, consulte

Desenvolvimento de modelos multiníveis para análise de contextualidade, heterogeneidade e mudança usando o MLwiN 3, volume 1 (atualizado em setembro de 2017); O volume 2 também está no RGate

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

Apêndice ao capítulo 10