<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Alexandre,<br><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">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.<br>
<br><span style="font-family:courier new,monospace">##-----------------------------------------------------------------------------<br><br>require(multcomp)<br>require(doBy)<br>require(latticeExtra)<br><br>## apc: build all pairwise comparisons matrix by a LSmatrix.<br>
source("<a href="http://dl.dropboxusercontent.com/u/48140237/apc.R">http://dl.dropboxusercontent.com/u/48140237/apc.R</a>")<br><br></span><span style="font-family:courier new,monospace"><span style="font-family:courier new,monospace">##-----------------------------------------------------------------------------<br>
<br></span>d <- data.frame(a=factor(sample(c('atacado','sadio'), 100, rep=TRUE)),<br>                b=factor(sample(c('48','72','96'), 100, rep=TRUE)),<br>                c= factor(sample(c('TY01','VM06'), 100, rep=TRUE)))<br>
<br>xtabs(~a+b+c, d)<br><br>d$y <- rpois(nrow(d), lambda=10)<br><br>## Interação.<br>d$abc <- interaction(d$a, d$b, d$c)<br>levels(d$abc)<br><br>## O padrão de funcionamento faz assim, o primeiro fator troca níveis<br>
## dentro do segundo que troca níveis dentro do terceiro, etc, veja.<br>do.call(rbind, strsplit(levels(d$abc), "\\."))<br><br>## Sabendo disso você pode passar da a, b e c na forma que lhe for mais<br>## conveniente. Por outro lado, é desperdício de tempo operar assim<br>
## porque eu considero mais fácil declarar o modelo fatorial, usar a<br>## LSmatrix para gerar a matriz de contraste e usá-la.<br><br>g0 <- glm(y~a*b*c, family="poisson", data=d)<br>g1 <- glm(y~0+abc, family="poisson", data=d)<br>
<br>M <- LSmatrix(g0, effect=c("a","b","c"))<br>M<br><br>## As estimativas das médias (na escala do preditor linear).<br>data.frame(g0=M%*%coef(g0), g1=coef(g1))<br><br>## Combinações de níveis correspondentes as médias.<br>
str(M)<br>grid <- attr(M, "grid")<br><br></span><span style="font-family:courier new,monospace">## Matriz de contrastes entre níveis de `focus` dentro dos níveis de<br>## `split`.<br>split <- c("b","c")<br>
focus <- "a"<br>spl <- interaction(grid[,split])<br>i <- 1:nrow(grid)<br>l <- split(i, f=spl)<br>contr <- lapply(l,<br>                function(row){<br>                    ## Matriz de contrastes par a par.<br>
                    a <- apc(M[row,], lev=levels(d[,focus]))<br>                    ## Prefixo no nome das linhas.<br>                    rownames(a) <- paste(spl[row[1]],<br>                                         rownames(a), sep="/")<br>
                    return(a)<br>                })<br>contr <- do.call(rbind, contr)<br>contr<br><br>## Constrastes.<br>summary(glht(g0, linfct=contr),<br>        test=adjusted(type="fdr"))<br><br>## Para ter a contra prova da equivalência dos procedimentos.<br>
cntrMat <- rbind("atacado.48.TY01-sadio.48.TY01"=c(<br>                     1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))<br>summary(glht(g1, linfct=cntrMat),<br>        test=adjusted(type="fdr"))<br><br>##-----------------------------------------------------------------------------<br>
## No caso de desdobrar interação dupla.<br><br>M <- LSmatrix(g0, effect=c("a","b"))<br>grid <- attr(M, "grid")<br><br>## Matriz de contrastes entre níveis de `focus` dentro dos níveis de<br>
## `split`.<br>split <- c("b")<br>focus <- "a"<br>spl <- interaction(grid[,split])<br>i <- 1:nrow(grid)<br>l <- split(i, f=spl)<br>contr <- lapply(l,<br>                function(row){<br>
                    ## Matriz de contrastes par a par.<br>                    a <- apc(M[row,], lev=levels(d[,focus]))<br>                    ## Prefixo no nome das linhas.<br>                    rownames(a) <- paste(spl[row[1]],<br>
                                         rownames(a), sep="/")<br>                    return(a)<br>                })<br>contr <- do.call(rbind, contr)<br>contr<br><br>## Constrastes.<br>summary(glht(g0, linfct=contr),<br>
        test=adjusted(type="fdr"))<br><br>##-----------------------------------------------------------------------------<br>## Representando em um gráfico média com IC.<br><br>grid <- attr(M, "grid")<br>
grid <- sapply(grid, function(x) if(is.character(x)) factor(x) else x)<br>means <- confint(glht(g0, linfct=M), calpha=univariate_calpha())<br>grid <- cbind(data.frame(grid), means$confint)<br>str(grid)<br><br>segplot(a~lwr+upr|b, data=grid, draw=FALSE, centers=Estimate)<br>
<br>##-----------------------------------------------------------------------------<br><br></span>À disposição.<br>Walmes.<br></div></div>