[R-br] Contraste com glm com distribuição de poisson em um fatorial 2x2

[walmes .]

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.
​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140508/6eabe626/attachment.html>