Fernando,
Por que você usa na função lmer o argumento REML=FALSE se na
função lme o default é exatamente a estimação REML
?
Você não está comparando coisas diferentes?
?lme
method: a character string. If ‘"REML"’ the model is fit by
maximizing the restricted log-likelihood. If ‘"ML"’ the
log-likelihood is maximized. Defaults to ‘"REML"’.
E os NA's produzidos não seriam devido você estar definindo tempo
como efeito fixo e aleatório?
modelo0 <-
lme(GAS~HIDRAT*DILU*TEMP,random=~1|TEMP,weights=varIdent(form=~1|TEMP),data=dados)
Ao gerar um gráfico com seus dados notei uma diferença na inclinação
da variável GAS ao longo do tempo:
library(lattice)
xyplot(GAS~TEMP|HIDRAT+DILU, groups = BLOC ,data=dados,type='b')
Como sugestão eu ajustaria um modelo considerando o intercepto e a
inclinação como efeito aleatório e sem o efeito fixo da variável
DILU:
# Intercepto
modelo0 <- lme(GAS~HIDRAT*as.numeric(TEMP),random=~1|BLOC,
weights=varIdent(form=~1|TEMP), data=dados)
# Inclinação
modelo0.1 <-
lme(GAS~HIDRAT*as.numeric(TEMP),random=~as.numeric(TEMP)-1|BLOC,
weights=varIdent(form=~1|TEMP), data=dados)
# Intercepto e Inclinaçãp
modelo0.2 <-
lme(GAS~HIDRAT*as.numeric(TEMP),random=~as.numeric(TEMP)|BLOC,
weights=varIdent(form=~1|TEMP), data=dados)
anova(modelo0,modelo0.1,modelo0.2)
Model df AIC BIC
logLik Test L.Ratio p-value
modelo0 1 13 814.6697 854.2343
-394.3349
modelo0.1 2 13 782.5680 822.1325
-378.2840
modelo0.2 3 15 773.8100 819.4613 -371.9050 2
vs 3 12.75806 0.0017
--
Atenciosamente
Felipe E. Barletta Mendes
Estatístico(UFPR) - Conre3 9766-A
Mestrando em Bioestatística(UEM)
+55 (41)-92077191
+55 (41)-33287216