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

[Walmes Zeviani]

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/50f0f90c/attachment.html>