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?
Respostas:
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:
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.
fonte
1-pchisq(logLik(model1) - logLik(model2), 1)
?