[R-br] Funções lme e gls - modelos mistos

Cássio Dessotti cassiodessotti em yahoo.com.br
Terça Dezembro 17 16:00:58 BRST 2013


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)
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131217/11d7dba4/attachment.html>


Mais detalhes sobre a lista de discussão R-br