[R-br] Blocos com algumas repetições do material

Sábado Outubro 20 20:59:57 BRT 2012

Éder,

Na simulação faltou o efeito de bloco. Bem, nas formulas para efeito
aleatório da nlme e lme4 nos temos que declarar usando termos que vão antes
de depois da barra em pé | (pipe). Antes da barra nos colocamos um fator de
efeito fixo e após um de efeito aleatório. Assim, essas expecificações tem
a seguinte interpretação

* (1|local) : 1 representa o parâmetro beta_0 ou intercepto e com isso
dizemos que os níveis de local são de efeito aleatório e incidem no
intercepto, ou seja, esses desvios somam-se ao intercepto e existe um para
cada nível de local;
* (1|local)+(x|local) : considerando x numérico (ex: dose, tempo) estamos
dizendo que o efeito aleatório de local incide em b_0 e b_1 que o
coeficiente angular, taxa de incremento por unidade de x, e isso é o modelo
de intercepto e inclinação aleatório, note que x tem efeito fixo.
* (1|local)+(trat|local): considerando trat como fator, isso significa que
o efeito aleatório de local causa desvios em cada um dos parâmetros
associados aos níveis de trat (efeito fixo). Dificilmente você ira declarar
uma formular com um fator antes do pipe, aqui será estimado uma variância
para bloco dentro de cada nível de trat, pouco usado. Na prática poderia-se
ter (ano|local), que significa que para ano o efeito de local tem uma
variância diferente (demora bastante para estimar).
* (1|local:bloco): com isso você declara que existe uma fonte de variação
cujos níveis são as combinações (produto cartesiano) dos níveis dos fatores
envolvidos e que o efeito deste incide no intercepto.

Bem, com isso o termo que você incluiu (0+Rep|Local) não tem sentido
olhando para sua descrição do experimento, embora coisas sejam estimadas.
Na realidade, se você repete níveis de um fator dentro do bloco você
aumenta a precisão na estimativa daqueles repetidos e aumenta a informação
para estimação do erro resídual. Você não precisa declarar essa fonte de
variação, ela "sai por diferença". Como é o ultimo termo de efeito
aleatório (que todo modelo tem) não precisa ser declarado. Então acredito
que você não esteja definindo corretamente o modelo que você descreveu. Se
você muda a estrutura do efeito aleatório algum impacto terá nas
estimativas dos efeitos fixos, por isso que m0 e m1 diferem. No CMR abaixo
eu simulo do modelo que entendi da sua descrição e faço ajustes de alguns
modelos alternativos.

#------------------------------------------------------------------------------------------
# Experimento com efeito de locais, locais/bloco e variedade
# dados artificiais (usando muitos níveis para estimativas ficarem próximas
dos valores
# usados na simulação, facilita depuração)
nl <- 20; nb <- 5; nc <- 12; nr <- 3; nu <- nc+nr
# efeito de local ~N(0, 2), 20 locais
l <- rnorm(nl, 0, 2)
# efeito de local/bloco ~N(0, 1), 5 blocos por local
lb <- rnorm(nl*nb, 0, 1)
# efeito de local/bloco/unidade ~N(0, 0.5), 8 unidades por bloco
lbu <- rnorm(nl*nb*nu, 0, 0.5)
# efeito de cultivar ~N(0, 1.5), 5 cultivares
c <- rnorm(nc, 0, 1.5)
da1 <- expand.grid(ano=gl(2,1), local=gl(nl,1), bloco=gl(nb,1),
cultivar=gl(nc,1), r=1)
da2 <- expand.grid(ano=gl(2,1), local=gl(nl,1), bloco=gl(nb,1),
cultivar=gl(nr,1), r=2)
da <- rbind(da1, da2)
Z <- with(da, cbind(model.matrix(~-1+local),
                    model.matrix(~-1+local:bloco),
                    model.matrix(~-1+cultivar)))
dim(Z)
ranef <- c(l, lb, c)
length(ranef)
X <- model.matrix(~ano, da)
y <- X%*%c(1,0)+Z%*%ranef+lbu

require(lme4)

# modelo declarado é o mesmo usado para simular
m0 <- lmer(y~ano+(1|local)+(1|local:bloco)+(1|cultivar), data=da)
summary(m0)

# local fixo, bloco dentro de local e cultivar aleatório
m1 <- lmer(y~ano+local+(1|local:bloco)+(1|cultivar), data=da)
summary(m1)

# cultivar fixo, bloco dentro de local aleatório
m2 <- lmer(y~ano+cultivar+(1|local/bloco), data=da) # local/bloco =
local+local:bloco
summary(m2) # cultivares mais repetidas, 1:3, tem menor erro padrão

# efeito de local/bloco com variâncias diferentes nos anos (omiti o outros
efeitos)
m3 <- lmer(y~(ano|local/local), data=da)
summary(m3) # cultivares mais repetidas, 1:3, tem menor erro padrão

#------------------------------------------------------------------------------------------

À disposição.
Walmes.

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes em ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================

2012/10/20 tiago souza marçal <tiagosouzamarcal em hotmail.com>

>  Não sei se eu entendi bem mas a sua duvida é saber qual dos modelos
> abaixo usar?
>
> Yijk = m + Gi + B/Ajk + Aj + GAij + Eijk
>
> Yijk = m + Gi + Bk + Aj + GAij + Eijk
>
> Att.
>
> Tiago.
>
> ------------------------------
> Date: Sat, 20 Oct 2012 15:36:31 -0300
> From: eder em leg.ufpr.br
> To: R-br em listas.c3sl.ufpr.br
> Subject: Re: [R-br] Blocos com algumas repetições do material
>
>
>
> Boa tarde Pessoal,
>
> Tenho experimentos em blocos(Locais) e dentro deles gostaria de por
> algumas de repetições de alguns tratamentos para estimar a variância dentro
> do bloco, que desconfio ser grande. Por exemplo, tenho 5 Cultivares
> gostaria de repetir 3 delas dentro do bloco, tendo ao todo 8 parcelas, o
> ideal seria ter uma repetição de tudo, porem na pratica isso não
> é possível.
> Minhas duvidas são:
> 1) Estou definindo de maneira correta essas repetições dentro do modelo?
> 2) Como o método  de calculo dos dois modelos abaixo é igual para os
> efeitos fixos, por que o intercepto esta com valores diferentes?
> 3) Neste dois casos(m0 e m1, abaixo) o BLUP seria a média mais o efeito da
> cultivar?
> Obrigado pela ajuda
>
> require(lme4)
> ### Simulando dados
> dados <-
> expand.grid(Variedade=factor(c(letters[1:3],letters[1:5])),Local=factor(1:5))
> dados$Rep <- factor(rep(c(2,2,2,1,1,1,1,1),t=5))
> dados$Resp <- sort(rnorm(nrow(dados)))
> dados
>
> ### Modelos
> # Sem efeito dentro do local
> m0 <- lmer(Resp~1+(1|Local)+(1|Variedade),dados)
> summary(m0)
> ranef(m0)
> # Com efeito dentro do local
> m1 <- lmer(Resp~1+(0+Rep|Local)+(1|Variedade),dados)
> summary(m1)
> ranef(m1)
>
>
> _______________________________________________ R-br mailing list
> R-br em 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 em 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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121020/81e62777/attachment.html>