<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Alexandre,<br><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Precisa tomar cuidado com a parametrização usualmente imposta em modelo lineares e generalizados para estimar os parâmetros. No R o primeiro nível, por padrão, é assumido como tendo efeito zero, é o contraste tipo tratamento (contr.treatment). É preciso reconhecer isso para montar as matrizes de contraste corretamente. No código abaixo ilustra um pouco disso.<br>
<br><span style="font-family:courier new,monospace">##-----------------------------------------------------------------------------<br><br>require(multcomp)<br>require(doBy)<br>require(lattice)<br><br>## Dados artificiais - Status (Atacado x Sadio) x Idade(48 x 72)<br>
d <- data.frame(a=factor(sample(c('atacado','sadio'), 100, rep=TRUE)),<br>                b=factor(sample(c('48','72'), 100, rep=TRUE)))<br>xtabs(~a+b, d)<br><br>## Crio a interação.<br>
d$ab <- interaction(d$a, d$b)<br><br>X <- model.matrix(~ab, d)<br>colnames(X)<br>beta <- c(1,1.2,1,-0.9)<br><br>d$y <- rpois(nrow(d), lambda=exp(X%*%beta))<br><br>xyplot(y~a, groups=b, data=d, jitter.x=TRUE, type=c("p","a"))<br>
<br>## Ajusto o GLM de Poisson.<br>glmRes <- glm(y~ab, family="poisson", data=d)<br>coef(glmRes)<br><br>## Para específicar contrastes tem que tomar cuidado com a<br>## parametrização.<br><br>## Específicos os contrastes desejados.<br>
cntrMat <-   rbind("atacado.48-sadio.48"=c(1, -1,  0,  0),<br>                   "atacado.72-sadio.72"=c(0,  0,  1, -1),<br>                   "atacado.48-atacado.72"=c(1, 0, 0, -1),<br>                   "sadio.48-sadio.72"=c(0, 1, 0,  -1))<br>
<br>## Testo os contrastes com função glht do pacote multcomp.<br>summary(glht(glmRes, linfct=mcp(ab=cntrMat)), test=adjusted("none"))<br><br>summary(glht(glmRes, linfct=cntrMat, alternative="two.sided"),<br>
         test=adjusted("none"))<br><br>## Pela sua lógica, usando essa matriz deveria retornar as médias (na<br>## escala do preditor linear) para cada tratamento. Mas não são as médias.<br>M <- diag(4)<br>summary(glht(glmRes, linfct=mcp(ab=M)), test=adjusted("none"))<br>
summary(glht(glmRes, linfct=M), test=adjusted("none"))<br>coef(glmRes)<br><br>## Matriz com o coeficientes, já respeitando a parametrização, para<br>## obter as médias ajustadas na escala do preditor linear<br>M <- LSmatrix(glmRes, effect="ab")<br>
summary(glht(glmRes, linfct=M), test=adjusted("none"))<br><br>## Estimado vs paramétrico.<br>M%*%coef(glmRes)<br>M%*%beta<br><br>## Se quer facilitar a especificação dos constrastes, declare o modelo<br>## sem restrições, por exemplo, usando y~0+ab na fórmula.<br>
<br>glmRes <- glm(y~0+ab, family="poisson", data=d)<br>coef(glmRes)<br><br>cntrMat <-   rbind("atacado.48-sadio.48"=c(1, -1,  0,  0),<br>                   "atacado.72-sadio.72"=c(0,  0,  1, -1),<br>
                   "atacado.48-atacado.72"=c(1, 0, 0, -1),<br>                   "sadio.48-sadio.72"=c(0, 1, 0,  -1))<br><br>summary(glht(glmRes, linfct=mcp(ab=cntrMat)), test=adjusted("none"))<br>
summary(glht(glmRes, linfct=cntrMat), test=adjusted("none"))<br><br>##-----------------------------------------------------------------------------<br></span><br>À disposição.<br>Walmes.<br></div>​</div>