Qual é a diferença entre usar aov () e lme () na análise de um conjunto de dados longitudinal?

13

Alguém pode me dizer a diferença entre usar aov()e lme()analisar dados longitudinais e como interpretar os resultados desses dois métodos?

Abaixo, eu analisar o mesmo conjunto de dados usando aov()e lme()e tem 2 resultados diferentes. Com aov(), obtive um resultado significativo no tempo por interação com o tratamento, mas, ao ajustar um modelo linear misto, o tempo por interação com o tratamento é insignificante.

> UOP.kg.aov <- aov(UOP.kg~time*treat+Error(id), raw3.42)
> summary(UOP.kg.aov)

Error: id
          Df  Sum Sq Mean Sq F value Pr(>F)
treat      1   0.142  0.1421  0.0377 0.8471
Residuals 39 147.129  3.7725               

Error: Within
            Df  Sum Sq Mean Sq  F value  Pr(>F)    
time         1 194.087 194.087 534.3542 < 2e-16 ***
time:treat   1   2.077   2.077   5.7197 0.01792 *  
Residuals  162  58.841   0.363                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> UOP.kg.lme <- lme(UOP.kg~time*treat, random=list(id=pdDiag(~time)), 
                    na.action=na.omit, raw3.42)
> summary(UOP.kg.lme)
Linear mixed-effects model fit by REML
 Data: raw3.42 
       AIC      BIC    logLik
  225.7806 248.9037 -105.8903

Random effects:
 Formula: ~time | id
 Structure: Diagonal
        (Intercept)      time  Residual
StdDev:   0.6817425 0.5121545 0.1780466

Fixed effects: UOP.kg ~ time + treat + time:treat 
                 Value Std.Error  DF   t-value p-value
(Intercept)  0.5901420 0.1480515 162  3.986059  0.0001
time         0.8623864 0.1104533 162  7.807701  0.0000
treat       -0.2144487 0.2174843  39 -0.986042  0.3302
time:treat   0.1979578 0.1622534 162  1.220053  0.2242
 Correlation: 
           (Intr) time   treat 
time       -0.023              
treat      -0.681  0.016       
time:treat  0.016 -0.681 -0.023

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.198315285 -0.384858426  0.002705899  0.404637305  2.049705655 

Number of Observations: 205
Number of Groups: 41 
biostat_newbie
fonte

Respostas:

18

Com base na sua descrição, parece que você tem um modelo de medidas repetidas com um único fator de tratamento. Como não tenho acesso ao conjunto de dados ( raw3.42), usarei os dados da Orthodont donlme pacote para ilustrar o que está acontecendo aqui. A estrutura dos dados é a mesma (medidas repetidas para dois grupos diferentes - neste caso, homens e mulheres).

Se você executar o seguinte código:

library(nlme)
data(Orthodont)

res <- lme(distance ~ age*Sex, random = ~ 1 | Subject, data = Orthodont)
anova(res)

você obterá os seguintes resultados:

            numDF denDF  F-value p-value
(Intercept)     1    79 4123.156  <.0001
age             1    79  122.450  <.0001
Sex             1    25    9.292  0.0054
age:Sex         1    79    6.303  0.0141

Se você executar:

res <- aov(distance ~ age*Sex + Error(Subject), data = Orthodont)
summary(res)

você vai ter:

Error: Subject
          Df Sum Sq Mean Sq F value   Pr(>F)   
Sex        1 140.46 140.465  9.2921 0.005375 **
Residuals 25 377.91  15.117                    

Error: Within
          Df  Sum Sq Mean Sq  F value  Pr(>F)    
age        1 235.356 235.356 122.4502 < 2e-16 ***
age:Sex    1  12.114  12.114   6.3027 0.01410 *  
Residuals 79 151.842   1.922                     

Observe que os testes F são exatamente idênticos.

Para lme(), você usou list(id=pdDiag(~time)), que não apenas adiciona uma interceptação aleatória ao modelo, mas também uma inclinação aleatória. Além disso, usando pdDiag, você está configurando a correlação entre a interceptação aleatória e a inclinação para zero. Este é um modelo diferente do que você especificou via aov()e, portanto, você obtém resultados diferentes.

Wolfgang
fonte
Obrigado @Wolfgang; sua explicação ajuda muito. Uma pergunta de acompanhamento que eu tenho então é esta. Na verdade, estou analisando um modelo de medidas repetidas com um único fator de tratamento. Cada sujeito é designado aleatoriamente para o tratamento A ou B. Em seguida, são medidos em 0 minutos, 15 minutos, 30 minutos, 60 minutos, 120 minutos e 180 minutos. Pelo meu entendimento, o tempo deve ser um fator aleatório, pois são apenas amostras do tempo de 0 a 180 minutos. Então, devo fazer: lme (UOP.kg ~ time * Treat, random = ~ time | id, raw3.42)?
Biostat_newbie 15/02/11
Sim, mas eu pensaria assim: Você está essencialmente permitindo que a interceptação e a inclinação da linha de regressão (de UOP.kg no tempo) diferam (aleatoriamente) entre os sujeitos do mesmo grupo de tratamento. Isto é o que random = ~ time | id fará. O que o modelo irá lhe dizer é a quantidade estimada de variabilidade nas interceptações e nas inclinações. Além disso, o tempo: termo de interacção mimo indica se a média inclinação é diferente de A e B.
Wolfgang
Obrigado @Wolfgang! Posso usar Error(Subject/age), já que procurei em algum tutorial, dizendo que isso /agesignifica medidas repetidas nesse fator? É o mesmo que o seu Error(Subject)? Outra pergunta é: para dados desequilibrados, aove lmepode ter resultados diferentes, certo?
breezeintopl
1

Gostaria de acrescentar que você pode instalar o carpacote e usá- Anova()lo, em vez de anova()por aov()e paralm() objetos, a baunilha anova()usa uma soma seqüencial de quadrados, que fornece o resultado errado para tamanhos de amostra desiguais, enquanto lme()usa o tipo -I ou a soma dos quadrados do tipo III, dependendo do typeargumento, mas a soma dos quadrados do tipo III viola a marginalidade - ou seja, trata as interações não de maneira diferente dos efeitos principais.

A lista de ajuda R não tem nada de bom a dizer sobre as somas de quadrados do tipo I e III, e, no entanto, essas são as únicas opções! Vai saber.

Edit: Na verdade, parece que o tipo II não é válido se houver um termo de interação significativo, e parece que o melhor que alguém pode fazer é usar o tipo III quando houver interações. Fui informado por uma resposta a uma de minhas próprias perguntas que, por sua vez, me indicaram esse post .

f1r3br4nd
fonte
0

Parece-me que você tem várias medidas para cada identificação de cada vez. Você precisa agregá-las para o Aov, porque aumenta injustamente o poder nessa análise. Não estou dizendo que fazer o agregado fará com que os resultados sejam os mesmos, mas deve torná-los mais semelhantes.

dat.agg <- aggregate(UOP.kg ~ time + treat + id, raw3.42, mean)

Em seguida, execute seu modelo Aov como antes, substituindo os dados por dat.agg.

Além disso, acredito que anova (lme) é mais o que você deseja fazer para comparar os resultados. A direção e magnitude de um efeito não é a mesma que a razão entre a variação do modelo e o erro.

(BTW, se você fizer a análise lme nos dados agregados, o que você não deve fazer e verificar anova (lme), obterá quase os mesmos resultados que o aov)

John
fonte