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

walmes . walmeszeviani em gmail.com
Quinta Dezembro 19 13:37:43 BRST 2013


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


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