Como estimar componentes de variação com o lmer para modelos com efeitos aleatórios e compará-los com os resultados do lme

14

Realizei um experimento em que criei famílias diferentes, provenientes de duas populações de fontes diferentes. Cada família recebeu um dos dois tratamentos. Após o experimento, medi várias características em cada indivíduo. Para testar um efeito do tratamento ou da fonte, bem como sua interação, usei um modelo linear de efeitos mistos com a família como fator aleatório, ou seja,

lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")

até agora tudo bem, agora eu tenho que calcular os componentes da variação relativa, ou seja, a porcentagem de variação que é explicada pelo tratamento ou pela fonte, bem como pela interação.

Sem um efeito aleatório, eu poderia facilmente usar as somas dos quadrados (SS) para calcular a variação explicada por cada fator. Mas para um modelo misto (com estimativa de ML), não há SS, portanto, pensei que poderia usar Tratamento e Origem como efeitos aleatórios também para estimar a variação, ou seja,

lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")

No entanto, em alguns casos, o lme não converge, portanto, usei o lmer do pacote lme4:

lmer(Trait~1+(Treatment*Source|Family),data=DATA)

Onde extraio as variações do modelo usando a função de resumo:

model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]

Eu recebo os mesmos valores da função VarCorr. Utilizo esses valores para calcular a porcentagem real de variação, tomando a soma como a variação total.

Onde estou lutando é com a interpretação dos resultados do modelo inicial de LME (com tratamento e fonte como efeitos fixos) e o modelo aleatório para estimar os componentes de variância (com tratamento e fonte como efeito aleatório). Na maioria dos casos, acho que a porcentagem de variação explicada por cada fator não corresponde à significância do efeito fixo.

Por exemplo, para a característica HD, o filme inicial sugere uma tendência para a interação e um significado para o tratamento. Usando um procedimento para trás, acho que o tratamento tem uma tendência próxima a significativa. No entanto, estimando os componentes da variação, acho que a Origem tem a variação mais alta, perfazendo 26,7% da variação total.

O lme:

anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
                                      numDF denDF  F-value p-value
(Intercept)                                1   426 0.044523  0.8330
as.factor(Treatment)                       1   426 5.935189  0.0153
as.factor(Source)                          1    11 0.042662  0.8401
as.factor(Treatment):as.factor(Source)     1   426 3.754112  0.0533

E o último:

summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
Linear mixed model fit by REML 
Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family) 
   Data: regrexpdat 
    AIC    BIC logLik deviance REMLdev
 -103.5 -54.43  63.75   -132.5  -127.5
Random effects:
 Groups   Name                                      Variance  Std.Dev. Corr                 
 Family   (Intercept)                               0.0113276 0.106431                      
          as.factor(Treatment)                      0.0063710 0.079819  0.405               
          as.factor(Source)                         0.0235294 0.153393 -0.134 -0.157        
          as.factor(Treatment)L:as.factor(Source)   0.0076353 0.087380 -0.578 -0.589 -0.585 
 Residual                                           0.0394610 0.198648                      
Number of obs: 441, groups: Family, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02740    0.03237  -0.846

Portanto, minha pergunta é: está correto o que estou fazendo? Ou devo usar outra maneira de estimar a quantidade de variação explicada por cada fator (ou seja, Tratamento, Origem e sua interação). Por exemplo, os tamanhos dos efeitos seriam um caminho mais apropriado?

KL_STKBK
fonte
O fator de tratamento tem 40x tantos graus de liberdade quanto o fator de origem (pseudo-replicação?). Sem dúvida, isso está diminuindo o valor P do tratamento.

Respostas:

1

Uma maneira comum de determinar a contribuição relativa de cada fator para um modelo é removê-lo e comparar a probabilidade relativa com algo como um teste qui-quadrado:

pchisq(logLik(model1) - logLik(model2), 1)

Como a maneira como as probabilidades são calculadas entre as funções pode ser um pouco diferente, normalmente as compararei apenas com o mesmo método.

Bill Denney
fonte
1
não deveria ser 1-pchisq(logLik(model1) - logLik(model2), 1)?
user81411