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: [email protected]/[email protected]
@
\end{signature}
________________________________
De: Walmes Zeviani <[email protected]>
Para: [email protected]; Ivan Bezerra Allaman <[email protected]>
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: [email protected]
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================_______________________________________________
R-br mailing list
[email protected]
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.
