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

Ivan Bezerra Allaman ivanalaman em yahoo.com.br
Quarta Maio 9 21:17:24 BRT 2012


Boa noite Walmes, obrigado pela pronta resposta!

É provável que não estejamos falando do mesmo tipo de ANCOVA (Análise de covariância). No entanto, como você está anos luz a frente, solucione algumas dúvidas primeiro para que eu não fale besteira. Eu não entendi a seguinte pedaço de frase:``você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F.'' O tipo de ANCOVA que estou falando, não tem sentido comparar níveis de um grupo para um determinado nível de covariável, afinal de contas ela é quantitativa e sua inclusão foi muito bem explicada por você nos trechos acima. O fato de incluirmos a interação entre a covariável e o fator é devido a intenção de avaliar se será necessário ajustar as médias considerandos slopes diferentes (caso dê significativo), ou se considera um slope comum na correção das médias. Até cheguei cogitar de no ítem `` médias populacionais marginais em x=0.5'', você ter considerado x=0.5 como um slope comum para ajuste
 das médias, mais pelos resultados, creio que não foi essa a intenção. Vamos pegar o gancho do seu CRM, e discutindo ao longo das análises para que eu possa me fazer entender.

Considerando já todo o enredo brilhantemente explicado por ti de acordo com o CRM temos o seguinte modelo:

ac0 <- glm(y~F*x, data=da, family=poisson)
anova(ac0, test="Chisq")

De acordo com os resultados, há a necessidade de se ajustar as médias considerando uma regressão diferente para cada grupo. E de fato, os slopes são diferentes estimados para cada grupo.

Graficamente, de fato:

require(lattice)
xyplot(y~x, groups=F, data=da, type="a")

Analiticamente, de fato:

slopes <- list()
for(i in 1:5){
    slopes[[i]] <- coef(glm(y ~ x, data=subset(da,F==i)))[2]
}
slopes

Segundo a literatura, o ideal é que ao se comparar as médias, estas devem estar ajustadas para as covariáveis (caso ela for significativa é claro). Então ao invés de eu comparar estas médias:

with(da, tapply(y, F, mean))

e devo comparar estas médias:

novamedia_F1 <- with(da,mean(y[F==1])) - (exp(as.numeric(slopes[1]))*(with(da,mean(x[F==1]))-mean(da$x)))
novamedia_F2 <- with(da,mean(y[F==2])) - (exp(as.numeric(slopes[2]))*(with(da,mean(x[F==2]))-mean(da$x)))
novamedia_F3 <- with(da,mean(y[F==3])) - (exp(as.numeric(slopes[3]))*(with(da,mean(x[F==3]))-mean(da$x)))
novamedia_F4 <- with(da,mean(y[F==4])) - (exp(as.numeric(slopes[4]))*(with(da,mean(x[F==4]))-mean(da$x)))
novamedia_F5 <- with(da,mean(y[F==5])) - (exp(as.numeric(slopes[5]))*(with(da,mean(x[F==5]))-mean(da$x)))

rbind(novamedia_F1,novamedia_F2,novamedia_F3,novamedia_F4,novamedia_F5)

E a grande pergunta: Como eu posso fazer isso no R?


 

\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: r-br em listas.c3sl.ufpr.br; Ivan Bezerra Allaman <ivanalaman em yahoo.com.br> 
Enviadas: Quarta-feira, 9 de Maio de 2012 16:05
Assunto: Re: [R-br] Teste de médias em GLMM com covariável!
 

Ivan,

Vamos primeiro primeiro apodar as arestas a fim de fazer a coisa rolar. Você tem um experimento com um fator de 5 níveis nominais, denominado de F. Cada animal é uma unidade experimental do qual se observa uma resposta y e uma covariável x. Acredita-se que x afeta y e portanto quer-se remover esse efeito para poder estudar a variação explicada por F. Pode-se supor que F e x tenham efeitos aditivos (sem interação) ou pode-se assumir que existe interação. Na ausência de interação as diferenças entre níveis de F são as mesmas sem importar o nível de x, mas o mesmo não vale para presença de interação, de forma que, você deve definir antecipadamente em que nível(is) de x quer comparar níveis de F. Não sei o CMR é exatamente o que você tem, pois o CMR não parece ser um caso de modelo misto porque animal é a unidade experimental e não um fator de agrupamento (embora para Poisson este seja um modelo estimável), enfim... vou
 considerar um glm fixo, daí é só fazer as devidas extensões para o glmm.

da <- data.frame(F=gl(5,20), x=runif(5*20))
da$Ey <- with(da, exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x))
da$y <- with(da, rpois(F, lambda=exp(1+as.numeric(F)/2+(as.numeric(F)-2)/3*x)))
with(da, tapply(y, F, mean))

require(lattice)
xyplot(y~x, groups=F, data=da, type="a")

ac0 <- glm(y~F*x, data=da, family=poisson)
anova(ac0, test="Chisq")
betas <- coef(ac0)

require(doBy)
# médias populacionais marginais em x=0.5
mpm <- popMatrix(ac0, effect="F", at=list(x=0.5))
mpm%*%betas
exp(mpm%*%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) mpm[vs[1],]-mpm[vs[2],]))
rownames(dmpm) <- rownames(comp)

apply(dmpm, 1, function(x) x%*%betas)

require(multcomp)
summary(glht(ac0, linfct=mpm))  # médias populacionais marginais
summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais

# na verdade isso não é a média, mas sim o exp() disso

E é isso aí Ivan. Esse procedimento aí que fiz foi o que deixou um determinado software famoso por dar esses resultados aos usuários. Trata-se de algo que, depois de fazer pela primeira vez e entender, passa a ser trivial.

À 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/20120509/c28e9405/attachment.html>


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