Arquivo

Posts Tagged ‘dbc’

Perda de observações em experimento em DBC

Muitas pessoas não sabem o que fazer quando tem observações não disponíveis no seu experimento – o termo observações perdidas dá a sensação de que um dia elas existiram mas por descuido se perderam. Outras pessoas sentem-se seguras porque têm a licença de determinado aplicativo para fazer análise que, segundo a crença popular, resolve tudo. Bem, o que se deve ter claro em mente, em qualquer condição experimental é: quais as hipóteses de interesse e quais os pressupostos.
Um experimento em delineamento de blocos ao acaso (os blocos não são sorteados, mas os níveis dos demais fatores que são sorteados em cada bloco) com um único fator, normalmente chamado fator tratamento, é representado pelo seguinte modelo estatístico

Y_{ij} = \mu + \beta_i + \tau_j + \epsilon_{ij} ,

em que \mu é uma constante incidente em todas as observações, \beta representa o efeito dos níveis do fator bloco, \tau representa o efeito dos níveis do fator tratamento e \epsilon é um termo que representa o erro. Os pressupostos desse modelo são:

  • os efeitos de bloco e tratamento são aditivos;
  • o erro é independente, de mesma variância e com distribuição normal (condição necessária para condução de inferências);

O fator bloco representa alguma fonte de variação identificável nas unidades experimentais (UE), como por exemplo idade, tamanho, vigor, topografia, umidade, procedência, posição, histórico, comportamento, sexo. Caso não exista uma fonte de variação identificável não há porque, nem como classificar/separar, as UE em blocos que sejam diferentes entre si. Há casos em que as UE são iguais mas ao longo do período experimental podem sofrer influências identificáveis devido, por exemplo, às pessoas que operam, as posições das UI no campo, às épocas de colheita, aos avaliadores. As diferenças entre as UE de um mesmo bloco não podem ser causadas por uma fonte identificável, e sim serem devidas somente ao acaso. Em cada nível de bloco, os níveis de tratamento são sorteados às UE. O experimento é conduzido e as variáveis resposta são observadas.

Às vezes, ocorre de algumas variáveis não serem medidas em algumas UE. Se os valores não disponíveis (NA, not available) surgiram ao acaso, a análise não é inviabilizada. Há casos em que NA depende do tratamento, por exemplo, vacas invadem a área experimental e comem as UI que são do melhor capim, uma doença não controlada mata as UE do genótipo não resistente, pessoas desavisadas colhem os frutos da variedade precoce do pomar experimental. Nesse caso, os NA não sugiram ao acaso e a análise não será confiável.

Nesse tipo de experimento, não interesse em avaliar as diferenças entre níveis de bloco. Já espera-se haver diferenças mas não se quer fazer inferências sobre esses valores. Os blocos são (e devem ser) considerados na análise porque são uma fonte identificável de variação. Por outro lado, para os níveis do fator tratamento se estabelece um conjunto de hipóteses de interesse. São comuns comparações entre os valores esperados (\mu+\tau ) dos níveis contra o nível de referência (testemunha/controle) e comparações duas a duas. Então é necessário que se estime os parâmetros para testar essas hipóteses. No caso balanceado, a média amostral para os níveis de tratamento é um ótimo estimador para (\mu+\tau ) mas caso desbalanceado, a média amostral é um estimador viesado.

O procedimento que vou apresentar aqui é matricial. É o procedimento geral para estimação das médias populacionais marginais (ou valor esperado, \mu+\tau ). As médias amostrais coincidem com esse estimador matricial no caso balanceado. Esse estimador é chamado de estimador de mínimos quadrados. No final apresento um ajuste considerando os efeitos de bloco aleatório mas que vou detalhar as propriedades num post futuro. O código abaixo está comentado para orientação do leitor. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados artificiais gerados a partir do modelo de 4 blocos e 3 tratatamentos

da <- expand.grid(bl=gl(4,1,la="b"), tr=gl(3,1,la="t"))

# vetor de parâmetros (ou efeitos)
theta <- c(mu=10, b1=-1, b2=-1, b3=0, b4=2, t1=-1, t2=0, t3=1)

# soma dos efeitos de bl e tr dá zero, isso é apenas para facilitar
sum(theta[2:5]); sum(theta[6:8])

# médias populacionais marginais dos níveis do fator tratamento
# obtidas fazendo-se a média ao longo dos níveis do fator bloco
#          mu b1   b2   b3   b4   t1 t2 t3
rbind(t1=c(1, 1/4, 1/4, 1/4, 1/4, 1, 0, 0),
      t2=c(1, 1/4, 1/4, 1/4, 1/4, 0, 1, 0),
      t3=c(1, 1/4, 1/4, 1/4, 1/4, 0, 0, 1))%*%theta

# matrizes de incidência dos efeitos dos termos do modelo
X <- mapply(model.matrix, object=list(mu=~1, bl=~-1+bl, tr=~-1+tr),
            MoreArgs=list(data=da))

# as matrizes são irrestritas, ou seja, não são omitidas colunas
X <- do.call(cbind, X)

# valores esperados para cada observação do experimento
da$eta <- X%*%theta

# valor observado para cada observação do experimento
# uma variância pequena foi usada para compararmos os resultados
da$y <- da$eta+rnorm(nrow(da), 0, 0.01)

with(da, tapply(eta, list(bl, tr), function(x) x)) # esperados
with(da, tapply(y, list(bl, tr), function(x) x))   # observados

#-----------------------------------------------------------------------------
# estimação dos efeitos do modelo

m0 <- lm(y~bl+tr, data=da)
coef(m0) # estimativas dos efeitos, apenas k-1 de cada fator são estimados
with(da, tapply(y, tr, mean)) # médias amostrais marginais para tratamentos

# obtidas fazendo-se a média ao longo dos níveis do fator bloco
# estamos usando as estimativas e k-1 efeitos para cada fator
#          mu b2   b3   b4   t2 t3
rbind(t1=c(1, 1/4, 1/4, 1/4, 0, 0),
      t2=c(1, 1/4, 1/4, 1/4, 1, 0),
      t3=c(1, 1/4, 1/4, 1/4, 0, 1))%*%coef(m0)

require(gmodels)

# estimable() é uma forma diferente de fazer essa conta
#                         mu b2   b3   b4   t2 t3
estimable(m0, rbind(tr1=c(1, 1/4, 1/4, 1/4, 0, 0),
                    tr2=c(1, 1/4, 1/4, 1/4, 1, 0),
                    tr3=c(1, 1/4, 1/4, 1/4, 0, 1)))

# deve se ter habilidade para manipilar essas funções lineres
# elas mudam com o tipo de contraste empregado

require(contrast)
# contrast() é outra forma de fazer a conta que não exige conhecimento
# da parametrização usada, type="average" é para fazer média nos blocos
contrast(m0, type="average", list(bl=levels(da$bl), tr="t1"))
contrast(m0, type="average", list(bl=levels(da$bl), tr="t2"))
contrast(m0, type="average", list(bl=levels(da$bl), tr="t3"))

c0 <- contrast(m0, type="average", list(bl=levels(da$bl), tr="t3"))
c0$X # vetor de coeficientes, foi o que passamos para a estimable()

# independência entre efeitos de blocos e tratamentos
round(vcov(m0)*100,3)

#-----------------------------------------------------------------------------
# teste de Tukey para contraste entre médias (populacionais marginais)

require(agricolae)

# esse teste é baseado nas médias amostrais
with(da, HSD.test(y, tr, df.residual(m0), deviance(m0)/df.residual(m0)))

#-----------------------------------------------------------------------------
# o que fazer quando perdemos uma observação? entrar em pânico?
# todos os procedimento matriciais acima continuam válidos

da$z <- da$y
da$z[1] <- NA
with(da, tapply(z, list(bl, tr), function(x) x)) # 1 valor perdido

m1 <- lm(z~bl+tr, data=da)
coef(m1) # os efeitos ainda são estimáveis
round(vcov(m1)*100,6) # perde a ortogonalidade dos efeitos
with(da, tapply(z, tr, mean, na.rm=TRUE)) # médias amostrais são viesadas

# médias de mínimos quadrados são estimadores não viesados
# o desbalanceamento só causou diferença de precisão nas médias
#                         mu b2   b3   b4   t2 t3
estimable(m1, rbind(tr1=c(1, 1/4, 1/4, 1/4, 0, 0),
                    tr2=c(1, 1/4, 1/4, 1/4, 1, 0),
                    tr3=c(1, 1/4, 1/4, 1/4, 0, 1)))

contrast(m1, type="average", list(bl=levels(da$bl), tr="t1"))
contrast(m1, type="average", list(bl=levels(da$bl), tr="t2"))
contrast(m1, type="average", list(bl=levels(da$bl), tr="t3"))

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

# esse teste é baseado nas médias amostrais que agora são viesadas
# média harmônica o número de repetições e usado
with(da, HSD.test(z, tr, df.residual(m0), deviance(m0)/df.residual(m0)))

# viés das médias e aproximação do número de repetições -> teste não confiável

#-----------------------------------------------------------------------------
# contraste entre médias populacionais marginais

contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t1"),
         list(bl=levels(da$bl), tr="t2"))
contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t1"),
         list(bl=levels(da$bl), tr="t3"))
contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t2"),
         list(bl=levels(da$bl), tr="t3"))

# essa sequência de testes de hipótese tem controle do nível nomínal de
# significância por teste, mas não controle global sobre todos os testes
# para isso deve se adotar o procedimento de Bonferroni ou outro
# não requer usar type="average" porque os efeitos de bloco se
# cancelam no contraste

#-----------------------------------------------------------------------------
# multcomp tem procedimentos que controlar o nível global de significância
# para o conjunto de hipóteses de interesse

require(multcomp)
summary(glht(m1, linfct=mcp(tr="Tukey")))

# Tukey não quer dizer teste de Tukey, e sim contrastes de Tukey, ou seja,
# todos as comparações duas a duas entre as médias

#-----------------------------------------------------------------------------
# níveis do fator bloco tem efeito aleatório

require(nlme)
m2 <- lme(z~tr, random=~1|bl, data=da, na.action=na.omit)

# não tem diferença de precisão nas estimativas das médias
estimable(m2, rbind(tr1=c(1,0,0), tr2=c(1,1,0), tr3=c(1,0,1)))
contrast(m2, list(tr=levels(da$tr)))

# contrastes devem ser feitos com as médias marginais sempre
# no caso de balanceamento, as médias amostrais são estimadores não viesados
# das médias populacionais marginais

#-----------------------------------------------------------------------------
Anúncios

Funções para análise de experimentos balanceados

ExpDes é um pacote não oficial do R desenvolvido por Denismar Alves Nogueira, Eric Batista Ferreira e Portya Piscitelli Cavalcanti. O pacote contém funções para análise de dados de experimentos balanceados. As funções são de fácil aplicação e apresentam como resultado o quadro de análise de variância, teste de normalidade para os resíduos, o teste de médias para comparar níveis dos fatores (quando qualitativos) e ajuste de modelos de regressão polinomial (quando quantitativo).

As funções disponíveis no pacote ExpDes (nomes em português) podem ser aplicadas aos seguintes tipos de delineamentos experimentais, sempre considerando fatores de efeito fixo:

função descrição
dic() Delineamento Inteiramente Casualizado balanceado com um só fator.
dbc() Delineamento em Blocos Casualizados balanceado com um só fator.
dql() Delineamento em Quadrado Latino balanceado com um só fator.
fat2.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo.
fat2.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo.
fat2.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo com um tratamento adicional.
fat2.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo com um tratamento adicional.
fat3.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo.
fat3.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo.
fat3.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo com um tratamento adicional.
fat3.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo com um tratamento adicional.
psub2.dic() Delineamento Inteiramente Casualizado balanceado em esquema de parcelas subdivididas.
psub2.dbc() Delineamento em Blocos Casualizados balanceado em esquema de parcelas subdivididas.

Para casos experimentais onde houve perda de observações (desbalanceados), experimentos com blocos incompletos, experimentos em parcela subsubvividida, experimentos com mais de um tratamento adicional, experimentos com emprego de regressão não linear, experimentos com fatores de efeito aleatório, experimentos com emprego de confundimento entre fatores (fatoriais fracionários), experimentos com resposta não normal (e.g. proporções, contagens) pode-se usar as funções lm(), aov(), nls(), e nlme::lme(), glm() do R, dentre outras. No entanto, o usuário deve ter um conhecimento mais aprofundado sobre o tipo de experimento em questão para evitar conclusões equivocadas dos seus resultados. Esses tipos de experimentos exigem mais conhecimento de modelos lineares, álgebra de matrizes, estimação e inferência.

O pacote está disponível nas versões português e inglês. Confesso que só consegui instalar a versão em inglês mas sei que pessoas que instalaram a versão português com sucesso. Sempre ao usar um pacote, faça a citação do mesmo. Use a função citation("ExpDes"). No link está disponível o manual em português do pacote.