source("config/_setup.R")
Nessa vignette vamos considerar BanzattoQd3.2.1
, já disponível no
pacote labestData.
library(labestData) ls("package:labestData")
help(BanzattoQd3.2.1, help_type = "html")
Este conjunto de dados é de um experimento instalado em delineamento
inteiramente casualizado (DIC) para estudar diferentes produtos para
controle de pulgão na cultura do pepino. Foram avaliados 4 produtos
(trat
), e mais uma testemunha, para controle do pulgão
(Aphis gossypii Glover). A variável resposta desse experimento foi o
número de pulgões encontrados, uma variável resposta do tipo contagem.
Em caso de não ter o pacote labestData, você pode ler a tabela de dados diretamente do repositório. No entanto, recomendamos que instale o pacote para ter acesso a todos os datasets e às documentações
Antes de iniciar a análise dos dados do experimento, é muito recomendável fazer uma análise exploratória que serve basicamente para
#----------------------------------------------------------------------- # Ler a partir do repositório do labestData. # url <- paste0("https://gitlab.c3sl.ufpr.br/pet-estatistica", # "/labestData/raw/devel/data-raw/BanzattoQd3.2.1.txt") # # BanzattoQd3.2.1 <- # read.table(file = url, sep = "\t", header = TRUE, # colClasses = c("factor", "integer", "integer")) #----------------------------------------------------------------------- # Análise exploratória. # Estrutura do objeto. str(BanzattoQd3.2.1) # Tabela de frequência para os tratamentos. xtabs(~trat, data = BanzattoQd3.2.1) # Dados desemplilhados. t(unstack(x = BanzattoQd3.2.1, form = pulgoes ~ trat)) # Média e desvio-padrão das observações em cada nível. aggregate(pulgoes ~ trat, data = BanzattoQd3.2.1, FUN = function(y) { c(mean = mean(y), dp = sd(y)) }) library(lattice) packageVersion("lattice") # Diagrama de dispersão. xyplot(pulgoes ~ reorder(trat, pulgoes), data = BanzattoQd3.2.1, xlab = "Produtos para controle de pulgão", ylab = "Número de pulgões 36hs após pulverização", scales = list(x = list(rot = 90)))
Pela análise exploratória temos que esse experimento é balanceado pois tem número igual de repetições de cada produto. O diagrama de dispersão mostra os valores observados do número de pulgões em função dos produtos. Os níveis foram ordenados de acordo com os números médios de pulgões. É imediato perceber a partir do gráfico que a variabilidade do número de pulgões aumenta com a média. Isso é importante porque adianta um possível problema com o pressuposto de homogeneidade de variâncias do modelo.
O modelo estatístico correspondente ao delineamento é
$$ \begin{aligned} Y|\text{produtos} &\,\sim \text{Normal}(\mu_i, \sigma^2) \newline \mu_i &= \mu+\tau_i \end{aligned} $$
em que (\tau_i) é o efeito do produto (i) sobre a variável resposta (Y), (\mu) é a média de (Y) na ausência do efeito dos produtos, (\mu_i) é a média para cada nível de produto e (\sigma^2) é a variância das observações ao redor desse valor médio (ou variância residual).
Para fazer a estimação desses parâmetros é necessário que seja considerada alguma restrição paramétrica. Para isso, no R por default considera-se o efeito do primeiro nível como zero ($\tau_1 = 0$). É possível mudar isso nas opções de contrastes do R. Outra parametrização, bem conhecida e usada, é a retrição do tipo soma zero dos efeitos ($\sum \tau_i = 0$) e, para essa situação, $\mu$ representa a média geral do experimento.
Existem duas funções no R para ajustar esse modelo: lm()
e
aov()
. Embora muitos acreditem que aov()
deva ser preferida, por
causa do nome AOV (Analysis Of Variance), a função lm()
é tão
apropriada quanto. Inclusive, objetos vindos da lm()
dispõem de mais
métodos que a da aov()
em pacotes para fazer comparações múltiplas, por
exemplo. Sendo assim, vamos usar a lm()
.
#----------------------------------------------------------------------- # Ajuste do modelo. m0 <- lm(pulgoes ~ trat, data = BanzattoQd3.2.1) # Estimativas dos efeitos. Restrição de zerar primeiro nível. cbind(coef(m0)) # Matriz de contrastes sob a retrição zerar primeiro nível. K <- cbind(1, contrasts(BanzattoQd3.2.1$trat)) K # Médias estimadas pelo modelo. K %*% coef(m0) # Médias amostrais. aggregate(pulgoes ~ trat, data = BanzattoQd3.2.1, FUN = mean)
Conforme comentado, devido à restrição usada, os coeficientes estimados
não são as médias em cada nível. Para a restrição do primeiro nível ter
efeito nulo, o coeficiente (Intercept)
é a média do primeiro nível, no
caso r levels(BanzattoQd3.2.1$trat)[1]
. Os demais coeficientes
são estimativas dos constrastes entre médias de cada um dos níveis
restantes com relação ao primeiro nível.
Como nesse exemplo existe uma testemunha, seria conveniente que a
testemunha fosse o primeiro nível para então termos os contrastes como
resultado do ajuste do modelo. No entanto, como a testemunha é o
último nível, pode-se mudar o tipo de contraste para contr.SAS
, cuja
restrição é zerar o efeito do último nível. Na sequência temos os
resultados das duas opções.
m1 <- update(m0, contrasts = list(trat = "contr.SAS")) coef(m1) BanzattoQd3.2.1 <- transform(BanzattoQd3.2.1, trat = relevel(trat, ref = "Testemunha")) levels(BanzattoQd3.2.1$trat) m0 <- update(m0, data = BanzattoQd3.2.1) coef(m0)
Para fazer inferência a partir do modelo ajustado aos dados, é necessário verificar se os pressupostos do modelo foram atendidos. Caso não sejam, as inferências produzidas não são fiéis, pois a distribuição amostral das estatísticas de teste não são as mesmas sob violação dos pressupostos.
#----------------------------------------------------------------------- # Exibe o quarteto fantástico da avaliação dos pressupostos. par(mfrow = c(2, 2)) plot(m0); layout(1)
O conjunto de gráficos acima exibe fuga dos pressupostos. A suposição de homogeneidade de variâncias é claramente violada. Vemos o padrão cone no gráfico 1-1 dessa matriz (Residuals vs Fitted). No gráfico 2-1 (Scale-Location), temos ênfase desse padrão porque a linha média mostra uma relação positiva com a dispersão, representada pela raíz dos valores absolutos dos resíduos padronizados, versus os valores ajustados (média). Ou seja, a relação média-variância não é nula, pois se fosse, a linha de tendência seria horizontal. Por fim, o gráfico 1-2 (Normal Q-Q), exibe fuga do pressuposto de normalidade dos erros, influenciada pela relação média-variância.
Da forma como está, o modelo não tem os pressupostos atendidos. Logo não é recomendado fazer inferências. Precisamos remediar a fuga dos pressupostos antes e a opção mais comum é procurar uma tranformação que estabilize a variância. Em casos como esse, as transformações do tipo potência são capazes de estabilizar a variância. Vamos considerar a transformação Box-Cox para isso.
A tranformação Box-Cox baseia-se na escolha do valor $\lambda$ de tal maneira que se aumente a compatibilidade dos dados com a suposição de normalidade. A transformação aplicada aos dados é
$$ \begin{aligned} y_{\text{transformado}} &= \frac{y^{\lambda}-1}{\lambda}, \quad \text{ se } \lambda\neq 0 \newline y_{\text{transformado}} &= \log(y), \quad \text{ se } \lambda=0. \end{aligned} $$
A transformação só pode ser aplicada para respostas positivas, (y>0) . Caso tenha-se zeros ou valores negativos podemos somar uma constante aos dados. Como o que interessa da transformação é parte não linear (a potenciação), na prática apenas se aplica a função potência sem subtrair 1 e dividir por $\lambda$. No entanto, essa é a transformação de fato aplicada para calcular o valor da log-verossimilhança. Se houver a curiosidade de fazer o cálculo da log-verrossimilhança passo a passo, lembre-se de considerar o jacobiano da transformação.
MASS::boxcox(m0) abline(v = c(0, 1/3), lty = 2, col = 2)
A função perfil de log-verossilhança em $\lambda$ tem intervalo de confiança (95%) que inclui os valores 0 e 1/3, correspondentes às transformações log e raíz cúbica, respectivamente. Não há necessidade de usar o valor exato correspondente ao máximo da função porque isso não vai mudar a normalidade da variável transformada e, como o valor exato é um valor quebrado, valores de $\lambda$ que correspondem a transformações conhecidas são convenientes, como log, raíz quadrada e raíz cúbica.
Encontrada a transformação a ser aplicada, devemos ajustar o modelo com a variável transformada e avaliar novamente os pressupostos. Não existe garantia de qua a transformação irá corrigir a fuga. O que ela faz é remediar o problema, mas o resultado pode ainda não ser satisfatório.
# Atualiza o modelo anterior passando o log na resposta. m0 <- update(m0, log(.) ~ .) par(mfrow = c(2, 2)) plot(m0); layout(1)
Note agora que os gráficos não exibem padrões de fuga de pressupostos. Pelo gráfico 2-1, a relação média variância agora está estável (ou nula). Não há fugas de normalidade de acordo com o gráfico 1-2. Portanto, foi bem sucedida a aplicação da transformação log ao número de pulgões para que o modelo tivesse os pressupostos atendidos. Agora é possível proceder com as inferências.
O objetivo de ter feito o experimento é saber se existe efeito dos produtos para o controle do pulgão na cultura do pepino. Ou seja, se o número médio de pulgões muda conforme o produto aplicado.
O efeito dos produtos é representado pelos coeficientes $\tau_i$ no modelo estatístico do experimento. Se não existir efeito dos produtos, os valores estimados para $\tau_{i}$ serão próximos a zero. A hipóteses nula de não haver efeito dos produtos é representada por
$$ \text{H}_{0}: \tau_i = 0, \text{para todo}\,i. $$
Essa hipótese é avaliada pelo estatística F do quadro de análise de variância.
anova(m0)
Pelo quadro de anova (analysis of variance), temos que o valor da estatística F produziu um $p$-valor bem pequeno, menor que qualquer nível de significancia adotado na prática. Dessa maneira, rejeitamos a hipótese nula de não haver efeito dos produtos no log do número de pulgões ou, por consequência, no número de pulgões.
Uma vez que existe efeito dos produtos, faz-se necessário aprofundar o conhecimento sobre eles. As hipóteses a serem testadas devem ser estabalecidas antes da coleta dos dados e ajuste do modelo. Como foram estudados 4 níveis de produto e uma testemunha, as hipóteses vão agora considerar diferenças entre cada produto e a testemunha. Como a testemunha é o primeiro nível, temos que testar as hipóteses individuais:
$$ \text{H}_{0}: \mu_i - \mu_1 = 0, \text{para cada}\,i \neq 1. $$
Da forma como ajustamos o modelo, temos as estimativas dos contrastes
disponíveis no vetor de coeficientes. O teste individual para cada
parâmetro é exibido com o summary()
do modelo.
summary(m0)
k <- nlevels(BanzattoQd3.2.1$trat) - 1 p <- 1 - (1 - 0.05)^k
Na tabela dos coeficientes tem-se o valor t correspondente a cada contraste.
O (Intercept)
é a média do log do número de pulgões para a
testemunha. Os demais representam contrastes dos produtos com relação à
testemunha. Baseado nos $p$-valores, todas as hipóteses de contraste nulo
com a testemunha são rejeitadas. Ou seja, todos os produtos apresentam número
médio diferente da testemunha, no caso um menor valor para média
se comparado à testemunha.
Deve-se mencionar que se o nível de significância adotado é de 5% para
cada teste. Assim, a probabilidade de rejeitar pelo menos uma hipótese é $1 -
(1 - 0,05)^{r k
} = r p
$, ou seja, superior a 5%. Com o aumento
do número de hipóteses testadas, é consequência aumentar o nível de
significância global (probabilidade de rejeitar pelo menos uma
hipótese). Existem diferentes formas de corrigir a elevação do nível
global de significancia, mas esse não é um objetivo desta vignette.
Uma adequada representação dos resultados é um gráfico com os valores médios e intervalo de confiança para os níveis dos produtos.
pred <- data.frame(trat = levels(BanzattoQd3.2.1$trat)) pred <- cbind(pred, as.data.frame(predict(m0, newdata = pred, interval = "confidence"))) pred$trat <- reorder(pred$trat, pred$fit) head(pred) library(latticeExtra) segplot(trat ~ lwr + upr, centers = fit, data = pred, draw = FALSE, horizontal = FALSE, xlab = "Produtos para controle de pulgão", ylab = "Log do número de pulgões", scales = list(x = list(rot = 90)), panel = function(x, y, z, centers, ...) { panel.segplot(x = x, y = y, z = z, centers = centers, ...) panel.text(x = as.numeric(z), y = centers, label = sprintf("%0.2f", centers), srt = 90, cex = 0.8, adj = c(0.5, -0.5)) })
É simples produzir o plano de um delineamento inteiramente
casualizado. Vamos gerar o plano para executar um experimento
semelhante ao apresentado, mas com 8 repetições e mais um produto,
o Novo
produto.
trat <- rep(x = c(levels(BanzattoQd3.2.1$trat), "Novo"), times = 8) # Opção 1: sorteie os níveis para as unidades experimentais ordenadas. data.frame(trat = sample(trat), ue = 1:length(trat)) # Opção 2: sorteie as unidades experimentais ordenadas para os níveis. data.frame(trat = trat, ue = sample(1:length(trat)))
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.