[R-br] Contraste com glm com distribuição de poisson em um fatorial 2x2
ASANTOS
alexandresantosbr em yahoo.com.br
Sexta Maio 9 00:00:30 BRT 2014
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 em yahoo.com.br
alexandre.santos em 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 em 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.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140508/9df5c663/attachment.html>
Mais detalhes sobre a lista de discussão R-br