Funções lme e gls - modelos mistos
Olá pessoal, boa tarde!! Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R). Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais. 1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R? 2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente. Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti. ### Segue o código (com dados fictícios) para exemplificar o meu problema: local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp) require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE) ### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados) ### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
Por que não considera os blocos como efeito fixo ? Você vai assumir que a estrutura do efeito espacial é sempre a mesma dentro dos blocos ? Em 17 de dezembro de 2013 16:00, Cássio Dessotti < cassiodessotti@yahoo.com.br> escreveu:
Olá pessoal, boa tarde!!
Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R).
Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais.
1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R?
2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente.
Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti.
### Segue o código (com dados fictícios) para exemplificar o meu problema:
local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp)
require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE)
### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)),
correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
Penso que 'Blocos' não pode ser considerado como efeito fixo. Blocos são agrupamentos particulares das repetições. Portanto são repetições. Para que os testes estatísticos sejam válidos as repetições tem que ser aleatórias, caso contrário o experimento não tem valor, posto que as repetições são tendenciosas e as conclusões não podem ser inferidas para a população. Além disto o resíduo (aleatório) contém a interação blocos*tratamentos, que, portanto, deve ser aleatória. Note que no delineamento em que os tratamentos são aleatorizados nos blocos o gl do resíduo é (t-1)*(b-1). O resíduo é sempre aleatório porque o bloco é aleatório. Luiz Roberto Martins Pinto Prof. Pleno/DCET/UESC Laboratório de Estatística Computacional Universidade Estadual de Santa Cruz Ilhéus-Bahia luizroberto.uesc@gmail.com skype: lrmpinto http://lattes.cnpq.br/2732314327604831 Em 18 de dezembro de 2013 10:37, Wagner Bonat <wbonat@gmail.com> escreveu:
Por que não considera os blocos como efeito fixo ?
Você vai assumir que a estrutura do efeito espacial é sempre a mesma dentro dos blocos ?
Em 17 de dezembro de 2013 16:00, Cássio Dessotti < cassiodessotti@yahoo.com.br> escreveu:
Olá pessoal, boa tarde!!
Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R).
Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais.
1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R?
2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente.
Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti.
### Segue o código (com dados fictícios) para exemplificar o meu problema:
local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp)
require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE)
### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)),
correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Luiz Roberto e Wagner Acredito que a utilização de blocos como sendo de efeito aleatório geralmente é mais adequada, porém também estou fazendo o modelo com blocos fixo a nível de comparação com os demais (e neste não tive problemas no R - função gls), o problema está quando quero fazer o modelo (M3) com efeito de blocos aleatório, e função exponencial para a matriz de covariâncias dos resíduos. Wagner a estrutura de efeito espacial é a mesma dentro de cada local. Abraços. Cássio Dessotti. Em Quarta-feira, 18 de Dezembro de 2013 13:07, Luiz Roberto Martins Pinto <luizroberto.uesc@gmail.com> escreveu: Penso que 'Blocos' não pode ser considerado como efeito fixo. Blocos são agrupamentos particulares das repetições. Portanto são repetições. Para que os testes estatísticos sejam válidos as repetições tem que ser aleatórias, caso contrário o experimento não tem valor, posto que as repetições são tendenciosas e as conclusões não podem ser inferidas para a população. Além disto o resíduo (aleatório) contém a interação blocos*tratamentos, que, portanto, deve ser aleatória. Note que no delineamento em que os tratamentos são aleatorizados nos blocos o gl do resíduo é (t-1)*(b-1). O resíduo é sempre aleatório porque o bloco é aleatório. Luiz Roberto Martins Pinto Prof. Pleno/DCET/UESC Laboratório de Estatística Computacional Universidade Estadual de Santa Cruz Ilhéus-Bahia luizroberto.uesc@gmail.com skype: lrmpinto http://lattes.cnpq.br/2732314327604831 Em 18 de dezembro de 2013 10:37, Wagner Bonat <wbonat@gmail.com> escreveu: Por que não considera os blocos como efeito fixo ?
Você vai assumir que a estrutura do efeito espacial é sempre a mesma dentro dos blocos ?
Em 17 de dezembro de 2013 16:00, Cássio Dessotti <cassiodessotti@yahoo.com.br> escreveu:
Olá pessoal, boa tarde!!
Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R).
Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais.
1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R?
2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente.
Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti.
### Segue o código (com dados fictícios) para exemplificar o meu problema: local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp)
require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE)
### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Segue CMR com comentários. require(lattice) require(nlme) ## Segue o código (com dados fictícios) para exemplificar o meu problema: local <- as.factor(c(rep(1,6), rep(2,6))); blo <- as.factor(rep(c(rep(1,3), rep(2,3)),2)) ## Para que criar como fator se vai usar como numérico? ## lat <- as.factor(rep(1:3,4)); ## long <- as.factor(rep(c(rep(1,3), rep(2,3)),2)) lat <- rep(1:3,4); long <- rep(c(rep(1,3), rep(2,3)),2); xyplot(lat~long|local, groups=blo) ## Dentro de um bloco não se tem variação no valor de longitude. Isso ## pode implicar em problemas de identificabilidade dependendo do modelo ## considerado para os dados. O Elias já comentou isso. ## blo_local <- as.factor(c(rep("1_1",3), rep("2_1",3), rep("1_2",3), rep("2_2",3))) blo_local <- interaction(blo, local) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local, trat, blo, lat, long, blo_local, resp) str(dados) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp~1+trat+local+local:trat, random=list(blo_local=pdIdent(~1)), method="REML", na.action=na.omit, data=dados, keep.data=FALSE) summary(M1) ### M2 - Exp H (sem efeito de blocos - função exponencial para R - ### heterogeneidade de variâncias nos locais) M2 <- gls(resp~1+trat+local+local:trat, weight=varComb(varIdent(form=~1|local)), correlation=corExp(form=~lat+long|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados) summary(M2) ### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) ## M3 <- lme(resp~1+trat+local+local:trat, random=list(blo_local=pdIdent(~1)), ## weight=varComb(varIdent(form=~1|local)), ## correlation=corExp(form=~lat+long|local, ## metric="euclidean", nugget=FALSE), ## method="REML", na.action=na.omit, data=dados) ## ## Erro em lme.formula(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), : ## incompatible formulas for groups in 'random' and 'correlation' ## O modelo acima apresenta problema de especificação. Teóricamente é um ## modelo válido mas pelo visto não declarável na lme(). Ela requer que ## a correlação seja dentro dos níveis do fator de efeito aleatório, no ## caso, os blocos dentro de locais, e por isso que não aceita o tempo ## '|local' par o argumento correlation. Abaixo segue um modelo válido ## que assume que existe correlção entre observações dentro do mesmo ## bloco, e portanto, observações e diferentes blocos são ## independentes. Se entre essas parcelas existe uma distância grande, ## até que é uma suposição plausível. M4 <- lme(resp~1+trat+local+local:trat, random=list(blo_local=pdIdent(~1)), weight=varComb(varIdent(form=~1|local)), correlation=corExp(form=~lat+long, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados) summary(M4) Todos os modelo considerados são válidos, muito embora eu considere o maior modelo muito grande, pelo menos para os dados usados como exemplo (especificação de efeitos fixos, efeitos aleatórios, correlação entre observações, variância separada por local). É possível simular dados para um modelo com todos esses termos mas nem sempre é possível estimar a partir de dados, principalmente com variâncias de efeito aleatório indo para zero. Modelos grandes assim começam apresentar fraca identificabilidade. Não foi dada uma descrição minuciosa do experimento e objetivo da análise mas eu acredito que o modelo mais adequado (embora grande) seja aquele no qual os blocos dentro de locais sejam aleatórios, exista correlação entre parcelas dentro um mesmo local (inclusive de blocos diferentes) e variância separada por local. Na minha experiência, como na pratica existe uma relação média variância, locais mais produtivos tem maior variabilidade também e declarar variância heterogênea pode fazer a diferença. A questão é se esse modelo é declarável. Quanto aos blocos terem ou não efeito efeito aleatório eu não vejo problemas em assumir efeito fixo não. A análise continua válida mas muda um pouco a interpretação (inferência ampla ou restrita, ver Schabenberger & Pearce (2002?)). Blocos com efeito aleatório tem um apelo maior, principalmente no caso desbalanceado. A verossimilhança é resultado de uma integral sobre o domínio dos efeitos aleatórios, o que pode ser visto como uma espécie de verossimilhança média. Semanas atrás o Éder Borges (dessa lista) apresentou para nós do LEG e o Havard Rue o ajuste de um modelo, via INLA, muito interessante. Se eu não estiver desmemoriado, era um experimento de melhoramento genético, um local, muitos blocos (linha coluna) de efeito aleatório, material genético também aleatório, correlação entre parcelas (que são bem próximas), havia muito desbalanceamento também, e cultivares testes e muitas das cultivares que competiam tinham uma só repetição. Talvez você pode considerar o INLA nas suas análises também. À disposição. Walmes.
Pessoal, nao entendo muito desses desenhos experimentais mais complicados... :) Entao, corrigindo o que eu havia comentado: dois efeitos aleatorios variando "nos mesmos indices" podem ser identificaveis se tiverem estruturas de correlacao diferentes. Um detalhe: no codigo elegante do Walmes latitude e longitude sao corretamente colocados como vetores numericos, necessario para calcular distancias... On 12/19/2013 04:37 PM, walmes . wrote:
## Dentro de um bloco não se tem variação no valor de longitude. Isso ## pode implicar em problemas de identificabilidade dependendo do modelo ## considerado para os dados. O Elias já comentou isso.
O modelo de acordo com o que entendi é apenas um intercepto aleatório pros blocos e o efeito espacial, se é possível estimar na função gls() eu não sei, pelo que o Walmes fala é possível estruturar os resíduos sem problemas mas colocar mais um efeito aleatório de bloco acho que não. O Walmes pode confirmar ... Agora sem olhar seus dados o modelo teoricamente é estimável como é tudo gaussiano, você pode escrever a verossimilhança e maximizar numericamente, deve ser um problema bem estável. No nosso livro do SINAPE tem um exemplo de geoestatístico pode te ajudar a começar. Tudo que precisa fazer é somar mais uma matriz na matriz de covariâncias marginal. www.leg.ufpr.br/mcie Em 18 de dezembro de 2013 13:07, Luiz Roberto Martins Pinto < luizroberto.uesc@gmail.com> escreveu:
Penso que 'Blocos' não pode ser considerado como efeito fixo. Blocos são agrupamentos particulares das repetições. Portanto são repetições. Para que os testes estatísticos sejam válidos as repetições tem que ser aleatórias, caso contrário o experimento não tem valor, posto que as repetições são tendenciosas e as conclusões não podem ser inferidas para a população. Além disto o resíduo (aleatório) contém a interação blocos*tratamentos, que, portanto, deve ser aleatória. Note que no delineamento em que os tratamentos são aleatorizados nos blocos o gl do resíduo é (t-1)*(b-1). O resíduo é sempre aleatório porque o bloco é aleatório.
Luiz Roberto Martins Pinto Prof. Pleno/DCET/UESC Laboratório de Estatística Computacional Universidade Estadual de Santa Cruz Ilhéus-Bahia
luizroberto.uesc@gmail.com skype: lrmpinto http://lattes.cnpq.br/2732314327604831
Em 18 de dezembro de 2013 10:37, Wagner Bonat <wbonat@gmail.com> escreveu:
Por que não considera os blocos como efeito fixo ?
Você vai assumir que a estrutura do efeito espacial é sempre a mesma dentro dos blocos ?
Em 17 de dezembro de 2013 16:00, Cássio Dessotti < cassiodessotti@yahoo.com.br> escreveu:
Olá pessoal, boa tarde!!
Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R).
Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais.
1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R?
2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente.
Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti.
### Segue o código (com dados fictícios) para exemplificar o meu problema: local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp)
require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE)
### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)),
correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
Elias, Walmes e Wagner, antes de mais nada muito obrigado pela ajuda. Elias muito bem observado o problema de identificabilidade acarretado pelo confundimento de longitude com blocos, já que os blocos "estão em linha". Estudarei o material sugerido por você, Wagner. Walmes muito obrigado pelas linhas de comando que resolveram o modelo 3, com meus dados reais, já que no CMR coloquei apenas alguns dados para possibilitar a reprodução. Disponho de um grupo de experimentos em blocos ao acaso com 6 tratamentos, 5 blocos e 4 locais. A variável resposta é altura de plantas de cana-de-açúcar. Os blocos estão "em linha" o que acaba confundindo "longitude" com "blocos" mas sua sugestão para meu M3 (seu M4) resolveu muito bem. Pesquisarei esta literatura que você citou a respeito da questão de blocos ser considerado fixo ou aleatório. Walmes, construí diversos modelos alternando blocos (fixo, aleatório ou sem efeito), variâncias de locais (a mesma ou diferentes) e função espacial para a matriz R (exponencial, esférica e gaussiana), e optei pelo modelo exponencial sem efeito de blocos e com homogeneidade de variâncias como o "melhor" observando os valores de AIC. Você acha que o critério de AIC (ou BIC - apesar de nem sempre apontarem o mesmo modelo) é suficiente para a escolha de um modelo? Seria necessário o uso do critério de validação cruzada neste caso? A partir do modelo selecionado (análise conjunta), posso extrair os resíduos condicionais estudentizados, separá-los segundo locais, para em cada local construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar krigagem em cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Muito grato pela força, Cássio Dessotti. Em Quinta-feira, 19 de Dezembro de 2013 18:35, Wagner Bonat <wbonat@gmail.com> escreveu: O modelo de acordo com o que entendi é apenas um intercepto aleatório pros blocos e o efeito espacial, se é possível estimar na função gls() eu não sei, pelo que o Walmes fala é possível estruturar os resíduos sem problemas mas colocar mais um efeito aleatório de bloco acho que não. O Walmes pode confirmar ... Agora sem olhar seus dados o modelo teoricamente é estimável como é tudo gaussiano, você pode escrever a verossimilhança e maximizar numericamente, deve ser um problema bem estável. No nosso livro do SINAPE tem um exemplo de geoestatístico pode te ajudar a começar. Tudo que precisa fazer é somar mais uma matriz na matriz de covariâncias marginal. www.leg.ufpr.br/mcie Em 18 de dezembro de 2013 13:07, Luiz Roberto Martins Pinto <luizroberto.uesc@gmail.com> escreveu: Penso que 'Blocos' não pode ser considerado como efeito fixo.
Blocos são agrupamentos particulares das repetições. Portanto são repetições. Para que os testes estatísticos sejam válidos as repetições tem que ser aleatórias, caso contrário o experimento não tem valor, posto que as repetições são tendenciosas e as conclusões não podem ser inferidas para a população. Além disto o resíduo (aleatório) contém a interação blocos*tratamentos, que, portanto, deve ser aleatória. Note que no delineamento em que os tratamentos são aleatorizados nos blocos o gl do resíduo é (t-1)*(b-1). O resíduo é sempre aleatório porque o bloco é aleatório.
Luiz Roberto Martins Pinto Prof. Pleno/DCET/UESC Laboratório de Estatística Computacional Universidade Estadual de Santa Cruz Ilhéus-Bahia
luizroberto.uesc@gmail.com skype: lrmpinto http://lattes.cnpq.br/2732314327604831
Em 18 de dezembro de 2013 10:37, Wagner Bonat <wbonat@gmail.com> escreveu:
Por que não considera os blocos como efeito fixo ?
Você vai assumir que a estrutura do efeito espacial é sempre a mesma dentro dos blocos ?
Em 17 de dezembro de 2013 16:00, Cássio Dessotti <cassiodessotti@yahoo.com.br> escreveu:
Olá pessoal, boa tarde!!
Estou trabalhando com um grupo de experimentos em blocos ao acaso, com dados georreferenciados (x,y), buscando comparar modelos mistos (usuais) com modelos mistos que levam em conta a dependência espacial, utilizando de funções geoestatísticas na matriz de resíduos (R).
Consigo construir o modelo M1 (BA) com efeito de blocos aleatório a partir da função "lme", além do modelo M2 (Exp-H) que desconsidera o efeito de blocos, e considera uma função exponencial na matriz de covariâncias para os resíduos, além de considerar variâncias diferentes para os 2 locais.
1) A primeira questão está no modelo M3 (BA-Exp-H) que deve "unir" as características dos dois anteriores, considerando efeito de blocos aleatório, função exponencial para a matriz de covariâncias, além de heterogeneidade de variâncias, porém não consegui tal realização nem por meio da função "lme", nem por "gls". Seria necessário o uso de uma outra função? Como eu poderia trabalhar com este modelo no R?
2) A segunda questão é o seguinte: a partir do modelo que eu selecionar (provavelmente por AIC), desejo extrair os resíduos condicionais estudentizados, e separá-los segundo seus locais, para que em cada local eu possa construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar a krigagem de cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos). Posso seguir esta linha? Acredito que é um ganho de informação trabalhar com todos os locais (análise conjunta) neste sentido, ao invés de fazer todo esse processo em cada local separadamente.
Muito obrigado desde já pela atenção. Abraço a todos. Cássio Dessotti.
### Segue o código (com dados fictícios) para exemplificar o meu problema: local <- as.factor(c(rep(1,6),rep(2,6))); blo <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) lat <- as.factor(rep(1:3,4)) ; long <- as.factor(rep(c(rep(1,3),rep(2,3)),2)) blo_local <- as.factor(c(rep("1_1",3),rep("2_1",3),rep("1_2",3),rep("2_2",3))) trat <- as.factor(c(2,3,1,1,3,2,1,3,2,3,2,1)) resp <- c(2,1.9,1.9,1.8,2,1.9,2,2.1,2.3,2,1.9,1.8) dados <- data.frame(local,trat,blo,lat,long,blo_local,resp)
require(nlme) ### M1 - BA (efeito de blocos aleatório) M1 <- lme(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), method = "REML", na.action = na.omit, data=dados, keep.data=FALSE)
### M2 - Exp H (sem efeito de blocos - função exponencial para R - heterogeneidade de variâncias nos locais) M2 <- gls(resp ~ 1 + trat + local + local:trat, weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form=~as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
### M3 - BA-Exp-H (PROBLEMA - blocos aleatórios - função exponencial - heterogeneidade de variâncias) M3 <- gls(resp ~ 1 + trat + local + local:trat, random = list(blo_local = pdIdent(~1)), weight=varComb(varIdent(form = ~ 1|local)), correlation=corExp(form = ~ as.numeric(as.character(lat))+as.numeric(as.character(long))|local, metric="euclidean", nugget=FALSE), method="REML", na.action=na.omit, data=dados)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
On 12/20/2013 09:53 AM, Cássio Dessotti wrote:
Você acha que o critério de AIC (ou BIC - apesar de nem sempre apontarem o mesmo modelo) é suficiente para a escolha de um modelo? Seria necessário o uso do critério de validação cruzada neste caso?
R2, AIC, BIC, DIC, WBIC basicamente medem quao bom seu modelo esta ajustado aos dados. Medidas baseadas em validacao cruzada (aqui entra diversas, inclusive WAIC) servem para medir quao bom de predicao seu modelo e' ou se ha observacoes discrepantes.
A partir do modelo selecionado (análise conjunta), posso extrair os resíduos condicionais estudentizados, separá-los segundo locais, para em cada local construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar krigagem em cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos).
porque estimar um modelo incluindo efeito espacial e depois verificar dependencia espacial nos residuoes desse modelo? Vc pode pensar sob a abordagem de model based geostatistics. nesta abordagem efeitos fixos (delineamento) e aleatorios (espacial neste caso) sao considerados conjuntamente. Se o modelo ja levou em conta a estrutura espacial os residuos devem ser iid. A menos que se tenha mais de uma estrutura espacial ou e a estrutura incluida no modelo nao seja suficiente.
Eu usaria apenas a verossimilhança, no caso de modelos aninhados teste de razão de verossimilhança. A maioria dos outros são apenas mudanças da verossimilhança. Tambem não entendi pq usar semivariograma nos resíduos, vc já não controlou o efeito espacial ??? Vc pode ver um mapa apenas do efeito espacial já filtrando efeito de tratamento e bloco apenas com seu modelo ajustado. Eu incluiria na sua comparação o modelo matern com parâmetro kappa fixado em 1.5 e 2.5 seriam modelos intermediários ao exponencial e gaussiano que são casos particulares da matern com kappa 0.5 e infinito. Em 20 de dezembro de 2013 07:16, Elias T. Krainski < eliaskrainski@yahoo.com.br> escreveu:
On 12/20/2013 09:53 AM, Cássio Dessotti wrote:
Você acha que o critério de AIC (ou BIC - apesar de nem sempre apontarem o mesmo modelo) é suficiente para a escolha de um modelo? Seria necessário o uso do critério de validação cruzada neste caso?
R2, AIC, BIC, DIC, WBIC basicamente medem quao bom seu modelo esta ajustado aos dados. Medidas baseadas em validacao cruzada (aqui entra diversas, inclusive WAIC) servem para medir quao bom de predicao seu modelo e' ou se ha observacoes discrepantes.
A partir do modelo selecionado (análise conjunta), posso extrair os resíduos condicionais estudentizados, separá-los segundo locais, para em cada local construir semivariogramas destes resíduos, verificando a "força" da dependência, e realizar krigagem em cada área (trabalhando com os resíduos para "limpar" os dados de seus efeitos de tratamentos e demais efeitos).
porque estimar um modelo incluindo efeito espacial e depois verificar dependencia espacial nos residuoes desse modelo? Vc pode pensar sob a abordagem de model based geostatistics. nesta abordagem efeitos fixos (delineamento) e aleatorios (espacial neste caso) sao considerados conjuntamente. Se o modelo ja levou em conta a estrutura espacial os residuos devem ser iid. A menos que se tenha mais de uma estrutura espacial ou e a estrutura incluida no modelo nao seja suficiente.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
Wagner e Elias já contribuíram no sentido dos testes e gráficos dos resíduos. Eu só tenho uma observação a fazer que á a seguinte. Será que o bloco pode ser omitido na análise? Bem, quando trabalhamos em modelos de regressão em situações que não são de experimentos planejamos, é comum testarmos e abandonarmos variáveis do modelo. Numa situação de experimento planejado um pouco de cuidado deve ser tomado. Nos meus cursos eu já fui questionado sobre isso. Vamos contextualizar. Imagine um experimento fatorial com dois fatores (A e B) em delineamento de blocos. Suponha que a interação seja não significativa. O modelo maior considerado para os dados é E(y) = constante+bloco+A+B+A:B, Se a interação é nula, então eu caio no modelo logo acima desse que é E(y) = constante+bloco+A+B. Predição de valores e testes de médias, por exemplo, seriam baseados nesse modelo menor. Os pesquisadores não foram favoráveis a isso, diziam que não se pode fazer isso porque a soma de quadrados (SQ) da interação se junta ao resíduo, deixando ele maior, perdendo poder. Não exatamente. Aumenta-se SQ, mas aumenta-se os graus de liberdade, o quadrado médio pouco se altera. Então se o efeito de blocos for nulo pode-se abandonar do modelo? Bem, discutindo com outras pessoas, conclui que isso não deve ser feito pelo fato de bloco não representar tão somente um efeito no modelo mas também uma restrição à casualização na hora de montar o experimento. Por isso acredito que esse termo não deveria ser abandonado. Na minha humilde opinião, em pesquisa em agronomia, os blocos têm sido usados mais como maneira de deixar o experimento mais organizado (no sentido de layout mesmo) do que cumprir o papel de blocos. Não que esse seja seu caso. Bem, de qualquer forma, minha sugestão é de manter os blocos tem todos os modelos. À disposição. Walmes.
Senhores, novamente meu muito obrigado pela ajuda, pelas diversas sugestões e discussões que foram levantadas, e pelo material indicado... tenho bastante conteúdo a analisar que foi abordado nestes emails. Wagner, este mapa que você cita que poderia ser feito com os resíduos, não é o mapa de krigagem né? Verificarei o modelo matern que você citou... muito bem lembrado!! Walmes, achei muito interessante seu ponto de vista a respeito de retirar o fator blocos do modelo. Realmente, se o experimento foi planejado com os blocos, eles possivelmente tem sua função no modelo e devem ser "significativos" e relevantes. Considerarei sua sugestão (observação), optando então pelo modelo - dentre os que possuem o efeito de blocos - com menor AIC. Elias, abusando um pouco mais da sua boa vontade, você sabe algum material interessante de validação cruzada que eu possa estudar? Existe alguma função no R que você poderia me indicar para analisar o meu modelo escolhido? E por fim, você considera mais interessante eu aplicar o critério de validação apenas no modelo escolhido, ou em todos os modelos trabalhados? Muito obrigado a todos por todo o auxílio e atenção. Grande abraço. Cássio Dessotti. Em Sexta-feira, 20 de Dezembro de 2013 11:52, walmes . <walmeszeviani@gmail.com> escreveu: Wagner e Elias já contribuíram no sentido dos testes e gráficos dos resíduos. Eu só tenho uma observação a fazer que á a seguinte. Será que o bloco pode ser omitido na análise? Bem, quando trabalhamos em modelos de regressão em situações que não são de experimentos planejamos, é comum testarmos e abandonarmos variáveis do modelo. Numa situação de experimento planejado um pouco de cuidado deve ser tomado. Nos meus cursos eu já fui questionado sobre isso. Vamos contextualizar. Imagine um experimento fatorial com dois fatores (A e B) em delineamento de blocos. Suponha que a interação seja não significativa. O modelo maior considerado para os dados é E(y) = constante+bloco+A+B+A:B, Se a interação é nula, então eu caio no modelo logo acima desse que é E(y) = constante+bloco+A+B. Predição de valores e testes de médias, por exemplo, seriam baseados nesse modelo menor. Os pesquisadores não foram favoráveis a isso, diziam que não se pode fazer isso porque a soma de quadrados (SQ) da interação se junta ao resíduo, deixando ele maior, perdendo poder. Não exatamente. Aumenta-se SQ, mas aumenta-se os graus de liberdade, o quadrado médio pouco se altera. Então se o efeito de blocos for nulo pode-se abandonar do modelo? Bem, discutindo com outras pessoas, conclui que isso não deve ser feito pelo fato de bloco não representar tão somente um efeito no modelo mas também uma restrição à casualização na hora de montar o experimento. Por isso acredito que esse termo não deveria ser abandonado. Na minha humilde opinião, em pesquisa em agronomia, os blocos têm sido usados mais como maneira de deixar o experimento mais organizado (no sentido de layout mesmo) do que cumprir o papel de blocos. Não que esse seja seu caso. Bem, de qualquer forma, minha sugestão é de manter os blocos tem todos os modelos. À disposição. Walmes. _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
participantes (5)
-
Cássio Dessotti -
Elias T. Krainski -
Luiz Roberto Martins Pinto -
Wagner Bonat -
walmes .