Walmes,

    Foi muito esclarecedor, mas suscitou novas dúvidas quando trabalho com GLM e são duas, sendo: 1) a ordem dos coeficientes em glm e a ordem dos níveis da interação e  2) se eu tivesse um fatorial 2x3x2, sendo as minhas indagações levantadas no CRM abaixo,

#1) Ordem dos coeficientes em glm e a ordem dos níveis da interação

#CRM

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)))

beta <- c(1,1.2,1,-0.9)
d$y <- rpois(nrow(d), lambda=exp(X%*%beta))

## Crio a interação.
d$ab <- interaction(d$a, d$b)

## Ordem dos níveis
levels(d$ab)

## Ajuste do GLM
glmRes <- glm(y~0+ab, family="poisson", data=d)
coef(glmRes)


Minha dúvida é sempre que ajustar um glm a ordem dos coeficientes, vai ser a mesma que a dos meus levels(d$ab) de forma que a parametrização de meu contraste for o for o primeiro nível "atacado.48"e o segundo "sadio.48" vou especificar o contrate (1, -1, 0, 0) e a matriz resultante vai retornar a média e para os outros níveis sucessivamente?


#2)Se eu tivesse um fatorial 2x3x2

## Dados artificiais - Status (Atacado x Sadio) x Idade(48 x 72 x 96) x Clone (TY01 e VM06)
d <- data.frame(a=factor(sample(c('atacado','sadio'), 100, rep=TRUE)),
                b=factor(sample(c('48','72','96'), 100, rep=TRUE)),
                c= factor(sample(c('TY01','VM06'), 100, rep=TRUE)))

xtabs(~a+b+c, d)
beta <- c(1,1.2,1,-0.9)
d$y <- rpois(nrow(d), lambda=exp(X%*%beta))

## Interação
d$abc<-interaction(d$a, d$b, d$c)

## Ajuste do GLM
glmRes2<- glm(y~0+abc, family="poisson", data=d)
coef(glmRes)

## Especifico os contrastes desejados, neste caso "atacado.48.TY01" e "sadio.48.TY01" (primeiro nível com o segundo)
cntrMat <-   rbind("atacado.48.TY01-sadio.48.TY01"=c(1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

summary(glht(glmRes2, linfct=mcp(abc=cntrMat)), test=adjusted("none"))


E pergunto se esta especificação de coeficientes esta correta no caso de um fatorial 2x3x2, imaginado a análise em one-way?


Novamente obrigado,

-- 
======================================================================
Alexandre dos Santos
Proteção Florestal 
IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
Campus Cáceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial 
Cáceres - MT                      CEP: 78.200-000
Fone: (+55) 65 8132-8112 (TIM)   (+55) 65 9686-6970 (VIVO)
e-mails:alexandresantosbr@yahoo.com.br 
        alexandre.santos@cas.ifmt.edu.br 
Lattes: http://lattes.cnpq.br/1360403201088680                   
======================================================================



Em 08/05/2014 08:27, walmes . escreveu:
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.
​


_______________________________________________
R-br mailing list
R-br@listas.c3sl.ufpr.br
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.