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

[walmes .]

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