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.