Usando lmer para modelo de efeito misto linear de medidas repetidas

EDIÇÃO 2: Originalmente, pensei que precisava executar uma ANOVA de dois fatores com medidas repetidas em um fator, mas agora acho que um modelo linear de efeito misto funcionará melhor para meus dados. Acho que quase sei o que precisa acontecer, mas ainda estou confuso com alguns pontos.

Os experimentos que preciso analisar são assim:

• Os indivíduos foram designados para um dos vários grupos de tratamento
• As medidas de cada sujeito foram realizadas em vários dias
• Tão:
  • O sujeito está aninhado dentro do tratamento
  • O tratamento é cruzado com o dia

(cada sujeito é atribuído a apenas um tratamento e são tomadas medidas em cada sujeito em cada dia)

Meu conjunto de dados contém as seguintes informações:

• Assunto = fator de bloqueio (fator aleatório)
• Dia = dentro do sujeito ou fator de medidas repetidas (fator fixo)
• Tratamento = entre o fator sujeito (fator fixo)
• Obs = variável medida (dependente)

ATUALIZAÇÃO OK, então fui falar com um estatístico, mas ele é um usuário do SAS. Ele acha que o modelo deve ser:

Tratamento + Dia + Assunto (Tratamento) + Dia * Assunto (Tratamento)

Obviamente, sua notação é diferente da sintaxe R, mas esse modelo deve explicar:

• Tratamento (fixo)
• Dia (fixo)
• a interação Tratamento * Dia
• Assunto aninhado no tratamento (aleatório)
• Dia cruzado com "Assunto dentro do tratamento" (aleatório)

Então, essa é a sintaxe correta para usar?

m4 <- lmer(Obs~Treatment*Day + (1+Treatment/Subject) + (1+Day*Treatment/Subject), mydata)

Estou particularmente preocupado se o dia cruzado com a parte "Assunto dentro do tratamento" está correto. Alguém está familiarizado com o SAS ou confiante de que entende o que está acontecendo em seu modelo, capaz de comentar se minha triste tentativa de sintaxe R corresponde?

Aqui estão minhas tentativas anteriores de criar um modelo e escrever sintaxe (discutidas nas respostas e comentários):

m1 <- lmer(Obs ~ Treatment * Day + (1 | Subject), mydata)

Como faço para lidar com o fato de o sujeito estar aninhado dentro do tratamento? Qual é a m1diferença de:

m2 <- lmer(Obs ~ Treatment * Day + (Treatment|Subject), mydata)
m3 <- lmer(Obs ~ Treatment * Day + (Treatment:Subject), mydata)

e são m2e m3equivalentes (e se não, por que)?

Além disso, preciso usar o nlme, em vez do lme4, se desejar especificar a estrutura de correlação (como correlation = corAR1)? De acordo com medidas repetidas , para uma análise de medidas repetidas com medidas repetidas em um fator, a estrutura de covariância (a natureza das correlações entre medidas do mesmo sujeito) é importante.

Quando estava tentando fazer uma ANOVA de medidas repetidas, decidi usar um SS Tipo II; isso ainda é relevante? Em caso afirmativo, como especifico isso?

Aqui está um exemplo de como são os dados:

mydata <- data.frame(
  Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
               34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
               19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
               40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
               29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
  Day       = c(rep(c("Day1", "Day3", "Day6"), each=28)), 
  Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
                      "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
  Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293, 
                8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689, 
                5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737, 
                7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293, 
                6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709, 
                7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196, 
                6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612, 
                7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383, 
                6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068, 
                7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754, 
                6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802, 
                7.764172, 7.789844, 7.616437, NA, NA, NA, NA))

Eu acho que sua abordagem está correta. Modelo m1especifica uma interceptação separada para cada sujeito. O modelo m2adiciona uma inclinação separada para cada assunto. Sua inclinação é de vários dias, pois os participantes participam apenas de um grupo de tratamento. Se você escrever o modelo da m2seguinte forma, é mais óbvio que você modela uma interceptação e uma inclinação separadas para cada sujeito

m2 <- lmer(Obs ~ Treatment * Day + (1+Day|Subject), mydata)

Isso é equivalente a:

m2 <- lmer(Obs ~ Treatment + Day + Treatment:Day + (1+Day|Subject), mydata)

Ou seja, os principais efeitos do tratamento, dia e interação entre os dois.

Eu acho que você não precisa se preocupar em aninhar, desde que não repita os IDs dos sujeitos nos grupos de tratamento. Qual modelo está correto, realmente depende da sua pergunta de pesquisa. Há razões para acreditar que as inclinações dos sujeitos variam além do efeito do tratamento? Você pode executar os dois modelos e compará-los com anova(m1,m2)para ver se os dados suportam qualquer um.

Não tenho certeza do que você deseja expressar com o modelo m3? A sintaxe de aninhamento usa a /, por exemplo (1|group/subgroup).

Não acho que você precise se preocupar com a autocorrelação com um número tão pequeno de pontos no tempo.

Não me sinto à vontade para comentar sobre o problema de erros autocorrelacionados (nem sobre as diferentes implementações no lme4 vs. nlme), mas posso falar com o resto.

Seu modelo m1é um modelo de interceptação aleatória, no qual você incluiu a interação em nível cruzado entre Tratamento e Dia (o efeito de Dia pode variar entre os grupos Tratamento). Para permitir que a mudança ao longo do tempo seja diferente entre os participantes (ou seja, para modelar explicitamente as diferenças individuais na mudança ao longo do tempo), você também precisa permitir que o efeito do Dia seja aleatório . Para fazer isso, você deve especificar:

m2 <- lmer(Obs ~ Day + Treatment + Day:Treatment + (Day | Subject), mydata)

Neste modelo:

• A interceptação se a pontuação prevista para a categoria de referência de tratamento no dia = 0
• O coeficiente para dia é a alteração prevista ao longo do tempo para cada aumento de 1 unidade em dias para a categoria de referência de tratamento
• Os coeficientes para os dois códigos fictícios para os grupos de tratamento (criados automaticamente por R) são a diferença prevista entre cada grupo de tratamento restante e a categoria de referência no dia = 0
• Os coeficientes para os dois termos de interação são a diferença no efeito do tempo (dia) nas pontuações previstas entre a categoria de referência e os demais grupos de tratamento.

As interceptações e o efeito do dia na pontuação são aleatórios (cada sujeito pode ter uma pontuação prevista diferente no dia = 0 e uma mudança linear diferente ao longo do tempo). A covariância entre interceptações e declives também está sendo modelada (eles podem covary).

Como você pode ver, a interpretação dos coeficientes para as duas variáveis ​​dummy é condicional no dia = 0. Eles informarão se a pontuação prevista no dia = 0 para a categoria de referência é significativamente diferente dos dois grupos de tratamento restantes. Portanto, onde você decide centralizar sua variável Day é importante. Se você centralizar no dia 1, os coeficientes informarão se a pontuação prevista para a categoria de referência no dia 1 é significativamente diferente da pontuação prevista dos dois grupos restantes. Dessa forma, você pode ver se existem diferenças pré-existentes entre os grupos . Se você centralizar no dia 3, os coeficientes informarão se a pontuação prevista para a categoria de referência no dia 3é significativamente diferente da pontuação prevista dos dois grupos restantes. Dessa forma, você pode ver se há diferenças entre os grupos no final da intervenção .

Por fim, observe que os sujeitos não estão aninhados no tratamento. Seus três tratamentos não são níveis aleatórios de uma população de níveis aos quais você deseja generalizar seus resultados - em vez disso, como você mencionou, seus níveis são fixos e você deseja generalizar seus resultados apenas para esses níveis. (Sem mencionar, você não deve usar modelagem multinível se tiver apenas três unidades de nível superior; consulte Maas & Hox, 2005.) Em vez disso, o tratamento é um preditor de nível 2, ou seja, um preditor que assume um valor único por dias (unidades de nível 1) para cada sujeito. Portanto, ele é meramente incluído como um preditor no seu modelo.

Referência:
Maas, CJM, & Hox, JJ (2005). Tamanhos de amostra suficientes para modelagem multinível. Metodologia: European Journal of Research Methods for the Behavioral and Social Sciences , 1 , 86-92.