Caros colegas. Estou com dados de um experimento realizado em esquema de parcela subdividida em delineamento inteiramente casualizado. Há dois fatores experimentais A ( 2 níveis) e B (dois níveis) e a resposta Y foi medida repetidamente em tempos diferentes ( 23 tempos), sendo o tempo a subparcela. Eu pretendo analisá-lo através de Modelo Misto, utilizando a função lme{nlme}. Ao analisá-los utilizando a função lme{ nlme}, a função converge, porém recebo mensagens de warnings , avisando que NA's foram produzidos modelo0 <- lme(GAS~HIDRAT*DILU*TEMP,random=~1|TEMP,weights=varIdent(form=~1|TEMP),data=dados) ℝ> anova(modelo0) numDF denDF F-value p-value (Intercept) 1 368 5594.089 <.0001 HIDRAT 1 368 246.749 <.0001 DILU 1 368 2.667 0.1033 TEMP 22 0 206.379 *NaN* HIDRAT:DILU 1 368 6.346 0.0122 HIDRAT:TEMP 22 368 4.199 <.0001 DILU:TEMP 22 368 0.374 0.9961 HIDRAT:DILU:TEMP 22 368 0.038 1.0000 Warning message: In pf(Fval[i], nDF[i], dDF[i]) : NaNs produzidos O mesmo *não* ocorre quando eu utilizo a função lmer{lmerTest} modelo <-lmer(GAS~HIDRAT*DILU*TEMP + (1|TEMP),REML=FALSE,data=dados) ℝ> anova(modelo) Analysis of Variance Table of type III with Satterthwaite approximation for degrees of freedom Sum Sq Mean Sq NumDF DenDF F.value Pr(>F) HIDRAT 2163 2163.24 1 150285 330.58 < 2.2e-16 *** DILU 52 52.32 1 150285 8.00 0.004689 ** TEMP 32798 1490.80 22 150285 227.82 < 2.2e-16 *** HIDRAT:DILU 21 20.66 1 150285 3.16 0.075610 . HIDRAT:TEMP 215 9.77 22 150285 1.49 0.064140 . DILU:TEMP 78 3.53 22 150285 0.54 0.960247 HIDRAT:DILU:TEMP 3 0.13 22 150285 0.02 1.000000 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 A questão é que meus dados são heterocedasticos e gostaria de utilizar a função lme{nlme} devido as facilidades implementadas para lidar com dados assim. Procurei na literatura e estas facilidades ainda não se encontram implementadas na função lmer. Segue CMR dados<-structure(list(HIDRAT = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("1", "12"), class = "factor"), DILU = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), .Label = c("1", "3"), class = "factor"), TEMP = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L), .Label = c("1", "2", "3", "4", "5", "6", "8", "10"), class = "factor"), BLOC = structure(c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("1", "2", "3", "4", "5"), class = "factor"), GAS = c(4.61, 1.91, 0.16, 2.33, 2.64, 0.51, 2.02, 0.48, 1.44, 2.74, 7.24, 5.03, 0.8, 3.61, 7.02, 1.15, 4.7, 1.82, 2.75, 6.8, 8.13, 7.54, 1.49, 3.77, 9.61, 1.48, 7.2, 2.94, 3.59, 9.18, 8.64, 9.84, 4.03, 4.35, 11.67, 1.75, 9.23, 5.81, 3.96, 10.68, 9.27, 13.17, 6.77, 5.72, 13.18, 2.12, 11.49, 7.9, 4.92, 11.84, 9.99, 15.55, 9.31, 7.53, 15.36, 2.42, 13.21, 9.67, 6.44, 13.17, 13.48, 18.07, 14.69, 10.65, 19.52, 4.92, 18.46, 15.12, 9.48, 16.46, 16.96, 23.01, 20.95, 13.55, 24.01, 8.67, 23.56, 19.49, 12.43, 5.7, 5.11, 3.4, 4.73, 4.72, 4.05, 6.88, 4.78, 5.05, 4.5, 9.25, 9.92, 6.04, 7.96, 8.62, 7.12, 11.49, 8.4, 9.63, 8.89, 10.88, 13.5, 8.78, 8.6, 10.38, 8.66, 15.73, 11.27, 10.73, 11.02, 12.32, 16.06, 13.49, 9.42, 11.79, 9.82, 18.54, 16.25, 11.85, 12.44, 14.08, 19.72, 17.72, 10.98, 13.25, 11.11, 21.55, 19.76, 13.42, 13.37, 15.35, 22.71, 20.94, 12.7, 15.21, 12.21, 24.21, 22.36, 15.45, 14.53, 21.57, 26.13, 28.41, 16.8, 19.97, 17.38, 29.57, 29.18, 19.7, 17.78, 28.34, 31.73, 36.12, 22.04, 24.64, 24.89, 35.99, 35.4, 24.09, 21.46)), .Names = c("HIDRAT", "DILU", "TEMP", "BLOC", "GAS"), row.names = c(NA, -159L), class = "data.frame") install.packages(c("nlme","lmerTest")) #instale caso não tenha instalado em seu sistema library(nlme) library (lmerTest) modelo0 <- lme(GAS~HIDRAT*DILU*TEMP,random=~1|TEMP,weights=varIdent(form=~1|TEMP),data=dados) modelo <-lmer(GAS~HIDRAT*DILU*TEMP + (1|TEMP),REML=FALSE,data=dados) ==============================
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
Olá Felipe, Obrigado por responder a minha questão. Eu a ajustei REML=FALSE na função lmer{lmerTest}, não para a função lme{nlme}. Na função lmer isto foi recomendado pela função nas mensagens "warnings" e ao fazer o modelo converge perfeitamente. Mas não tenho costume em utilizar esta função e não estou seguro se este procedimento é correto. Além disto a fumção lmer não flexibiliza para trabalhar com dados heterocedastico por isso quero utilizar a lme{nlme} o qual estou acostumado. Sobre as suas soluções eu tenho interesse em avaliar se há efeito de interação entre Hidratação (HIDRAT) e Diluição (DILU) ao longo do tempo, por isso acho não ser possível excluir DILU do modelo. O que penso é que a parcela neste experimento, seria as unidades experimentais que receberam as combinações dos níveis dos fatores HIDRAT e DILU. A subparcela seria o tempo (TEMP). Eu tenho dúvidas se utilizo mesmo a abordagem do modelo misto. visto que neste experimento houve somente 1 aleatorização , já que o TEMPO não pode ser aleatorizado. Talvez algo assim seja mais apropridado. lme(GAS~HIDRAT*DILU,random=~1+as.numeric(TEMP)|BLOC,weights=varIdent(form=~1|T EMP),data=dados) Att ======================================================================= Fernando Souza Zootecnista, DSc. Produção e Alimentação Animal celular: (31)99796-8781 (Vivo) / (31)97358-4685 (Tim) e-mail:nandodesouza@gmail.com Lattes: http://lattes.cnpq.br/6519538815038307 blog: https://producaoanimalcomr.wordpress.com/ ======================================================================== On Mai 16 2016, at 4:00 pm, Felipe <felipe.e.barletta@gmail.com> wrote:
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
participantes (2)
-
Felipe -
Fernando Antonio de souza