Como escolher a estrutura de efeitos aleatórios e fixos em modelos lineares mistos?

Considere os seguintes dados de uma maneira bidirecional no design de assuntos:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

Eu gostaria de analisar isso usando modelos lineares mistos. Considerando todos os possíveis efeitos fixos e aleatórios, existem vários modelos possíveis:

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

1. Qual é a maneira recomendada de selecionar o modelo de melhor ajuste nesse contexto? Ao usar testes de razão de verossimilhança, qual é o procedimento recomendado? Gerando modelos para cima (do modelo nulo para o modelo mais complexo) ou para baixo (do modelo mais complexo para o modelo nulo)? Inclusão ou exclusão gradual? Ou é recomendável colocar todos os modelos em um teste de razão de verossimilhança e selecionar o modelo com o menor valor de p? Como comparar modelos que não estão aninhados?

2. É recomendável primeiro encontrar a estrutura apropriada de efeitos fixos e depois a estrutura apropriada de efeitos aleatórios ou o contrário (eu encontrei referências para ambas as opções ...)?

3. Qual é a maneira recomendada de relatar resultados? Relatar o valor-p do teste da razão de verossimilhança de log comparando o modelo misto completo (com o efeito em questão) ao modelo reduzido (sem o efeito em questão). Ou é melhor usar o teste da razão de verossimilhança de log para encontrar o melhor modelo de ajuste e, em seguida, usar o lmerTest para relatar valores de p dos efeitos no modelo de melhor ajuste?

Não tenho certeza se há realmente uma resposta canônica para isso, mas vou tentar.

Qual é a maneira recomendada de selecionar o modelo de melhor ajuste nesse contexto? Ao usar testes de razão de verossimilhança, qual é o procedimento recomendado? Gerando modelos para cima (do modelo nulo para o modelo mais complexo) ou para baixo (do modelo mais complexo para o modelo nulo)? Inclusão ou exclusão gradual? Ou é recomendável colocar todos os modelos em um teste de razão de verossimilhança e selecionar o modelo com o menor valor de p? Como comparar modelos que não estão aninhados?

Depende de quais são seus objetivos.

• Em geral, você deve ter muito , muito cuidado com a seleção de modelos (veja, por exemplo, esta resposta ou este post , ou apenas o Google "Harrell stepwise" ...).
• Se você estiver interessado em ter seus valores de p ser significativo (ou seja, você está fazendo testes de hipóteses de confirmação), você deve não fazer seleção de modelos. No entanto : não está tão claro para mim se os procedimentos de seleção de modelos são tão ruins se você estiver fazendo a seleção de modelos em partes não focais do modelo , por exemplo, fazendo a seleção de modelos nos efeitos aleatórios se o seu interesse principal for a inferência nos efeitos fixos.
• Não existe algo como "colocar todos os modelos em um teste de razão de verossimilhança" - o teste da razão de verossimilhança é um procedimento pareado. Se você quiser fazer a seleção do modelo (por exemplo) nos efeitos aleatórios, provavelmente recomendaria uma abordagem "de uma só vez" usando critérios de informação como neste exemplo - que pelo menos evite alguns dos problemas das abordagens passo a passo (mas não de seleção de modelos de maneira mais geral).
• Barr et al. O Journal of Memory and Language de 2013 "keep it maximal" (doi: 10.1016 / j.jml.2012.11.001) recomendaria o uso do modelo máximo (apenas).
• Shravan Vasishth discorda , argumentando que esses modelos serão insuficientes e, portanto, problemáticos, a menos que o conjunto de dados seja muito grande (e a relação sinal / ruído seja alta)
• Outra abordagem razoavelmente defensável é ajustar-se a um modelo grande, mas razoável e, se o ajuste for singular, remova os termos até que não seja mais
• Com algumas ressalvas (enumeradas nas Perguntas frequentes do GLMM ), é possível usar critérios de informação para comparar modelos não aninhados com diferentes efeitos aleatórios (embora Brian Ripley discorde: veja o final da p. 6 aqui )

É recomendável primeiro encontrar a estrutura apropriada de efeitos fixos e depois a estrutura apropriada de efeitos aleatórios ou o contrário (eu encontrei referências para ambas as opções ...)?

Eu acho que ninguém sabe. Consulte a resposta anterior sobre a seleção de modelos de maneira mais geral. Se você pudesse definir seus objetivos com clareza suficiente (o que poucas pessoas fazem), a pergunta pode ser respondida. Se você tiver referências para as duas opções, seria útil editar sua pergunta para incluí-las ... (Para o que vale a pena, este exemplo (já citado acima) usa critérios de informação para selecionar a parte de efeitos aleatórios, depois evita a seleção no parte de efeito fixo do modelo.

Qual é a maneira recomendada de relatar resultados? Relatar o valor-p do teste da razão de verossimilhança de log comparando o modelo misto completo (com o efeito em questão) ao modelo reduzido (sem o efeito em questão). Ou é melhor usar o teste da razão de verossimilhança de log para encontrar o melhor modelo de ajuste e, em seguida, usar o lmerTest para relatar valores de p dos efeitos no modelo de melhor ajuste?

Esta é (infelizmente) outra questão difícil. Se você relatar os efeitos marginais como relatado por lmerTest, você tem que se preocupar com a marginalidade (por exemplo, se as estimativas dos principais efeitos Ae Bsão significativos quando há uma A-by- Binteração no modelo); essa é uma enorme lata de worms, mas é um pouco atenuada se você usar contrasts="sum"como recomendado por afex::mixed(). Projetos equilibrados também ajudam um pouco. Se você realmente deseja cobrir todas essas rachaduras, acho que recomendo afex::mixed, o que fornece uma saída semelhante a lmerTest, mas tenta lidar com esses problemas.

Atualização de maio de 2017 : Acontece que muito do que escrevi aqui é meio errado . Algumas atualizações são feitas ao longo do post.

Eu concordo muito com o que já foi dito por Ben Bolker (obrigado pela mensagem afex::mixed()), mas deixe-me acrescentar algumas idéias mais gerais e específicas sobre esse assunto.

Concentre-se em efeitos fixos versus efeitos aleatórios e em como reportar resultados

Para o tipo de pesquisa experimental representada no conjunto de dados de exemplo de Jonathan Baron, você usa a pergunta importante: geralmente se um fator manipulado tem ou não um efeito geral. Por exemplo, encontramos um efeito principal geral ou interação de Task? Um ponto importante é que nesses conjuntos de dados geralmente todos os fatores estão sob controle experimental completo e são aleatoriamente designados. Consequentemente, o foco de interesse geralmente é nos efeitos fixos.
Por outro lado, os componentes de efeitos aleatórios podem ser vistos como parâmetros "incômodos" que capturam variação sistemática (isto é, diferenças interindividuais no tamanho do efeito) que não são necessariamente importantes para a questão principal. Desse ponto de vista, sugere-se o uso da estrutura máxima de efeitos aleatórios, como preconizado por Barr et al. segue um pouco naturalmente. É fácil imaginar que uma manipulação experimental não afeta todos os indivíduos da mesma maneira e queremos controlar isso. Por outro lado, o número de fatores ou níveis geralmente não é muito grande, de modo que o risco de sobreajuste parece comparativamente pequeno.

Consequentemente, eu seguiria a sugestão de Barr et al. e especifique uma estrutura máxima de efeitos aleatórios e relate os testes dos efeitos fixos como meus principais resultados. Para testar os efeitos fixos, sugiro também usarafex::mixed() pois ele relata testes de efeitos ou fatores (em vez de testes de parâmetros) e calcula esses testes de uma maneira um tanto sensível (por exemplo, usa a mesma estrutura de efeitos aleatórios para todos os modelos em que um o efeito único é removido, usa contrastes soma-zero, oferece métodos diferentes para calcular valores- p , ...).

E os dados de exemplo

O problema com os dados de exemplo que você forneceu é que, para esse conjunto de dados, a estrutura máxima de efeitos aleatórios leva a um modelo super saturado, pois há apenas um ponto de dados por célula do design:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

Consequentemente, lmer engasga com a estrutura máxima de efeitos aleatórios:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

Infelizmente, até onde sei, não existe uma maneira acordada de lidar com esse problema. Mas deixe-me esboçar e discutir algumas:

1. Uma primeira solução poderia ser remover a maior inclinação aleatória e testar os efeitos para este modelo:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67

   No entanto, essa solução é um pouco ad-hoc e não excessivamente motivada.

   Atualização de maio de 2017: esta é a abordagem que estou adotando atualmente. Veja esta postagem no blog e o rascunho do capítulo em que sou coautor. , seção "Estruturas de efeitos aleatórios para projetos tradicionais de ANOVA".

2. Uma solução alternativa (e que poderia ser vista como defendida pela discussão de Barr et al.) Poderia ser sempre remover as inclinações aleatórias para o menor efeito. Porém, isso tem dois problemas: (1) Que estrutura de efeitos aleatórios usamos para descobrir qual é o menor efeito e (2) R reluta em remover um efeito de ordem inferior, como um efeito principal, se efeitos de ordem superior, como um efeito interação deste efeito está presente (consulte aqui ). Como conseqüência, seria necessário configurar essa estrutura de efeitos aleatórios manualmente e passar a matriz do modelo assim construído para a chamada mais próxima.

3. Uma terceira solução poderia ser usar uma parametrização alternativa da parte de efeitos aleatórios, a saber, uma que corresponda ao modelo RM-ANOVA para esses dados. Infelizmente (?), lmerNão permite "variações negativas", portanto essa parametrização não corresponde exatamente ao RM-ANOVA para todos os conjuntos de dados ; veja a discussão aqui e em outros lugares (por exemplo, aqui e aqui ). O "lmer-ANOVA" para esses dados seria:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76

Considerando todos esses problemas, eu simplesmente não usaria lmerpara ajustar conjuntos de dados para os quais existe apenas um ponto de dados por célula do design, a menos que esteja disponível uma solução mais acordada para o problema da estrutura de efeitos aleatórios máxima.

4. Em vez disso, eu também poderia usar a ANOVA clássica. Usar um dos wrappers para car::Anova()nos afexresultados seria:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

   Observe que afexagora também permite retornar o modelo equipado com o aovqual pode ser passado para lsmeanstestes post-hoc (mas, para testes de efeitos, os relatados por car::Anovaainda são mais razoáveis):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95