Contraste com glm com distribuição de poisson em um fatorial 2x2
Boa noite Pessoal, Estou tentando realizar um contraste com glm com distribuição de poisson em um fatorial 2x2 e ao final testar se esta tudo OK, mas esbarei em algumas dúvidas a respeito da interpretação de dados com distribuição normal x poisson, em meu exemplo tenho: #--------------------------------------------------------------------------------------------- require(multcomp) ##Dados artificiais - Status (Atacado x Sadio) x Idade(48 x 72) d <- data.frame(a=factor(sample(c('atacado','sadio'), 100, rep=T)), b=factor(sample(c('48','72'), 100, rep=T))) d$y <- as.numeric(d$a)*rpois(100,5)+ as.numeric(d$b)*rpois(100,20)+ as.numeric(d$a)*as.numeric(d$b)*rpois(100, 15)+rpois(100, 40); ##Crio a interação d$ab <- interaction(d$a, d$b) ##Ajusto o GLM de poisson glmRes <- glm(y ~ ab, family="poisson", data=d) ##Especificos 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), alternative="two.sided"), test=adjusted("none")) #-------------------------------------------------------------------------------------------------------- Mas agora gostaria de testar manualmente os resultados do primeiro contraste, se fossem para uma variável resposta com distribuição normal aplicaria o teste de t e calcularia o p valor, mas para dados com distribuição de erros de Poisson? Alguém poderia me ajudar? 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 ======================================================================
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.
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.
Alexandre, Existe uma lógica que pode ser considerada para montar esses contrastes corretamente. A lógica não é difícil, por outro lado, montar as matrizes de contraste "na mão" é demorado, sujeito à erros, não prontamente adaptável à outros experimentos. Com a experiência eu desenvolvi um jeito mais mínimo reproduzível de fazer. Para experimento fatoriais com perda de caselas (combinação de níveis) esse código deve ser modificado. ##----------------------------------------------------------------------------- require(multcomp) require(doBy) require(latticeExtra) ## apc: build all pairwise comparisons matrix by a LSmatrix. source("http://dl.dropboxusercontent.com/u/48140237/apc.R") ##----------------------------------------------------------------------------- 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) d$y <- rpois(nrow(d), lambda=10) ## Interação. d$abc <- interaction(d$a, d$b, d$c) levels(d$abc) ## O padrão de funcionamento faz assim, o primeiro fator troca níveis ## dentro do segundo que troca níveis dentro do terceiro, etc, veja. do.call(rbind, strsplit(levels(d$abc), "\\.")) ## Sabendo disso você pode passar da a, b e c na forma que lhe for mais ## conveniente. Por outro lado, é desperdício de tempo operar assim ## porque eu considero mais fácil declarar o modelo fatorial, usar a ## LSmatrix para gerar a matriz de contraste e usá-la. g0 <- glm(y~a*b*c, family="poisson", data=d) g1 <- glm(y~0+abc, family="poisson", data=d) M <- LSmatrix(g0, effect=c("a","b","c")) M ## As estimativas das médias (na escala do preditor linear). data.frame(g0=M%*%coef(g0), g1=coef(g1)) ## Combinações de níveis correspondentes as médias. str(M) grid <- attr(M, "grid") ## Matriz de contrastes entre níveis de `focus` dentro dos níveis de ## `split`. split <- c("b","c") focus <- "a" spl <- interaction(grid[,split]) i <- 1:nrow(grid) l <- split(i, f=spl) contr <- lapply(l, function(row){ ## Matriz de contrastes par a par. a <- apc(M[row,], lev=levels(d[,focus])) ## Prefixo no nome das linhas. rownames(a) <- paste(spl[row[1]], rownames(a), sep="/") return(a) }) contr <- do.call(rbind, contr) contr ## Constrastes. summary(glht(g0, linfct=contr), test=adjusted(type="fdr")) ## Para ter a contra prova da equivalência dos procedimentos. cntrMat <- rbind("atacado.48.TY01-sadio.48.TY01"=c( 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) summary(glht(g1, linfct=cntrMat), test=adjusted(type="fdr")) ##----------------------------------------------------------------------------- ## No caso de desdobrar interação dupla. M <- LSmatrix(g0, effect=c("a","b")) grid <- attr(M, "grid") ## Matriz de contrastes entre níveis de `focus` dentro dos níveis de ## `split`. split <- c("b") focus <- "a" spl <- interaction(grid[,split]) i <- 1:nrow(grid) l <- split(i, f=spl) contr <- lapply(l, function(row){ ## Matriz de contrastes par a par. a <- apc(M[row,], lev=levels(d[,focus])) ## Prefixo no nome das linhas. rownames(a) <- paste(spl[row[1]], rownames(a), sep="/") return(a) }) contr <- do.call(rbind, contr) contr ## Constrastes. summary(glht(g0, linfct=contr), test=adjusted(type="fdr")) ##----------------------------------------------------------------------------- ## Representando em um gráfico média com IC. grid <- attr(M, "grid") grid <- sapply(grid, function(x) if(is.character(x)) factor(x) else x) means <- confint(glht(g0, linfct=M), calpha=univariate_calpha()) grid <- cbind(data.frame(grid), means$confint) str(grid) segplot(a~lwr+upr|b, data=grid, draw=FALSE, centers=Estimate) ##----------------------------------------------------------------------------- À disposição. Walmes.
Walmes, Muito mais eficiente que fazer na mão, problema resolvido e abordagem adotada. Obrigado, Alexandre Em 09/05/2014 11:14, walmes . escreveu:
##----------------------------------------------------------------------------- ## No caso de desdobrar interação dupla.
M <- LSmatrix(g0, effect=c("a","b")) grid <- attr(M, "grid")
## Matriz de contrastes entre níveis de `focus` dentro dos níveis de ## `split`. split <- c("b") focus <- "a" spl <- interaction(grid[,split]) i <- 1:nrow(grid) l <- split(i, f=spl) contr <- lapply(l, function(row){ ## Matriz de contrastes par a par. a <- apc(M[row,], lev=levels(d[,focus])) ## Prefixo no nome das linhas. rownames(a) <- paste(spl[row[1]], rownames(a), sep="/") return(a) }) contr <- do.call(rbind, contr) contr
## Constrastes. summary(glht(g0, linfct=contr), test=adjusted(type="fdr"))
##----------------------------------------------------------------------------- ## Representando em um gráfico média com IC.
grid <- attr(M, "grid") grid <- sapply(grid, function(x) if(is.character(x)) factor(x) else x) means <- confint(glht(g0, linfct=M), calpha=univariate_calpha()) grid <- cbind(data.frame(grid), means$confint) str(grid)
segplot(a~lwr+upr|b, data=grid, draw=FALSE, centers=Estimate)
##-----------------------------------------------------------------------------
-- ====================================================================== 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 ======================================================================
participantes (2)
-
ASANTOS -
walmes .