
Alexandre, 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. ##----------------------------------------------------------------------------- require(multcomp) require(doBy) require(lattice) ## Dados artificiais - Status (Atacado x Sadio) x Idade(48 x 72) d <- data.frame(a=factor(sample(c('atacado','sadio'), 100, rep=TRUE)), b=factor(sample(c('48','72'), 100, rep=TRUE))) xtabs(~a+b, d) ## Crio a interação. d$ab <- interaction(d$a, d$b) X <- model.matrix(~ab, d) colnames(X) beta <- c(1,1.2,1,-0.9) d$y <- rpois(nrow(d), lambda=exp(X%*%beta)) xyplot(y~a, groups=b, data=d, jitter.x=TRUE, type=c("p","a")) ## Ajusto o GLM de Poisson. glmRes <- glm(y~ab, family="poisson", data=d) coef(glmRes) ## Para específicar contrastes tem que tomar cuidado com a ## parametrização. ## Específicos os contrastes desejados. cntrMat <- rbind("atacado.48-sadio.48"=c(1, -1, 0, 0), "atacado.72-sadio.72"=c(0, 0, 1, -1), "atacado.48-atacado.72"=c(1, 0, 0, -1), "sadio.48-sadio.72"=c(0, 1, 0, -1)) ## Testo os contrastes com função glht do pacote multcomp. summary(glht(glmRes, linfct=mcp(ab=cntrMat)), test=adjusted("none")) summary(glht(glmRes, linfct=cntrMat, alternative="two.sided"), test=adjusted("none")) ## Pela sua lógica, usando essa matriz deveria retornar as médias (na ## escala do preditor linear) para cada tratamento. Mas não são as médias. M <- diag(4) summary(glht(glmRes, linfct=mcp(ab=M)), test=adjusted("none")) summary(glht(glmRes, linfct=M), test=adjusted("none")) coef(glmRes) ## Matriz com o coeficientes, já respeitando a parametrização, para ## obter as médias ajustadas na escala do preditor linear M <- LSmatrix(glmRes, effect="ab") summary(glht(glmRes, linfct=M), test=adjusted("none")) ## Estimado vs paramétrico. M%*%coef(glmRes) M%*%beta ## Se quer facilitar a especificação dos constrastes, declare o modelo ## sem restrições, por exemplo, usando y~0+ab na fórmula. glmRes <- glm(y~0+ab, family="poisson", data=d) coef(glmRes) cntrMat <- rbind("atacado.48-sadio.48"=c(1, -1, 0, 0), "atacado.72-sadio.72"=c(0, 0, 1, -1), "atacado.48-atacado.72"=c(1, 0, 0, -1), "sadio.48-sadio.72"=c(0, 1, 0, -1)) summary(glht(glmRes, linfct=mcp(ab=cntrMat)), test=adjusted("none")) summary(glht(glmRes, linfct=cntrMat), test=adjusted("none")) ##----------------------------------------------------------------------------- À disposição. Walmes.