Valeu Walmes, era exatamente isso que eu procurava. Agora entendi.
Muito obrigado!
\begin{signature}
<<>>=
Prof. Dr. Ivan Bezerra Allaman
Universidade Estadual de Santa Cruz
Departamento de Ciências Exatas e Tecnológicas
Ilhéus/BA - Brasil
Fone: +55 73 3680-5596
E-mail: ivanala...@yahoo.com.br/ivanala...@gmail.com
@
\end{signature}
________________________________
De: Walmes Zeviani <walmeszevi...@gmail.com>
Para: Ivan Bezerra Allaman <ivanala...@yahoo.com.br>
Enviadas: Quinta-feira, 10 de Maio de 2012 0:14
Assunto: Re: [R-br] Teste de médias em GLMM com covariável!
Ivan,
Bem, de fato não estou seguro sobre ter entendido o que você quer, mas vamos
lá. Partindo do fato que o nosso modelo aponta interação entre F e x nos
teremos retas não paralelas (na escala do preditor linear), portando diferenças
de médias não são constantes a medida que variamos o valor de x. Podemos
comparar as médias em um valor fixo de x. Por exemplo, podemos comparar o peso
final de suínos do sexo macho vs fêmea (F) supondo que os animais tem mesma
idade (x=100 dias). Outra possibilidade, que deve ser a que você quer, é
comparar nas médias amostrais de x nos níveis de F. Assim, comparamos o peso
final dos machos com idade média de 100 dias com as fêmeas com idade média
final de 95 dias. Se isso é uma comparação justa ou não é outra história e
depende do contexto. Eu fiz o CMR desse último caso. É isso que você procura?
Eu prefiro sempre que possível evitar as comparações (principalmente seguidas
de letras) e dar um
gráfico com curvas ajustadas acompanhadas das bandas de confiança. Assim o
leitor pode ver como são as diferenças em cada nível de x, entender a relação
de x com F e apontar para quais intervalos em x existe ou se pronuncia a
diferença entre os níveis de F. É sem dúvida mais informativo. Espero ter
ajudado.
pred <- expand.grid(F=levels(da$F), x=seq(0,1,l=30))
pred$eta <- predict(ac0, newdata=pred)
pred$y <- predict(ac0, newdata=pred, type="response")
require(lattice); require(latticeExtra)
# retas não paralelas. em que ponto vou compara-las?
# em cada lugar dá uma conclusão diferente!
xyplot(eta~x, groups=F, data=pred, type="l")
# você pode comparar todas em um valor fixo, como x=0.5
# ou pode comparar ela no valor médio de x para cada nível de F
xm <- with(da, tapply(x, F, mean)); xm
mpm==0.5
mpm2 <- mpm
for(i in 1:length(xm)) mpm2[i,mpm2[i,]==0.5] <- xm[i]
mpm2
mpm2%*%betas
exp(mpm2%*%betas)
# diferença nas médias populacionais marginais duas à duas
nlevels(da$F)
comp <- t(combn(1:nlevels(da$F), 2))
rownames(comp) <- apply(comp, 1, paste, collapse="-")
dmpm <- t(apply(comp, 1, function(vs) mpm2[vs[1],]-mpm2[vs[2],]))
rownames(dmpm) <- rownames(comp)
apply(dmpm, 1, function(x) x%*%betas)
require(multcomp)
summary(glht(ac0, linfct=mpm2)) # médias populacionais marginais
summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais
marginais
xyplot(y~x, groups=F, data=da, type="p")+
as.layer(xyplot(y~x, groups=F, data=pred, type="l"))
À 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: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================_______________________________________________
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.
