Início > delineamento > Perda de observações em experimento em 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

#-----------------------------------------------------------------------------
About these ads
  1. Nenhum comentário ainda.
  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 51 outros seguidores

%d blogueiros gostam disto: