[R-br] Matrix de contrastes em fatorial DBC com a função glht!

Walmes Zeviani walmeszeviani em gmail.com
Terça Dezembro 15 18:59:22 BRST 2015


Ivan,

Você forneceu todo o código necessário para reprodução da sua dúvida mas
não deu as informações suficientes com relação ao experimento, etc. Na
possibilidade de fracassar com tantas suposições à respeito do que não foi
dito, eu adaptei código de scripts que eu já tinha sobre desdobramento de
interação com a glht()

Bem, um ponto é que desbalanceado ou não, os coeficientes considerados para
as médias ajustadas (e os contrates) são os mesmos. O desbalanceamento já é
considerado no tamanho dos erros padrões dos efeitos, onde àqueles de erro
padrão menor são efeitos com mais repetições, nesse caso. Existe outros
tipos de "ponderação" para obter essas médias, mas só use se realmente você
se sentir convencido da necessidade ou motivado por algum artigo (e não
documentação de $oftware$).

##-------------------------------------------

## As coisas ficam mais fáceis quando se compreende a estrutura
## experimental.

xtabs(~sujeito, data=gnexpo)
xtabs(~celula+tratamento, data=gnexpo)
ftable(xtabs(~celula+tratamento+sujeito, data=gnexpo))

sum(!complete.cases(gnexpo))

## Distribuição dos casos perdidos.
addmargins(
    with(gnexpo,
         tapply(ybox, list(celula, tratamento),
                FUN=function(x){
                    sum(is.na(x))
                })))

m0 <- lm(ybox~sujeito+celula*tratamento, data=gnexpo)
par(mfrow=c(2,2)); plot(m0); layout(1)

drop1(m0, scope=.~., test="F")

## Exista desbalanceamento. Será que todos os efeitos foram estimados?
sum(is.na(coef(m0)))==0

##-------------------------------------------
## Estudo da interação. Contrastes entre níveis de tratamento separado
## por níveis de célula.

library(doBy)

## Médias ajustadas para o efeito de sujeito.
lsm <- LSmeans(m0, effect=c("tratamento", "celula"))
lsm

L <- by(data=lsm$K, INDICES=lsm$grid$celula, FUN=as.matrix)
str(L)

## Contrastes de Tukey significa todos contra todos (all pairwise). Com
## 4 níveis são 6 constrastes 2 a 2.

x <- L[[1]]
str(x)

## Veja outra versão em:
## http://git.leg.ufpr.br/leg/legTools/blob/master/R/apcMatrix.R
apcMatrix <- function(x){
    ## --- All Pairwise Contrasts Matrix ---
    ## x: é uma matriz onde cada linha assumesse ter os coeficientes
    ## para estimar uma média (ou qualquer outra função linear dos
    ## parâmetros). Ela returna a matriz com as funções lineares, em
    ## cada linha, representando cada contraste.
    cbn <- combn(x=nrow(x), m=2)
    K <- apply(cbn, MARGIN=2, FUN=function(l){
        x[l[1], ]-x[l[2], ]
    })
    return(t(K))
}

K <- lapply(L, FUN=apcMatrix)

## Abordagem A: Testes separados por celula. Cada um deles tem cobertura
## nominal das 6 hipóteses que é de 95% (default).
lapply(K,
       FUN=function(l){
           summary(glht(m0, linfct=l), test=adjusted(type="fdr"))
       })

KK <- do.call(rbind, K)
nrow(KK)

## Abordagem B: testes todos juntos. Nivel global geral de 95%.
summary(glht(m0, linfct=KK), test=adjusted(type="fdr"))

## Para ficar mais fácil de identificar os contrastes, é necessário
## atribuir rownames.

O seu CMR foi fundamental para aplicar uma solução para os seus dados.
Porém, conhecendo os seus recursos, você poderia ter hospedado os dados na
sua página pessoal (diretório na nuvem), pasta pública do dropbox ou
repositório no GithHub. Com isso incluiria apenas o link para leitura dos
dados dentro do read.table() poupando espaço na mensagem. Mensagens com
corpo grande são descartadas pelo serviço de lista. Me espantou que a sua
passou, então procure sempre evitar.

À disposição.
Walmes.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151215/86656663/attachment.html>


Mais detalhes sobre a lista de discussão R-br