[R-br] Teste de médias em GLMM com covariável!

Ivan Bezerra Allaman ivanalaman em yahoo.com.br
Quinta Maio 10 08:13:20 BRT 2012


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: ivanalaman em yahoo.com.br/ivanalaman em gmail.com
@
\end{signature}


________________________________
 De: Walmes Zeviani <walmeszeviani em gmail.com>
Para: Ivan Bezerra Allaman <ivanalaman em 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: walmes em ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120510/9aadfb5b/attachment.html>


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