2.1 Introdução


Em um delineamento inteiramente casualizado com um fator de tratamento, \(n\) unidades experimentais são divididas aleatoriamente em \(t\) grupos. Cada grupo está então sujeito a um dos níveis ou valores exclusivos do fator de tratamento. Se \(n = tr\) for um múltiplo de \(t\), então cada nível do fator será aplicado a \(r\) unidades experimentais únicas e haverá \(r\) réplicas de cada execução com o mesmo nível do fator de tratamento. Se \(n\) não for múltiplo de \(t\), então haverá um número desigual de réplicas de cada nível de fator. Todas as outras variáveis independentes conhecidas são mantidas constantes para não distorcerem os efeitos. Este delineamento deve ser utilizado quando há apenas um fator em estudo e as unidades experimentais são homogêneas.

Por exemplo, numa experiência para determinar o efeito do tempo de crescimento na altura da massa de pão, um lote homogéneo de massa de pão seria dividido em \(n\) formas de pão com uma quantidade igual de massa em cada uma. As formas de massa seriam então divididas aleatoriamente em \(t\) grupos. Cada grupo poderia crescer por um tempo único e a altura da massa crescida seria medida e registrada para cada pão.

O fator de tratamento seria o tempo de subida, a unidade experimental seria um pão individual e a resposta seria a altura medida. Embora se saiba que outros fatores, como a temperatura, afetam a altura da massa de pão levedada, eles seriam mantidos constantes e cada pão poderia crescer sob as mesmas condições, exceto pelos diferentes tempos de crescimento.


2.2 Replicação e randomização


A replicação e a randomização foram popularizadas por R.A. Fisher. Estas são as primeiras técnicas que se enquadram na categoria de controle de erros explicada na Seção 1.4.

A técnica de replicação determina que \(r\) pães sejam testados em cada um dos \(t\) tempos de crescimento, em vez de um único pão em cada tempo de crescimento. Ao ter unidades experimentais replicadas em cada nível do fator de tratamento, a variância do erro experimental pode ser calculada a partir dos dados e esta variância será comparada com os efeitos do tratamento.

Se a variabilidade entre as médias dos tratamentos não for maior que a variância do erro experimental, as diferenças nos tratamentos provavelmente se devem a diferenças nas unidades experimentais atribuídas a cada tratamento. Sem replicação é impossível dizer se as diferenças de tratamento são reais ou apenas uma manifestação aleatória das unidades experimentais específicas utilizadas no estudo. Subamostras ou medições duplicadas, descritas no Capítulo 1, não podem substituir réplicas.

A divisão aleatória de unidades experimentais em grupos é chamada de randomização e é o procedimento pelo qual a validade do experimento é garantida contra vieses causados por outras variáveis ocultas. No experimento de crescimento do pão, a randomização evitaria que variáveis ocultas, como a variabilidade no fermento de pão para pão e tendências na técnica de medição ao longo do tempo, distorcessem o efeito do tempo de crescimento.

Quando unidades experimentais são randomizadas para níveis de fatores de tratamento um teste exato da hipótese de que o efeito do tratamento é zero pode ser realizado utilizando um teste de randomização e um teste de parâmetros no modelo linear geral normalmente utilizado na análise de dados experimentais é uma boa aproximação para o teste de randomização.

Uma maneira simples de construir um formulário de coleta de dados aleatório, dividindo \(n\) unidades experimentais em \(t\) grupos de tratamento, pode ser realizada usando comandos base R. Por exemplo, no experimento de crescimento de pão, se o pesquisador quiser examinar três tempos de crescimento diferentes: 35 minutos, 40 minutos e 45 minutos e testar quatro réplicas de pães em cada tempo de crescimento, o código a seguir criará a lista .

set.seed(7638)
f = factor( rep( c(35, 40, 45 ), each = 4))
fac = sample( f, 12 )
eu = 1:12
plan = data.frame( loaf=eu, time=fac ) # loaf = pão
head(plan)
##   loaf time
## 1    1   40
## 2    2   40
## 3    3   35
## 4    4   45
## 5    5   35
## 6    6   45


O comando R factor cria um vetor dos níveis do fator para o tempo de subida e o armazena na variável f. Há também um comando ordered em R que cria um fator que se presume ter níveis numéricos igualmente espaçados. R lida com fatores criados pelos comandos factor e ordered de maneira diferente ao fazer comparações de tratamentos após ajustar um modelo. Haverá mais discussão sobre isso na Seção 2.8.

A função sample randomiza a ordem dos níveis dos fatores e armazena o vetor randomizado na variável fac. A função seq cria um vetor numérico de unidades experimentais, ou seja, pão guardados no objeto R eu. Em seguida, a função data.frame combina os dois vetores eu, fac como colunas que são armazenadas em plan, um objeto do quadro ou conjunto de dados com títulos de coluna loaf (pão) e time (tempo). Finalmente, a função head mostra as primeiras 6 linhas no quadro de dados.

Estes dados nos mostra que o primeiro pão ou unidade experimental, deve crescer 40 minutos, o segundo pão 45 minutos, etc. Se você executar os mesmos comandos em R repetidamente, obterá a mesma ordem aleatória por causa do comando set.seed. Remova esta instrução para obter uma ordem aleatória diferente.

Além dos comandos básicos do R mostrados acima, vários pacotes R podem criar listas aleatórias de experimentos, que podem ser convenientemente convertidas em formulários eletrônicos de coleta de dados. Entretanto, esses pacotes serão ilustrados para a criação de projetos mais complicados nos próximos capítulos e não serão mostrados agora.


2.3 Um exemplo histórico


Para ilustrar a lista de verificação para planear uma experiência descrita na Secção 1.6, considere um exemplo histórico retirado do Relatório da Estação Experimental de Rothamstead de 1937 (desconhecido, 1937). Isto ilustra alguns dos primeiros trabalhos realizados por R.A. Fisher no desenvolvimento de ideias de desenho experimental e análise de variância para uso em experimentos agrícolas na estação de pesquisa.

Objetivos

O objetivo do estudo foi comparar as épocas de plantio e os métodos de aplicação de fertilizantes artificiais mistos (NPK) antes do plantio, no rendimento da beterraba sacarina. Normalmente o fertilizante é aplicado e as sementes plantadas assim que o solo pode ser trabalhado.

Unidades experimentais

As unidades experimentais foram as parcelas de solo combinadas com sementes específicas a serem plantadas em cada parcela de solo.

Resposta ou variável dependente

A variável dependente seria a rendimento de beterraba sacarina medido em cwt por acre.

Variáveis independentes e variáveis ocultas

As variáveis independentes de interesse foram o tempo e o método de aplicação de fertilizantes artificiais mistos. Quatro níveis do fator de tratamento foram escolhidos conforme listado abaixo:

    1. nenhum fertilizante artificial aplicado
    1. artificiais aplicados em janeiro (arados)
    1. artigos aplicados em janeiro (transmissão)
    1. artigos aplicados em abril (transmissão)

As variáveis ocultas que poderiam causar diferenças nos rendimentos da beterraba sacarina entre parcelas eram diferenças na fertilidade das próprias parcelas, diferenças nas sementes de beterraba usadas em cada parcela, diferenças entre parcelas no nível de infestação de ervas daninhas, diferenças nas práticas de cultivo de desbaste beterraba e colheita manual da beterraba.

Testes Piloto

A beterraba sacarina foi cultivada rotineiramente em Rothamstead e fertilizantes artificiais foram usados tanto na lavoura quanto na lavoura para muitas plantas cultivadas; portanto, sabia-se que a variável independente poderia ser controlada e que a resposta era mensurável.

Escolha do delineamento experimental

O delineamento completamente aleatório (CRD) foi escolhido de modo que as diferenças nas variáveis ocultas entre as parcelas dificilmente correspondessem a mudanças nos níveis dos fatores listados acima.

Determinar o número de réplicas

Uma diferença no rendimento de 6 cwt por acre foi considerada de importância prática e, com base em estimativas históricas de variabilidade nos rendimentos da beterraba sacarina em Rothamstead, quatro ou cinco réplicas foram determinadas como suficientes.

Randomização de unidades experimentais para níveis de tratamento

Dezoito parcelas foram escolhidas para o experimento e uma lista aleatória foi construída atribuindo quatro ou cinco parcelas para cada nível de fator.


2.4 Modelo linear para o delineamento completamente aleatório


O modelo matemático para os dados de um delineamento completamente aleatório (CRD) ou delineamento completamente casualizado, com um número desigual de repetições para cada nível de fator pode ser escrito como: \[\begin{equation*} \tag{2.1} Y_{ij}=\mu_i+\epsilon_{ij}, \end{equation*}\] onde \(Y_{ij}\) é a resposta para a \(j\)-ésima unidade experimental sujeita ao \(i\)-ésimo nível do fator de tratamento, \(i = 1,\cdots,t\), \(j = 1,\cdots,r_i\) e \(r_i\) é o número de unidades experimentais ou repetições no \(i\)-ésimo nível do fator de tratamento.

Isso às vezes é chamado de modelo de médias celulares com uma média diferente \(\mu_i\), para cada nível do fator de tratamento. A distribuição dos erros experimentais \(\epsilon_{ij}\) são mutuamente independentes devido à randomização e assumidas como normalmente distribuídas.

Este modelo está representado graficamente na Figura 2.2. Figura 2.2 Figura 2.1 Modelo de médias celulares.

Uma maneira alternativa de escrever um modelo para os dados é \[\begin{equation*} \tag{2.2} Y_{ij}=\mu+\tau_i+\epsilon_{ij}\cdot \end{equation*}\] Isso é chamado de modelo de efeitos e os \(\tau_i\)’s são chamados de efeitos. \(\tau_i\) representa a diferença entre a média de longo prazo de todos os experimentos possíveis no \(i\)-ésimo nível do fator de tratamento e a média geral.

Com a suposição de normalidade \(Y_{ij}\sim N(\mu+\tau_i,\sigma^2)\) ou \(\epsilon_{ij}\sim N(0,\sigma^2)\). Para igual número de réplicas, as médias amostrais dos dados no \(i\)-ésimo nível do fator de tratamento são representadas por \[\begin{equation*} \tag{2.3} \overline{y}_{i\bullet}=\dfrac{1}{r_i}\sum_{j=1}^{r_i} y_{ij} \end{equation*}\] e a grande média ou média global é dada por \[\begin{equation*} \tag{2.4} \overline{y}_{\bullet\bullet}=\dfrac{1}{t}\sum_{i=1}^{t} \overline{y}_{i\bullet}=\dfrac{1}{n}\sum_{i=1}^t\sum_{j=1}^{r_i}y_{ij}, \end{equation*}\] onde \(n=sum_i r_i\).

Utilizando o método de máxima verossimilhança, que é equivalente ao método dos mínimos quadrados com essas suposições, as estimativas das médias das células são encontradas escolhendo-as para minimizar o erro da soma dos quadrados \[\begin{equation*} \tag{2.5} SSE = \sum_{i=1}^{t}\sum_{j=1}^{r_i} \big(y_{ij}-\mu_i\big)^2\cdot \end{equation*}\]

Isso é feito tomando derivadas parciais de \(SSE\) em relação à média de cada célula, definindo os resultados iguais a zero e resolvendo cada equação \[ \dfrac{\partial SSE}{\partial \mu_i} = -2\sum_{i=1}^t \sum_{j=1}^{r_i} \big(y_{ij}-\mu_i\big)=0\cdot \] Isso resulta nas estimativas: \[ \widehat{\mu}_i = \overline{y}_{i\bullet} \]


2.4.1 Representação matricial


Considere um CRD com \(t = 3\) níveis de fator e \(r_i = 4\) réplicas para \(i = 1,\cdots,t\). Podemos escrever o modelo de efeitos de forma concisa usando notação matricial como: \[\begin{equation*} \tag{2.6} \pmb{y} = \pmb{X}\beta+\epsilon, \end{equation*}\] onde \[ \pmb{y}=\begin{pmatrix} y_{11} \\ y_{12} \\ y_{13} \\ y_{14} \\ y_{21} \\ y_{22} \\ y_{23} \\ y_{24} \\ y_{31} \\ y_{32} \\ y_{33} \\ y_{34} \end{pmatrix}, \pmb{X}=\begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{pmatrix}, \beta = \begin{pmatrix} \mu \\ \tau_1 \\ \tau_2 \\ \tau_3 \end{pmatrix}, \epsilon = \begin{pmatrix} \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{13} \\ \epsilon_{14} \\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{23} \\ \epsilon_{24} \\ \epsilon_{31} \\ \epsilon_{32} \\ \epsilon_{33} \\ \epsilon_{34} \end{pmatrix}, \] onde \(\epsilon\sim MVN(0,\sigma^2\pmb{I})\).

Os estimadores de mínimos quadrados são a solução das equações normais \(\pmb{X}^\top\pmb{X} = \pmb{X}y\). O problema com as equações normais é que se \(\pmb{X}^\top \pmb{X}\) é singular não pode ser invertido. Usando a codificação de tratamento para um fator não ordenado criado com o comando factor, a função R lm torna a classificação completa da matriz \(\pmb{X}\) eliminando a coluna que corresponde ao primeiro nível do fator, conforme mostrado abaixo: \[ \pmb{X}=\begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{pmatrix}\cdot \]

Esta codificação de tratamento torna o primeiro nível do fator o padrão e todos os outros níveis do fator são comparados a ele. Para o exemplo com \(t = 3\) níveis de fator, a solução para as equações normais é \[ \big(\pmb{X}^\top\pmb{X} \big)^{-1} \pmb{X}^\top y = \widehat{\beta} = \begin{pmatrix} \widehat{\mu}+\widehat{\tau}_1 \\ \widehat{\tau}_2-\widehat{\tau}_1 \\ \widehat{\tau}_3-\widehat{\tau}_1\end{pmatrix}\cdot \]


2.4.2 Cálculos de mínimos quadrados com a função R lm


A Tabela 2.1 mostra os dados de um projeto CRD para o experimento de crescimento de pão descrito anteriormente neste capítulo.

\[ \begin{array}{cc}\hline \mbox{Tempo de subida} & \mbox{Alturas do pão} \\\hline 35 \mbox{ minutos} & 4.5, 5.0, 5.5, 6.75 \\ 40 \mbox{ minutos} & 6.5, 6.5, 10.5, 9.5 \\ 45 \mbox{ minutos} & 9.75, 8.75, 6.5, 8.25 \\\hline \end{array} \] Tabela 2.1: Dados do experimento de crescimento de pão


Usando esses dados temos \[ \pmb{X}^\top \pmb{X}=\begin{pmatrix} 12 & 4 & 4 \\ 4 & 4 & 0 \\ 4 & 0 & 4 \end{pmatrix}, \quad \pmb{X}^\top\pmb{y}=\begin{pmatrix} 88.0 \\ 33.0 \\ 33.25 \end{pmatrix} \] e \[ \big(\pmb{X}^\top \pmb{X}\big)^{-1}=\begin{pmatrix} 0.25 & -0.25 & -0.25 \\ -0.25 & 0.50 & 0.25 \\ -0.25 & 0.25 & 0.50 \end{pmatrix}, \]

\[ \widehat{\beta}=\big(\pmb{X}^\top \pmb{X}\big)^{-1}\pmb{X}^\top\pmb{y}=\begin{pmatrix} \widehat{\mu}+\widehat{\tau}_1 \\ \widehat{\tau}_2-\widehat{\tau}_1 \\ \widehat{\tau}_3-\widehat{\tau}_1\end{pmatrix}= \begin{pmatrix} 5.4375 \\ 2.8125 \\ 2.8750 \end{pmatrix}\cdot \]

A biblioteca ou pacote de dados e comandos R daewr contém este conjunto de dados. A execução do código a seguir abre o conjunto de dados mostrado na Tabela 2.1 e calcula as estimativas.

library(daewr)
data("bread")
head(bread)
##   loaf time height
## 1    1   35   4.50
## 2    2   35   5.00
## 3    3   35   5.50
## 4    4   35   6.75
## 5    5   40   6.50
## 6    6   40   6.50
mod0 = lm( height ~ time, data = bread )
summary( mod0 )
## 
## Call:
## lm(formula = height ~ time, data = bread)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.812 -1.141  0.000  1.266  2.250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4375     0.7655   7.104 5.65e-05 ***
## time40        2.8125     1.0825   2.598   0.0288 *  
## time45        2.8750     1.0825   2.656   0.0262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.531 on 9 degrees of freedom
## Multiple R-squared:  0.5056, Adjusted R-squared:  0.3958 
## F-statistic: 4.602 on 2 and 9 DF,  p-value: 0.042


O comando lm testa o modelo linear e armazena os resultados no objeto mod0, o comando summary imprime os resultados, uma parte dos quais é mostrado utilizando o comando head.

Como a variável time é um fator no quadro de dados bread e a codificação de tratamento padrão foi usada pela função lm, as estimativas descritas acima são produzidas.


2.4.3 Estimação de \(\sigma^2\) e distribuição de formas quadráticas


A estimativa da variância do erro experimental \(\sigma^2\) é \(ssE/(n-t)\). Somente é possível estimar esta variância quando existem experiências replicadas em cada nível do fator de tratamento. Quando medições em subamostras ou medições duplicadas na mesma unidade experimental são tratadas como réplicas, esta estimativa pode ser seriamente distorcida.

Na forma matricial, \(SSE\) pode ser escrita como \[ ssE = \pmb{y}^\top \pmb{y}-\widehat{\beta}^\top \pmb{X}^\top \pmb{y} = \pmb{y}^\top \Big(I-\pmb{X}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{X}^\top \Big)^{-1}\pmb{y}, \] e a partir da teoria dos modelos lineares pode-se mostrar que a razão entre \(SSE\) e a variância do erro experimental \(\sigma^2\), segue uma distribuição qui-quadrado com \(n-t\) graus de liberdade, ou seja, \[ ssE/(n-t) \sim \chi^2_{n-t}\cdot \]


2.4.4 Funções estimáveis


Uma análise de variância padrão ou ANOVA, fornece um teste \(F\), que reflete todas as diferenças possíveis entre as médias dos grupos analisados. No entanto, a maioria dos experimentadores quer tirar conclusões mais precisas do que “a manipulação experimental tem efeito no comportamento dos participantes”.

Conclusões precisas podem ser obtidas a partir da análise de contraste, porque um contraste expressa uma questão específica sobre o padrão de resultados de uma ANOVA. Especificamente, um contraste corresponde a uma previsão suficientemente precisa para ser traduzida num conjunto de números chamados coeficientes de contraste que refletem a previsão.

A correlação entre os coeficientes de contraste e as médias do grupo observado avalia diretamente a similaridade entre a previsão e os resultados. Ao realizar uma análise planejada envolvendo vários contrastes, precisamos avaliar se esses contrastes são mutuamente ortogonais ou não.

Dois contrastes são ortogonais quando seus coeficientes de contraste não estão correlacionados, ou seja, seu coeficiente de correlação é zero. O número de contrastes ortogonais possíveis é um a menos que o número de níveis da variável independente.

Uma combinação linear das médias das células é chamada de função estimável se puder ser expressa como o valor esperado de uma combinação linear das respostas, ou seja, \[\begin{equation*} \tag{2.7} \sum_{i=1}^t b_i(\mu+\tau_i) = \mbox{E}\left( \sum_{i=1}^t \sum_{j=1}^{r_i} a_{ij}Y_{ij}\right)\cdot \end{equation*}\] A partir desta definição, pode-se ver que os efeitos \(\tau_i\), não são estimáveis, mas uma média de célula \(\mu+\tau_i\) ou um contraste de efeitos, \(\sum c_i \tau_i\), onde \(\sum c_i=0\), é estimável.

Em notação matricial \(\pmb{L}\beta\) é um conjunto de funções estimáveis se cada linha de \(\pmb{L}\) for uma combinação linear das linhas de \(\pmb{X}\) e \(\pmb{L}\widehat{\beta}\) é seu estimador não viciado. \(\pmb{L}\widehat{\beta}\) segue a distribuição normal multivariada com matriz de covariância \(\sigma^2 \pmb{L}^\top (\pmb{X}^\top\pmb{X})^{-1}\pmb{L}\) e o estimador da matriz de covariância é \[ \widehat{\sigma}^2 \pmb{L}^\top (\pmb{X}^\top\pmb{X})^{-1}\pmb{L}\cdot \]

Por exemplo, usando os dados do experimento de crescimento de pão acima, \[\begin{equation*} \tag{2.8} \pmb{L}=\begin{pmatrix} 0 & 1 & -1 & 0 \\ 0 & 1 & 0 & -1 \end{pmatrix}, \; \pmb{L}\beta = \begin{pmatrix} \tau_1-\tau_2 \\ \tau_1-\tau_3 \end{pmatrix} \quad \mbox{e} \quad \pmb{L}\widehat{\beta} = \begin{pmatrix} \widehat{\tau}_1-\widehat{\tau}_2 \\ \widehat{\tau}_1-\widehat{\tau}_3 \end{pmatrix}=\begin{pmatrix} -2.8125 \\ -2.8750 \end{pmatrix} \end{equation*}\] é um vetor de contrastes dos efeitos.

Quando uma hipótese de pesquisa é precisa, é possível expressá-la como um contraste. Uma hipótese de pesquisa, em geral, pode ser expressa como uma forma, uma configuração ou uma classificação dos meios experimentais. Em todos estes casos, podemos atribuir números que refletirão os valores previstos das médias experimentais. Esses números são chamados de coeficientes de contraste quando sua média é zero.

O número de graus de liberdade, ou número de contrastes de efeitos linearmente independentes em um delineamento completamente aleatório é sempre o número de níveis do fator de tratamento menos um, ou seja, \(t − 1\). Sempre que houver um conjunto de \(t − 1\) linearmente independentes contrastes dos efeitos, eles são chamados de conjunto saturado de contrastes estimáveis.

Pode ser mostrado que \((\pmb{L}\widehat{\beta})^\top \Big(\pmb{L}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{L}^\top\Big)^{-1}(\pmb{L}\widehat{\beta})\) segue a distribuição qui-quadrado não central \(\chi^2_p(\lambda)\) onde o parâmetro de não centralidade é \[ \lambda = (\sigma^2)^{-1}(\pmb{L}{\beta})^\top \Big(\pmb{L}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{L}^\top\Big)^{-1}(\pmb{L}{\beta}) \] e \(\pmb{L}\) é a matriz de coeficientes para um contraste estimável como (2.8), os graus de liberdade \(p\) são iguais ao posto de \(\pmb{L}\).

Contrastes estimáveis podem ser obtidos na função fit.contrast no pacote R gmodels, Warnes (2012). Esta função é usada para estimar a diferença média nas médias das células para o primeiro e segundo níveis do fator de tratamento, \((\mu+\tau_1)-(\mu+\tau_2)=\tau_1-\tau_2\).

library(gmodels)
fit.contrast( mod0, "time", c(1,-1,0) )
##                   Estimate Std. Error   t value   Pr(>|t|)
## time c=( 1 -1 0 )  -2.8125   1.082532 -2.598076 0.02882905
## attr(,"class")
## [1] "fit_contrast"


Na chamada de função acima, o mod0 é o nome de um modelo anteriormente ajustado com a função R lm, o termo entre aspas é o nome do fator no modelo cujas médias das células são comparadas e o vetor c(1 -1, 0) são os coeficientes de contraste, \(c_i\). Isso produz o resultado -2.8125, que é o negativo da segunda estimativa produzida na saída R da Seção 2.4.2 usando a codificação de tratamento padrão no modelo mod0.


2.4.5 Teste de hipótese de ausência de efeitos do tratamento


No modelo para o delineamento completamente aleatório, a hipótese estatística de interesse é \[ H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_t \qquad \mbox{ou} \qquad \tau_1 = \tau_2 = \cdots = \tau_t \] contra a alternativa \[ H_1 \, : \, \mbox{pelo menos dois dos } \tau's \mbox{ diferem}\cdot \]

Se a hipótese nula for verdadeira, o modelo \(y_{ij}=\mu_i+\epsilon_{ij}=\mu+\tau_i+\epsilon_{ij}\) simplifica para \(y_{ij} = \mu + \epsilon_{ij}\), que pode ser representado como uma única distribuição normal com média \(\mu\) e variância \(\sigma^2\) em vez de múltiplas distribuições normais como as mostradas na Figura 2.2.

A soma dos quadrados sobre a média é \[ ssTotal=\sum_{i=1}^t \sum_{j=1}^{r_i} (y_{ij}-\overline{y}_{\bullet\bullet})^2=\dfrac{\pmb{y}^\top\pmb{y}-(\pmb{1}^\top \pmb{y})^2}{\pmb{1}^\top\pmb{1}}, \] onde \(\overline{y}_{\bullet\bullet}\) é a média geral e \(\pmb{1}\) é um vetor coluna de unidades. Esta soma de quadrados pode ser particionada como: \[\begin{equation*} \tag{2.9} ssTotal = ssT+ssE, \end{equation*}\] onde \[ ssT=\dfrac{\widehat{\beta}^\top\pmb{X}^\top\pmb{y}-(\pmb{1}^\top\pmb{y})^2}{\pmb{1}^\top\pmb{1}}=(\pmb{L}{\beta})^\top \Big(\pmb{L}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{L}^\top\Big)^{-1}(\pmb{L}{\beta}) \] e \(\pmb{L}\) é a matriz de coeficientes para um conjunto saturado de contrastes estimáveis. Essa quantidade é chamada de soma de quadrados de tratamento.

Sob a hipótese nula \(H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_t\), tanto \(ssT\) quanto \(ssE\) seguem a distribuição qui-quadrado. Estas somas de quadrados e os seus correspondentes quadrados médios, que são formados pela divisão de cada soma de quadrados pelos seus graus de liberdade, são apresentados numa análise de variância ou numa tabela ANOVA como a mostrada simbolicamente na Tabela 2.2.

\[ \begin{array}{ccccc}\hline \mbox{Fonte} & \mbox{df} & \mbox{Soma de Quadrados} & \mbox{Quadrados Médios} & F \\ \hline \mbox{Tratamentos} & t-1 & ssT & msT & F=msT/msT \\ \mbox{Erro} & n-t & ssE & msE & \\ \mbox{Total} & n=1 & ssTotal & msTotal & \\ \hline \end{array} \] Tabela 2.2: Tabela de Análise de Variância.

Sob a hipótese nula, a razão \(F\) \(msT/msE\) segue a distribuição \(F\) com \(t-1\) e \(n-t\) graus de liberdade e sob a alternativa segue a distribuição \(F\) não central com parâmetro de não centralidade \[ \lambda = (\sigma^2)^{-1}(\pmb{L}{\beta})^\top \Big(\pmb{L}\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{L}^\top\Big)^{-1}(\pmb{L}{\beta})=\dfrac{r}{\sigma^2}\sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet} \big)^2\cdot \]

É a estatística de teste de razão de verossimilhança generalizada para \(H_0\) e é o método formal de comparação dos efeitos do tratamento com a variância do erro experimental descrita na Seção 2.2.

As somas dos quadrados, quadrados médios, graus de liberdade na tabela ANOVA e estatística do teste \(F\) associada podem ser calculadas pela função aov no R. As entradas para a função aov são as mesmas daquelas para a função lm mostrada anteriormente, mas o resumo de um objeto criado pela função aov é a tabela ANOVA e não as estimativas produzidas pela função lm.

O código para produzir a tabela ANOVA para o experimento de crescimento de massa de pão é:

mod1 = aov( height ~ time, data = bread )
summary(mod1)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## time         2  21.57  10.786   4.602  0.042 *
## Residuals    9  21.09   2.344                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Nesta tabela, \(ssT\) e \(msT\) e os graus de liberdade (df) associados estão na linha denominada time, o \(ssE\) está na linha denominada Residuais e o ssTotal otal pode ser calculado adicionando o ssT ao ssE.

O valor \(F\) é a razão \(msT/msE\) e a última coluna denominada Pr(>F) é a probabilidade de exceder o valor \(F\) calculado se a hipótese nula for verdadeira. Isto é chamado de \(p\)-valor e é ilustrado graficamente na Figura 2.3. Se o experimentador escolher o nível de significância \(\alpha\), para seu teste de hipótese, ele rejeitará a hipótese se o valor de Pr(>F) na saída aov for menor que o valor escolhido de \(\alpha\).

Figure 2.3: Pr > F.

Para o experimento de crescimento de pão existem diferenças significativas entre as alturas médias da massa levedada para cada tempo de levedação no nível de significância \(\alpha = 0.05\), já que 0.042 < 0.05.

Quando vários contrastes são avaliados, vários testes estatísticos são realizados no mesmo conjunto de dados e isso aumenta a probabilidade de um erro Tipo I, ou seja, rejeitar a hipótese nula quando esta for verdadeira. Para controlar o erro Tipo I ao nível do conjunto, também conhecido como família de contrastes, é necessário corrigir o nível \(\alpha_L\) usado para avaliar cada contraste.

Esta correção para múltiplos contrastes pode ser feita usando a equação de \(\check{\mbox{S}}\)id\(\grave{a}\)k, a desigualdade de Bonferroni, também conhecida como Boole ou Dunn ou a técnica de Monte Carlo. A probabilidade de cometer pelo menos um erro Tipo I, para uma família de \(C\) contrastes ortogonais, isto é, \(C\) contrastes estatisticamente independentes é \[ \alpha_L = 1-(1-\alpha)^C \] sendo \(\alpha_L\) o erro Tipo I para a família de contrastes \(L\) e \(\alpha\) sendo o erro Tipo I por contraste, tipicamente os números mágicos 0.01, 0.05 ou 0.10. Esta equação pode ser reescrita como \[ \alpha = 1-(1-\alpha_L)^{1/C}\cdot \]

Esta fórmula, conhecida como equação de \(\check{\mbox{S}}\)id\(\grave{a}\)k mostra como corrigir os valores de \(\alpha\) usados para cada contraste. Como esta equação envolve uma potência fracionária, pode-se usar uma aproximação conhecida como desigualdade de Bonferroni, que relaciona \(\alpha\) a \(\alpha_L\) por \[ \alpha \approx \dfrac{\alpha_L}{C}\cdot \]

\(\check{\mbox{S}}\)id\(\grave{a}\)k e Bonferroni estão relacionados pela desigualdade \[ \alpha=1-(1-\alpha_L)^{1/C}\geq \dfrac{\alpha_L}{C}\cdot \]

Eles são, em geral, muito próximos um do outro. Como pode ser visto, a desigualdade de Bonferroni é uma estimativa pessimista. Conseqüentemente, \(\check{\mbox{S}}\)id\(\grave{a}\)k deve ser preferido. No entanto, a desigualdade de Bonferroni é mais conhecida e, portanto, é usada e citada com mais frequência

A técnica de Monte-Carlo também pode ser usada para corrigir múltiplos contrastes. A técnica de Monte Carlo consiste em realizar diversas vezes um experimento simulado utilizando dados aleatórios, com o objetivo de obter um padrão de resultados que mostre o que aconteceria apenas com base no acaso.

Esta abordagem pode ser usada para quantificar \(\alpha_L\), a inflação do erro Tipo I devido a múltiplos testes. A equação \(\alpha = 1-(1-\alpha_L)^{1/C}\) pode então ser usada para definir \(\alpha_L\) a fim de controlar o valor geral do erro Tipo I.

Como ilustração, suponha que 6 grupos com 100 observações por grupo sejam criados com dados amostrados aleatoriamente de uma população normal. Por construção, a hipóteses nula \(H_0\) é verdadeira, ou seja, todas as médias populacionais são iguais. Agora, construa 5 contrastes independentes desses 6 grupos. Para cada contraste, calcule um teste \(F\).

Se a probabilidade associada ao índice estatístico for menor que \(\alpha = 0.05\), diz-se que o contraste atinge significância, ou seja, \(\alpha\) é usado. Em seguida, faça com que um computador refaça o experimento 10.000 vezes. Em suma, existem 10.000 experimentos, 10.000 famílias de contrastes e \(5\times 10.000 = 50.000\) contrastes.

Os resultados desta simulação são apresentados na tabela a seguir: \[ \begin{array}{ccc} \mbox{Número de famílias com} & X: \mbox{Número de erros Tipo} & \mbox{Número de erros do Tipo I} \\ X \mbox{ erros Tipo I} & \mbox{ 1 por família} & \\\hline 7.868 & 0 & 0 \\ 1.907 & 1 & 1.907 \\ 192 & 2 & 384 \\ 20 & 3 & 60 \\ 13 & 4 & 52 \\ 0 & 5 & 0 \\\hline 10.000 & & 2.403 \\\hline \end{array} \] Tabela 2.3: Resultados de uma simulação de Monte-Carlo. Números de erros do Tipo I ao realizar \(C=5\) contrastes para 10.000 análises de variância realizadas em um projeto de 6 grupos quando \(H_0\) é verdadeiro. Como ler a tabela? Por exemplo, 192 famílias nas 10.000 têm 2 erros do Tipo I, o que dá \(2\times 192 = 384\) erros do Tipo I.

A Tabela 2.3 mostra que o \(H_0\) é rejeitado para 2.403 contrastes nos 50.000 contrastes realmente realizados. A partir desses dados, uma estimativa de \(\alpha\) é calculada como: \[ \alpha = \dfrac{\mbox{Número de contrastes que atingiram significância}}{\mbox{Número total de contrastes}} = \dfrac{2.403}{50.000} = 0.0479\cdot \]

Este valor fica próximo do valor teórico de \(\alpha = 0.05\). Pode-se observar também que em 7.868 experimentos nenhum contraste alcançou significância. Equivalentemente, para 2.132 experimentos (10.000-7.868), pelo menos um erro Tipo I foi cometido.

A partir desses dados, \(\alpha_L\) pode ser estimado como: \[ \alpha_L = \dfrac{\mbox{Número de famílias com pelo menos 1 erro Tipo I}}{\mbox{Número total de famílias}}= = \dfrac{2.132}{10.000} = 0.2132\cdot \]

Este valor se aproxima do valor teórico \[ \alpha_L = 1-(1-\alpha)^C = 1-(1-0.05)^5 = 0.226\cdot \]

Comentamos que contrastes devem ser linearmente independentes, o quê é isso? Dois contrastes são ortogonais ou independentes se seus coeficientes de contraste não estiverem correlacionados. Lembre-se de que os coeficientes de contraste têm soma zero e, portanto, média zero. Portanto, dois contrastes \(L_i\) e \(L_j\), cujos \(t\) coeficientes de contraste, denotados por \(L_{i}(k)\) e \(L_{j}(k)\), \(k=1,\cdots,t\) serão ortogonais se e somente se: \[ \sum_{k=1}^t L_{i}(k) L_{j}(k)=0, \quad i,j =1,\cdots, C, \quad i\neq j \cdot \] Os contrastes ortogonais planejados são equivalentes a perguntas independentes feitas aos dados. Devido a essa independência, o procedimento atual é agir como se cada contraste fosse o único contraste testado. Isso equivale a não usar uma correção para vários testes.

Por exemplo, vamos verificar se os contrastes na equação (2.8) são ortogonais. Primeiro observamos que \(C=2\), \(t=4\), \[ L_{1}(1) =0, \; L_{1}(2)=1, \; L_{1}(3)=-1, \; L_{1}(4)=0 \] e \[ L_{2}(1) =0, \; L_{2}(2)=1, \; L_{2}(3)=0, \; L_{2}(4)=-1\cdot \] Então \[ \begin{array}{l} \sum_{k=1}^4 L_{1}(k) L_{2}(k) & = & 0\times 0 + 1\times 1 + -1\times 0 + 0\times -1 \; = \; 1,\\ \end{array} \] logo, não são ortogonais.

Os contrastes ortogonais são relativamente simples porque cada contraste pode ser avaliado por si só. Os contrastes não ortogonais, entretanto, são mais complexos. O principal problema é avaliar a importância de um determinado contraste em conjunto com os outros contrastes.

Para solucionar este problema existem uas abordagens mencionadas anetriormente. A abordagem clássica corrige múltiplos testes estatísticos, por exemplo, usando uma correção \(\check{\mbox{S}}\)id\(\grave{a}\)k ou Bonferroni, mas essencialmente avalia cada contraste como se viesse de um conjunto de contrastes ortogonais.

Significa que, para a situação na equação (2.8), a implementação R seria

library(gmodels)
fit.contrast( mod0, "time", matrix(c(1,-1,0,1,0,-1),nrow=2,ncol=3,byrow = T) )
##                   Estimate Std. Error   t value   Pr(>|t|)
## time c=( 1 -1 0 )  -2.8125   1.082532 -2.598076 0.02882905
## time c=( 1 0 -1 )  -2.8750   1.082532 -2.655811 0.02622521
## attr(,"class")
## [1] "fit_contrast"

e o o nível de significância corrigido por Bonferroni seria \(2*\alpha \approx \alpha_L\) ou \(\alpha_L/2\approx \alpha\), sendo então que, para concluir acerca da significância do teste em cada contraste utilizamos \(\alpha\approx 0.05/2=0.025\). Assim, nenhum dos contrastes é significativo, isso significa que não pode-se afirmar existam diferenças estatisticamente significativas no crescimento provocadas por deixar mais tempo o pão no forno.


2.4.6 Uma palavra de cautela


Quando se realiza um delineamento completamente aleatório em um fator, o modelo para análise é a Equação (2.1) ou (2.2) e a análise correta é através da análise de variância conforme mostrado simbolicamente na Tabela 2.2. O uso de linguagens de programação como o R facilita a análise de dados e a obtenção de conclusões; entretanto, se o experimento não fosse conduzido adequadamente, mesmo uma análise sofisticada dos dados poderia ser inútil.

O termo de erro \(\epsilon_{ij}\) no modelo (2.1) ou (2.2) e suas somas de quadrados associadas \(ssE\), representam unidades experimentais replicadas. Em muitos casos, os experimentadores não possuem unidades experimentais replicadas em cada nível do fator de tratamento e as substituem por subamostras ou duplicatas na análise. Em outros casos as unidades experimentais não estão devidamente randomizados para níveis de fator de tratamento. Quando esta for a situação, realizar a análise como se o projeto tivesse sido conduzido adequadamente pode ser completamente errado e enganoso.

Podem ser tiradas conclusões erradas que não resistem a um exame posterior e uma má reputação é injustamente atribuída a experiências estatisticamente concebidas e a análises estatísticas de dados.

Por exemplo, considere uma experiência em que um professor gostaria de determinar o efeito dos métodos de ensino nas notas dos testes dos alunos. Se ele usar um método de ensino para as aulas da manhã e outro para as aulas da noite e tratar os resultados dos testes de alunos individuais como réplicas, os resultados de sua análise poderão estar totalmente errados. Esta situação é semelhante aos estudos randomizados por cluster descritos no exercício 5 do Capítulo 1.

A unidade experimental é a turma, uma vez que aplicou o método de ensino a toda uma turma simultaneamente e os alunos individuais são subamostras ou unidades observacionais; uma vez que deve testar os alunos individualmente e não a turma como um todo, O efeito do tratamento deve ser julgado em relação à variabilidade nas unidades ou classes experimentais.

A variabilidade entre os alunos de uma turma pode ser muito diferente da variabilidade entre turmas. Deve-se calcular a média das observações da subamostra antes da análise, conforme explicado na Seção 1.3. Se isso fosse feito, ele teria apenas uma observação por aula e por método de ensino e nenhuma réplica para uso no cálculo de \(ssE\) na Tabela 2.2. Não há denominador para calcular a estatística do teste \(F\) para método de ensino. Se ele usar a variabilidade nos alunos dentro de uma turma para calcular o \(ssE\), ela pode ser muito grande ou muito pequena, fazendo com que ele atinja a conclusão errada sobre o significado do efeito do tratamento.

Além disso, se ele não randomizou qual método de ensino foi usado nas aulas da manhã e da noite e se ele não tiver aulas replicadas que foram ministradas com o mesmo método de ensino, sua análise estará aberta a preconceitos. Os alunos das aulas da manhã podem ser fundamentalmente diferentes dos alunos das aulas da noite e qualquer diferença nas pontuações médias entre os dois métodos de ensino pode ser inteiramente devida a diferenças entre os dois grupos de alunos.

Na verdade, se o professor souber que há diferenças entre os alunos do período matutino e vespertino, ele poderá utilizar propositalmente o método de ensino que deseja promover nas melhores turmas, arruinando assim a objetividade de sua pesquisa.


2.5 Verificando as suposições do modelo linear


Duas suposições necessárias para a validade da análise baseada no modelo linear apresentado na última seção são a constância da variância do erro experimental \(\sigma^2\), em todos os níveis do fator de tratamento e a normalidade dos erros experimentais.

Para verificar essas suposições, gráficos simples podem ser construídos. Um gráfico de dispersão dos resíduos do modelo versus os níveis dos fatores pode mostrar se a variabilidade observada em cada nível do fator é aproximadamente igual. Os resíduos do modelo são as diferenças das respostas \(y_{ij}\) em cada célula ou nível do fator e as médias das células \(\widehat{\mu}_i\). Quando a variância difere entre os níveis dos fatores, muitas vezes é porque a variabilidade na resposta tende a aumentar quando o nível médio da resposta aumenta.

Um gráfico que pode revelar essa tendência é um gráfico dos resíduos do modelo versus as médias das células ou valores previstos. Finalmente, a normalidade dos erros experimentais pode ser verificada fazendo um gráfico de probabilidade normal dos resíduos do modelo. A suposição mais crítica que justifica a análise baseada no modelo linear é a independência dos termos de erro experimental \(\epsilon_{ij}\). Esta suposição é justificada se a randomização adequada das unidades experimentais para os níveis dos fatores de tratamento tiver sido realizada e réplicas verdadeiras forem incluídas.

Um gráfico de dispersão dos resíduos do modelo versus o número da unidade experimental pode revelar inadequações na randomização. Se houver um padrão crescente, decrescente ou cíclico neste gráfico, isso pode indicar que a randomização não equilibrou unidades experimentais heterogêneas entre os níveis do fator. Os quatro gráficos usados para verificar as suposições do modelo linear podem ser facilmente realizados usando R.

O código abaixo produz esses gráficos.

par( mfrow = c(2,2), pch = 19, mar = c(1,4,3,1) )
plot(mod1, which=5); grid()
plot(mod1, which=1); grid()
plot(mod1, which=2); grid()
plot(residuals(mod1) ~ loaf, main="Residuals vs Exp. Unit", font.main=1,data=bread)
abline(h = 0, lty = 2); grid()

Figura 2.4: Gráficos para verificação das premissas do modelo linear

Neste código, o comando R par(mfrow=c(2,2), pch = 19, mar = c(1,4,3,1)) divide a região do gráfico em quatro sub-regiões, pch=19 indica um símbolo ou um único caractere a ser usado como padrão no gráfico de pontos, neste caso um ponto preto cheio e mar = c(1,4,3,1) fornece o número de linhas de margem a serem especificadas nos quatro lados do gráfico, permite que os gráficos ocupem ao máximo cada janela. Os gráficos resultantes estão organizados em linhas na Figura 2.4. O comando plot(mod1, which=5) produz um gráfico dos resíduos padronizados versus os níveis dos fatores no canto superior esquerdo.

O comando plot(mod1, which=1) produz o gráfico dos resíduos versus as médias das células ou valores estimados no canto superior direito. O comando plot(mod1, which=2) produz o gráfico de probabilidade normal ou QQ-plot dos resíduos padronizados no canto inferior esquerdo. A declaração final do gráfico produz o gráfico dos resíduos versus números de unidades experimentais no canto inferior direito. Nesta instrução de gráfico, residuals(mod1) extrai os resíduos do objeto mod1 que foi calculado pela função aov. Uma lista completa das quantidades calculadas pela função aov pode ser obtida digitando names(mod1) no console R.

Se a variabilidade dos resíduos diferir entre os níveis dos fatores no gráfico no canto superior esquerdo da Figura 2.4, isso indicaria que a variância dos \(\epsilon_{ij}\) não é constante. Com apenas quatro réplicas em cada célula, isto é difícil de determinar.

O gráfico dos resíduos versus médias das células, mostrado no canto superior direito da Figura 2.4, pode indicar que a variabilidade nos resíduos aumenta à medida que a média das células aumenta, mas não é claro. Uma maneira melhor de determinar se este é o caso será mostrada na próxima seção. O gráfico de probabilidade normal ou QQ-plot dos resíduos no canto inferior esquerdo justifica a suposição de normalidade relativa aos \(\epsilon_{ij}\)’s se os pontos caírem ao longo de uma linha reta. Quando há mais dados e mais pontos no gráfico, os pontos devem estar mais próximos de uma linha reta para justificar esta suposição. No gráfico normal da Figura 2.4, os pontos afastam-se da linha no canto inferior esquerdo e no canto superior direito, possivelmente indicando caudas curtas na distribuição dos resíduos, mas novamente é difícil determinar com apenas 12 pontos de dados no gráfico.

A suposição de variância igual é mais crítica do que a suposição de normalidade, mas às vezes elas andam de mãos dadas. Quando a suposição de variância igual é violada, a suposição de normalidade também é frequentemente violada e as medidas corretivas usadas para modificar a análise quando há heterogeneidade de variância muitas vezes corrigirão ambos os problemas.


2.6 Estratégias de análise quando as suposições são violadas


Uma causa comum de heterogeneidade de variâncias entre os níveis do fator de tratamento é uma relação não linear entre a resposta e o estímulo ou tratamento. Por exemplo, na metade superior da Figura 2.5, pode-se observar que a resposta aumenta de forma não linear em função dos níveis dos fatores.

As funções de densidade, desenhadas nas laterais em três níveis de tratamento, representam como a não linearidade freqüentemente afeta a distribuição da resposta. À medida que a média ou centro da distribuição aumenta, a variância ou spread na distribuição também aumenta e as distribuições têm caudas longas à direita. Uma forma de corrigir esta situação é transformar os dados de resposta antes da análise.

Figura 2.5: Representação do efeito das não linearidades na distribuição da resposta.

A metade inferior da Figura 2.5 mostra o resultado potencial de uma transformação estabilizadora de variância. Na escala transformada, a variância parece constante em diferentes níveis de fator e as distribuições parecem mais normais.


2.6.1 Transformações de potência Box-Cox


Uma maneira de reconhecer a necessidade de uma transformação estabilizadora de variância é examinar o gráfico de resíduos versus médias celulares descrito na última seção.

Se a dispersão nos resíduos tende a aumentar proporcionalmente em função da média da célula, como possivelmente indicado no canto superior direito da Figura 2.4, geralmente pode ser encontrada uma transformação, \(Y = f(y)\), que resultará em uma análise mais sensível.

Box e Cox (1964) propuseram uma série de transformações de potência \(Y = y^\lambda\) que normalmente funcionam bem. Se a variância tende a aumentar à medida que a média aumenta, escolha um valor de \(\lambda\) menor que um e se a variância tende a diminuir à medida que a média aumenta, escolha um valor de \(\lambda\) maior que um. A Tabela 2.3 resume algumas transformações de potência Box-Cox comuns. Uma situação comum em que \(\sigma \varpropto \mu\) é quando a resposta é na verdade uma medida de variabilidade, como a variância amostral \(s^2\).

\[ \begin{array}{ccc}\hline \mbox{Relação entre } \sigma \mbox{ e } \mu & \lambda & \mbox{Transformação} \\\hline \sigma \varpropto \mu^2 & -1 & \mbox{Recíproca} \\ \sigma \varpropto \mu^{3/2} & -1/2 & \mbox{Raiz quadrada do recíproco} \\ \sigma \varpropto \mu & 0 & \mbox{Logaritmo} \\ \sigma \varpropto \mu^{1/2} & 1/2 & \mbox{Raiz quadrada} \\\hline \end{array} \] Tabela 2.4: Transformações de potência Box-Cox

Em um delineamento completamente aleatório com experimentos replicados em cada nível do fator de tratamento, uma maneira de determinar o valor mais apropriado de \(\lambda\) para usar na transformação Box-Cox é traçar o máximo da função de log verossimilhança, que é proporcional ao recíproco da soma dos quadrados dos erros na ANOVA, versus o valor de \(\lambda\) usado na transformação dos dados.

O valor de \(\lambda\) que maximiza a log-verossimilhança ou que minimiza a soma dos erros dos quadrados, seria o mais apropriado. Este gráfico é chamado de gráfico Box-Cox. A função boxcox no pacote R MASS faz esse gráfico automaticamente.

No exemplo mostrado abaixo, a função boxcox é usada com o objeto lm de nome mod1, que foi ajustado para os dados do experimento de crescimento de pão. O gráfico é mostrado na Figura 2.6 e \(\lambda= -0.5050\) maximiza a log-verossimilhança.

library(MASS)
bc = boxcox(mod1)
lambda = bc$x[which.max(bc$y)]
grid()

Figura 2.6: Gráfico Box-Cox para o experimento de crescimento de pão.


lambda
## [1] -0.5050505


Na Figura 2.6, os valores de \(\lambda\) diretamente abaixo dos pontos onde a linha horizontal pontilhada denominada 95% cruza a curva são limites de confiança de 95% para \(\lambda\).

Neste exemplo, o intervalo de confiança é amplo e inclui \(\lambda= 1\), o que implica sem transformação e \(\lambda = -1\), que implica transformação recíproca. Isto mostra que há uma incerteza considerável sobre a heterogeneidade das variâncias com apenas quatro experiências replicadas em cada nível do factor. Entretanto, para fins ilustrativos, será mostrada a análise utilizando o valor ótimo \(\lambda = −0.5050505\). Segue o código R para adicionar a transformação ao quadro de dados tbread e ajustar o modelo aos dados transformados.

tbread = transform(bread, theight = height^(-.5050505))
mod2 = aov( theight~time, data = tbread )
summary(mod2)
##             Df  Sum Sq  Mean Sq F value Pr(>F)  
## time         2 0.01732 0.008662   6.134 0.0209 *
## Residuals    9 0.01271 0.001412                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


A tabela ANOVA resultante acima mostra que o \(p\)-valor para o fator tempo diminuiu para 0.0209 em relação ao valor 0.042 mostrado na ANOVA anterior dos dados não transformados. Portanto, a transformação tornou a análise um pouco mais sensível.

Para experimentos onde a heterogeneidade da variância é mais pronunciada, a transformação Box-Cox pode aumentar bastante a sensibilidade na detecção dos efeitos do tratamento. Os gráficos para verificar os pressupostos da análise dos dados transformados são mostrados na Figura 2.7.

par( mfrow = c(2,2), pch = 19, mar = c(1,4,3,1) )
plot( mod2, which = 5 ); grid ()
plot( mod2, which = 1 ); grid ()
plot( mod2, which = 2 ); grid ()
plot( residuals(mod2) ~ loaf, main = "Residuals vs Exp. Unit", 
      font.main = 1, data = tbread)
grid ()
abline( h=0, lty = 2 )

Figura 2.7: Gráfico de resíduos versus médias celulares após transformação \(y^\lambda\) para o experimento de crescimento de pão.

Pode-se ver nesta figura que a dispersão ou variabilidade dos resíduos é quase a mesma para cada valor do valor previsto ou média celular das respostas elevadas à potência de -0.505.


2.6.2 Transformações baseadas na distribuição


A suposição de distribuição para o modelo de efeitos para o delineamento completamente aleatório descrito na Seção 2.3 foi \(Y_{ij}\sim N(\mu+\tau_i,\sigma^2)\). Porém, se for sabido que os dados seguem alguma distribuição diferente da distribuição normal, como Binomial, Poisson ou Lognormal, então também se saberia que o desvio padrão não seria constante.

Por exemplo, se a resposta \(Y\), fosse uma contagem binomial do número de sucessos em \(n\) tentativas, então, devido ao Teorema do Limite Central, \(Y\) seria aproximadamente normal, mas \[ \mu_Y = np \qquad \mbox{e} \qquad \sigma_Y = \sqrt{np(1-p)}, \] onde \(p\) é a probabilidade de sucesso. Em situações como esta, onde se sabe que a distribuição da resposta segue alguma forma específica, então uma transformação apropriada pode ser encontrada para estabilizar a variância.

A tabela mostrada a continuação apresenta as transformações recomendadas para situações comuns frequentemente encontradas.

\[ \begin{array}{ccc}\hline \mbox{Distribuuição da resposta} & \mbox{Variuância em termos da média } \mu & \mbox{Transformação } f(y) \\\hline \mbox{Binomial} & \dfrac{\mu(1-\mu)}{n} & \sin^{-1}\Big(\sqrt{y/n}\Big) \; \mbox{(em radianos)} \\ \mbox{Poisson} & \mu & \sqrt{y} \quad \mbox{ou} \quad \sqrt{y+\frac{1}{2}} \\ \mbox{Lognormal} & c \, \mu^2 & \log(y) \\\hline \end{array} \] Tabela 2.5: Transformações baseadas na distribuição da resposta.


2.6.3 Alternativas à análise de mínimos quadrados


Quando a variância do erro experimental não é constante para todos os níveis do fator de tratamento, mas não está relacionada com as médias das células, uma transformação não será uma forma adequada de equalizar ou estabilizar as variâncias.

Uma solução mais geral para o problema é usar mínimos quadrados ponderados. Usando mínimos quadrados ponderados, \(\widehat{\beta}\) é a solução para as equações normais \[ \pmb{X}^\top\pmb{W}\pmb{X}\beta = \pmb{X}^\top\pmb{W}y, \] onde \(\pmb{W}\) é uma matriz diagonal cujos elementos diagonais são os recíprocos do desvio padrão dentro de cada nível de tratamento.

Como ilustração deste método, considere o código R abaixo para analisar os dados do experimento de crescimento de pão.

with(bread, { std <- tapply(height, time, sd)
    weights <- rep( 1/std, each = 4 )
    mod3 <- lm( height ~ time, weights = weights, data = bread )
    anova( mod3 )
 })
## Analysis of Variance Table
## 
## Response: height
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## time       2 18.209  9.1047  6.2263 0.02006 *
## Residuals  9 13.161  1.4623                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### alternate method of analysis using the gls() function in the nlme Package
## An additional example of this method is shown in Section 12.6
library(nlme)
mod_gls = gls(height~time, data=bread, weights=varIdent(form=~1|time))
anova(mod_gls)
## Denom. DF: 9 
##             numDF  F-value p-value
## (Intercept)     1 325.7029  <.0001
## time            2   7.3476  0.0128


Neste exemplo, a função with(bread, f…g) faz com que todas as instruções dentro dos colchetes { } usem as variáveis do quadro de dados bread. A função tapply(height, time, var) é usada para calcular a variância da resposta em cada nível do fator time, isso porque tapply aplica a função var a cada grupo (não vazio) de valores dados height por uma combinação única dos níveis de certos fatores, time nesta situação.

Os pesos são calculados como o inverso dos desvios padrão e a função rep() é usada para expandir o vetor de pesos para o número de linhas no quadro de dados bread. A função lm calcula as estimativas de mínimos quadrados ponderados e a função anova imprime a tabela ANOVA.

Com estes resultados, pode-se observar que o teste \(F\) dos mínimos quadrados ponderados é mais sensível do que os mínimos quadrados não ponderados e o \(p\)-valor é semelhante ao obtido com a transformação Box-Cox mostrada na Seção 2.6.1.

Quando a distribuição do erro não é normal, uma alternativa à análise de uma transformação da resposta é utilizar um modelo linear generalizado (ver McCullagh and Nelder, 1989). Ao ajustar um modelo linear generalizado, o usuário deve especificar a distribuição do erro e uma função de ligação além do modelo. O método de máxima verossimilhança é utilizado para estimar os parâmetros do modelo e os testes da razão de verossimilhança generalizada são utilizados para testar as hipóteses. Quando a função de ligação for a identidade e a distribuição for normal, a análise do modelo linear generalizado resultará no método dos mínimos quadrados e no teste \(F\) ANOVA. Existem várias funções R para ajustar os modelos lineares generalizados e calcular as estatísticas de teste de razão de verossimilhança apropriadas.

Para ilustrar o uso de uma dessas funções para analisar dados experimentais, considere o exemplo a seguir. Um professor queria comparar três métodos de ensino diferentes para determinar como os alunos perceberiam o curso. O fator de tratamento foi o método de ensino, a unidade experimental foi uma turma de alunos e a resposta foi o resumo das avaliações dos alunos do curso. O docente lecionou duas turmas do curso durante três semestres consecutivos resultando num total de seis unidades experimentais ou aulas.

Ele construiu uma lista aleatória para que duas turmas fossem atribuídas a cada método de ensino. Isso reduziria a chance de que outras diferenças nas aulas ou diferenças na execução dos métodos de ensino influenciassem os resultados. Ao final de cada semestre, os alunos foram solicitados a avaliar o curso em uma escala de cinco pontos, sendo 1 o pior e 5 o melhor. Portanto, a resposta de cada classe não foi uma resposta única e normalmente distribuída, \(y\), mas uma resposta vetorial \((y_1,\cdots,y_5)\) que seguiu a distribuição multinomial.

Os dados resumidos do experimento são mostrados na Tabela 2.5.

\[ \begin{array}{ccccccc}\hline \mbox{Classe} & \mbox{Método} & 1 & 2 & 3 & 4 & 5 \\\hline 1 & 1 & 2 & 14 & 12 & 8 & 6 \\ 2 & 3 & 1 & 11 & 15 & 15 & 10 \\ 3 & 2 & 3 & 8 & 18 & 14 & 10 \\ 4 & 3 & 1 & 9 & 17 & 15 & 12 \\ 5 & 2 & 4 & 12 & 19 & 9 & 7 \\ 6 & 1 & 3 & 16 & 13 & 10 & 4 \\\hline \end{array} \] Tabela 2.6: Contagens de pontuações de avaliação dos alunos.


library(daewr);library(MASS)
modf = polr( score ~ method, weight = count, data=teach)
modr = polr( score ~ 1, weight = count, data = teach)
anova(modf,modr)
## Likelihood ratio tests of ordinal regression models
## 
## Response: score
##    Model Resid. df Resid. Dev   Test    Df LR stat.     Pr(Chi)
## 1      1       294   885.9465                                  
## 2 method       292   876.2986 1 vs 2     2  9.64787 0.008035108


Esses dados são armazenados no quadro ou arquivo de dados teach no pacote daewr. O código R acima disponibiliza esses dados e usa a função polr do pacote R MASS (Venables and Ripley, 2002) para ajustar o modelo completo e o modelo reduzido.

A função polr por padrão usa a função de ligação logística e a distribuição multinomial. A resposta score e method é o fator de tratamento no ensino do quadro de dados, enquanto score é uma variável numérica que contém as contagens das várias pontuações de avaliação dos alunos. A fórmula do modelo completo modf, inclui o fator de tratamento, enquanto a fórmula do modelo reduzido strong>modr, inclui apenas o intercepto.

A função anova exibe o teste da razão de verossimilhança da significância do fator de tratamento, conforme mostrado acima.

par( mar = c(1,1,1,1) )
meth = cut(as.numeric(teach$method), c(-Inf, 1, 2, 3), 
           labels = c('Method 1', 'Method 2', 'Method 3'))
class = cbind(teach[ , 3:4], meth)
library(lattice)
barchart(count ~ score | meth, data = class, layout = c(1,3), 
         horizontal = FALSE, xlab = "Score", col = "grey")

Figura 2.8: Porcentagem de notas dos alunos por método de ensino

O \(p\)-valor para a estatística qui-quadrado da razão de verossimilhança é pequeno, indicando que há uma diferença significativa entre os métodos de ensino. O método de ensino 1 teve nota média de 2.98, o método de ensino 2 teve nota média de 3.22 e o método de ensino 3 pareceu ser o melhor com nota média de 3.47.

Isto também pode ser visualizado nos gráficos de barras da Figura 2.8, que mostram que a percentagem de pontuações altas atribuídas aumenta para os métodos de ensino 2 e 3.


2.7 Determinando o número de réplicas


O nível de significância \(\alpha\) do teste \(F\) ANOVA de ausência de efeito de tratamento é a probabilidade de rejeição da hipótese nula \[ H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_t, \] quando for verdade.

O poder do teste é a probabilidade de rejeitar a hipótese nula quando esta é falsa. A estatística de teste \(msT/msE\) segue a distribuição \(F\) quando a hipótese nula é verdadeira, mas quando a hipótese nula é falsa ela segue a distribuição \(F\) não central. A distribuição \(F\) não central tem uma distribuição mais ampla do que a distribuição \(F\) central, como mostra a Figura 2.9.

Figura 2.9: Distribuição \(F\) central e não central.

A propagação na distribuição \(F\) não central e a probabilidade que exceder o limite crítico da distribuição \(F\) central é uma função crescente do parâmetro de não centralidade \(\lambda\). Quando a distribuição é \(F\) não central, a probabilidade de exceder o limite crítico da distribuição \(F\) central é chamada de potência.

O poder é maior que o nível de significância \(\alpha\) quando a hipótese nula é falsa tornando o parâmetro de não centralidade maior que zero. O poder pode ser calculado para qualquer cenário de médias diferentes, se os valores das médias da célula, a variância do erro experimental e o número de repetições por nível de fator forem especificados.

Para uma diferença constante entre as médias das células, representada por \[ \sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet}\big)^2, \] o parâmetro de não centralidade e o poder aumentam à medida que o número de repetições aumenta.

Quando as diferenças entre as médias celulares são grandes o suficiente para ter importância prática, o experimentador gostaria de ter alto poder ou probabilidade de rejeitar a hipótese de ausência de efeitos do tratamento. Quando a diferença entre as médias tem importância prática para o pesquisador, chamamos isso de significado prático. A significância prática nem sempre corresponde à significância estatística determinada pelo teste \(F\) da ANOVA.

Às vezes, o número de repetições no experimento é muito pequeno e a probabilidade ou o poder de detectar uma diferença de significância prática é muito baixo. A significância estatística pode ser feita para coincidir com a significância prática, determinando o número apropriado de réplicas que resultam no poder desejado. Fazer isso é a segunda técnica que se enquadra na categoria de controle de erros discutida no Capítulo 1. A ideia de que aumentar o número de repetições aumenta a sensibilidade do experimento também se deve a R.A. Fisher (1935).

Por exemplo, se houver uma diferença entre as médias das células de modo que a soma dos quadrados corrigida \(css =\sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet}\big)^2\) seja maior que zero, então o poder ou probabilidade de rejeitar \[ H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_t \] é dado por \[\begin{equation*} \tag{2.10} \pi(\lambda)=\int_{}^\infty F(x,t-1,t(r-1),\lambda)\mbox{d}x \end{equation*}\] onde \(F_{t-1,t(r-1),\alpha}\) é o \(\alpha\)-ésimo percentil da distribuição \(F\) central com \(t-1\) e \(t(r − 1)\) graus de liberdade, \(F(x, t-1, t(r-1),\lambda)\) é a distribuição \(F\) não central com parâmetro de não centralidade \[ \lambda = \dfrac{r}{\sigma^2} \sum_{i=1}^t \big( \mu_i-\overline{\mu}_{\bullet}\big)^2\cdot \]

Para um valor fixo de \(\frac{1}{\sigma^2}\sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet}\big)^2\), o parâmetro de não centralidade e o poder aumentam em função do número de repetições \(r\). Esta probabilidade pode ser calculada para vários valores de \(r\) até que um valor com potência adequada seja encontrado. Desta forma, o número apropriado de réplicas pode ser determinado.

A função Fpower no pacote R daewr facilita esses cálculos.

library(daewr)
rmin = 2 # smallest number of replicates
rmax = 6 # largest number of replicates
alpha = rep(0.05, rmax - rmin +1)
sigma = sqrt(2.1)
nlev = 3
nreps = rmin:rmax
Delta = 3
power = Fpower1(alpha, nlev, nreps, Delta, sigma)  
power
##      alpha nlev nreps Delta    sigma     power
## [1,]  0.05    3     2     3 1.449138 0.1947995
## [2,]  0.05    3     3     3 1.449138 0.4041857
## [3,]  0.05    3     4     3 1.449138 0.5903406
## [4,]  0.05    3     5     3 1.449138 0.7328895
## [5,]  0.05    3     6     3 1.449138 0.8329923


No experimento de crescimento de pão, suponha que uma diferença de menos de 3 polegadas (7 cm) na altura da massa levedada não tenha consequências. No entanto, se alterar o tempo de crescimento de 35 minutos para 45 minutos causar uma diferença de mais de 3 polegadas (7 cm) na altura da massa levedada, o experimentador gostaria de saber sobre isso, porque ele precisará monitorar de perto o tempo de crescimento no futuro para produzir pães de altura consistente.

Neste caso, podemos considerar \(\Delta= 3.0\) como uma diferença prática nas médias das células. O menor \(css = \sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet}\big)^2\) poderia ser, com pelo menos duas médias de células diferindo em \(\Delta\), seria o caso quando a média de uma célula fosse \(\Delta/2\) maior que a média geral, uma segunda fosse \(\Delta/2\) a menos que a média geral e um terço era igual à média geral.

Isto resultaria em \[ css = \sum_{i=1}^t \big(\mu_i-\overline{\mu}_{\bullet}\big)^2 = \left(\dfrac{\Delta}{2}\right)^2+0^2+\left(-\dfrac{\Delta}{2}\right)^2 = \dfrac{\Delta^2}{2}+\dfrac{3^2}{2}=4.5\cdot \]

Supondo que a variância do erro experimental \(\widehat{\sigma}^2 = 2.1\) foi estimada a partir da variância amostral nas alturas da massa crescida em um experimento piloto onde vários pães foram autorizados a crescer pelo mesmo período de tempo, então o fator de não centralidade pode ser calculado como \(\lambda=\frac{r}{2.1}\times (4.5)\).

A potência é calculada para \(r = 2,\cdots,6\) usando o código R mostrado acima. Este código ilustra o uso da função Fpower1 que toma como argumentos alfa=\(\alpha\), nlev=\(t\), o número de níveis do fator, nreps=\(r\), Delta=\(\Delta\) e sigma=\(\sigma\).

Ao usar um argumento vetorial para nreps, a função produz um vetor correspondente de valores de power. A partir disso podemos ver que com \(r = 5\) réplicas haveria 73% de chance de detectar uma diferença nas médias celulares tão grande quanto 3.0 e com \(r = 6\) há 83% de chance. Com menos de cinco réplicas, há pelo menos 40% de probabilidade de esta diferença ser perdida. Como regra geral, o número de réplicas que resulta em poder entre 0.80 e 0.90 é geralmente suficiente para a maioria dos estudos experimentais.


2.8 Comparação de tratamentos após o teste \(F\)


Quando o teste \(F\) para a hipótese nula \[ H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_t \] é rejeitado, isso nos diz que existem diferenças significativas entre pelo menos duas das médias das células, mas se houver vários níveis de Se considerarmos o fator tratamento, isso não significa necessariamente que todas as médias celulares sejam significativamente diferentes umas das outras.

Quando a hipótese nula é rejeitada, investigações adicionais devem ser realizadas para descobrir exatamente quais médias de células diferem. Em alguns casos, o investigador terá comparações pré-planejadas que gostaria de fazer; em outras situações, ele pode não ter ideia de quais diferenças procurar.

Contrastes planejados ou a priori são selecionados antes da execução do experimento. Em geral, eles refletem as hipóteses que o experimentador queria testar e geralmente são poucos. Os contrastes não planejados ou a posteriori, após o fato, são decididos após a execução do experimento. O objetivo dos contrastes a posteriori é garantir que resultados inesperados sejam confiáveis.


2.8.1 Comparações pré-planejadas


Considerando os níveis dos factores de tratamento na experiência de produção de beterraba sacarina realizada em Rothamstead em 1937 e descrita na Secção 2.3, algumas comparações pré-planeadas que podem ter sido de interesse são:

  1. \[ H_0 \, : \, \mu_1=\dfrac{1}{3}\big(\mu_2+\mu_3+\mu_4\big) \]

  2. \[ H_0 \, : \, \mu_2=\mu_3 \]

  3. \[ H_0 \, : \, \mu_3=\mu_4 \]

A primeira comparação levanta a questão: uma mistura de fertilizantes artificiais altera o rendimento? A segunda comparação levanta a questão: existe uma diferença nos rendimentos entre a aplicação de fertilizante artificial arado e a lanço? A terceira comparação levanta a questão: o momento da aplicação altera o rendimento?

Todas essas hipóteses podem ser expressas na forma geral \[ H_0 \, : \, \sum_{i=}^t c_i \mu_i = 0, \] onde \(\sum_{i=1}^t c_i = 0\)

Dado que \(\sum_{i=1}^t c_i \mu_i=0\) são funções estimáveis, cada uma das hipóteses anteriores pode ser testada calculando a única função estimável \(\pmb{L}\widehat{\beta}\) e seu erro padrão é \[ s_{\pmb{L}\widehat{\beta}} = \sqrt{\widehat{\sigma}^2 \pmb{L}^\top \big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{L}}\cdot \] A razão entre a função estimável e seu erro padrão segue a distribuição \(t\)-Student.

A função fit.contrast no pacote R gmodels realiza esse teste. Para o experimento da beterraba sacarina, o código abaixo carrega os dados, configura e imprime a matriz de contraste \(\pmb{L}\).

library(daewr)
mod4 = aov( yield ~ treat, data = sugarbeet )
con = matrix(c(1, -1/3, -1/3, -1/3, 0, 1, -1, 0, 0, 0, 1, -1 ), 4, 3 )
L = t(con)
rownames(L) = c("-fertilizer effect", "-plowed vs. broadcast", "-January vs. April")
L
##                       [,1]       [,2]       [,3]       [,4]
## -fertilizer effect       1 -0.3333333 -0.3333333 -0.3333333
## -plowed vs. broadcast    0  1.0000000 -1.0000000  0.0000000
## -January vs. April       0  0.0000000  1.0000000 -1.0000000


A chamada de função abaixo imprime os resultados. A instrução options controla o número de dígitos após o decimal para impressão.

options(digits = 3)
library(gmodels)
fit.contrast( mod4, "treat", L)
##                            Estimate Std. Error t value Pr(>|t|)
## treat-fertilizer effect        -8.8      0.825 -10.664 4.19e-08
## treat-plowed vs. broadcast     -3.8      0.975  -3.897 1.61e-03
## treat-January vs. April         0.1      0.919   0.109 9.15e-01
## attr(,"class")
## [1] "fit_contrast"


Os \(p\)-valores na coluna rotulada Pr(>|t|), na saída acima, podem ser interpretados da mesma forma que os \(p\)-valores para a estatística \(F\) foram interpretados e podemos ver que:

  1. os fertilizantes artificiais aumentam o rendimento,

  2. a aplicação a lanço resulta em rendimentos mais elevados do que a aplicação arada, e

  3. há nenhuma diferença significativa no rendimento entre o período de aplicação de abril e janeiro.

Quando os níveis dos fatores são quantitativos, como o tempo de crescimento no experimento de crescimento da massa de pão, as comparações pré-planejadas geralmente envolvem a busca pela significância de tendências polinomiais lineares ou de ordem superior na resposta.

Em alguns experimentos, os tratamentos estão associados a níveis numéricos como dose do medicamento, tempo de cozimento ou temperatura de reação. Iremos nos referir a tais níveis como doses, não importa quais sejam realmente, e o valor numérico da dose para tratamento \(i\) será denotado por \(z_i\).

Quando temos doses numéricas, podemos reexpressar a média do tratamento em função da dose \(z_i\) como \[ \mu+\tau_i = f(z_i;\theta), \] onde \(\theta\) é algum parâmetro desconhecido da função. Por exemplo, poderíamos expressar o peso médio de mudas em função do pH da chuva ácida.

As funções \(f\) mais comumente usadas são polinômios na dose \(z_i\) \[ \mu+\tau_i = \theta_0+\theta_1 z_i +\theta_2 z_i^2+\cdots+\theta_{g-1} z_i^{g-1}\cdot \]

Usamos a potência \(g-1\) porque as médias em \(g\) doses diferentes determinam um polinômio de ordem \(g-1\). Polinômios são usados com frequência porque são simples e fáceis de entender; nem sempre são a escolha mais adequada.

Se conhecermos os coeficientes polinomiais \(\theta_0,\theta_1,\cdots, \theta_{g-1}\), então podemos determinar as médias de tratamento \(\mu+\tau_i\) e vice-versa. Se conhecermos os coeficientes polinomiais exceto a constante \(\theta_0\), então podemos determinar os efeitos do tratamento \(\tau_i\), e vice-versa.

Os \(g-1\) parâmetros \(\theta_1\) a \(\theta_{g-1}\) neste modelo polinomial completo correspondem aos \(g-1\) graus de liberdade entre os grupos de tratamento. Assim, os polinômios na dose não são inerentemente melhores ou piores que o modelo de efeitos do tratamento, apenas outra forma de descrever as diferenças entre as médias.

Por experiência, consideramos a modelagem polinomial é útil em dois contextos. Primeiro, se apenas alguns dos coeficientes polinomiais forem necessários, ou seja, os outros podem ser definidos como zero sem diminuir significativamente a qualidade do ajuste, então este modelo polinomial reduzido representa uma redução na complexidade do nosso modelo. Por exemplo, aprender que a resposta é linear ou quadrática em dose é útil, enquanto um polinômio de grau seis ou sete será difícil de compreender ou vender para qualquer outra pessoa.

Em segundo lugar, se quisermos estimar a resposta numa dose diferente da utilizada na experiência, o modelo polinomial fornece um mecanismo para gerar as estimativas. Observe que essas estimativas podem ser ruins se extrapolarmos além da faixa de doses em nosso experimento ou se o grau do polinômio for alto. Polinômios de ordem superior se ajustarão exatamente às nossas médias de tratamento observadas, mas esses polinômios de ordem superior podem ter um comportamento bizarro longe de nossos pontos de dados.

Considere uma sequência de modelos de regressão para nossos dados, regredindo as respostas em dose, dose ao quadrado e assim por diante. O primeiro modelo inclui apenas a constante \(\theta_0\); isto é, cabe um único valor para todas as respostas. O segundo modelo inclui a constante \(\theta_0\) e um termo linear \(\theta_1 z_i\); este modelo ajusta as respostas como uma regressão linear simples em dose. O terceiro modelo inclui a constante \(\theta_0\), um termo linear \(\theta_1 z_i\) e o termo quadrático \(\theta_2 z_i^2\); este modelo ajusta as respostas como uma função quadrática (parábola) da dose. Modelos adicionais incluem potências adicionais de dose até \(g-1\).

Um método para escolher um modelo polinomial é escolher a menor ordem, de modo que nenhum termo significativo seja excluído. Existem métodos de seleção de modelos mais sofisticados. É importante saber que os coeficientes estimados θbi dependem de quais termos estão no modelo quando o modelo é estimado.

Assim, se decidirmos que só precisamos de \(\theta_0\), \(\theta_1\) e \(\theta_2\) quando \(g\) for 4 ou mais, devemos reajustar usando apenas esses termos para obter estimativas de parâmetros apropriadas. Um truque adicional a lembrar ao construir um modelo dose-resposta é que podemos transformar ou reexpressar a dose \(z_i\). Ou seja, podemos construir modelos utilizando o logaritmo da dose ou a raiz quadrada da dose tão simplesmente como podemos utilizar a dose. Para alguns dados é muito mais simples construir um modelo em função de uma transformação da dose.

Coeficientes de contraste para testar tendências polinomiais ortogonais, podem ser obtidos a partir da função R contr.poly. A entrada necessária para esta função é o número de níveis do fator. O resultado é uma matriz ortogonal com os coeficientes de contraste desejados. Por exemplo, para o experimento de crescimento de massa de pão, os comandos mostrados a seguir constroem e imprimem a matriz de contraste.

A matriz de contraste resultante abaixo possui coeficientes para os contrastes lineares e quadráticos.

contrasts(bread$time) = contr.poly(3)
contrasts(bread$time)
##           .L     .Q
## 35 -7.07e-01  0.408
## 40 -7.85e-17 -0.816
## 45  7.07e-01  0.408


O código que usa as funções R aov e summary.lm mostradas abaixo calcula os contrastes e exibe os resultados. Nos resultados a seguir, podemos ver que há uma tendência linear significativa, no nível \(\alpha= 0.05\), mas nenhuma tendência quadrática significativa.

mod3 = aov( height ~ time, bread )
summary.lm(mod3)
## 
## Call:
## aov(formula = height ~ time, data = bread)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -1.81  -1.14   0.00   1.27   2.25 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.333      0.442   16.59  4.7e-08 ***
## time.L         2.033      0.765    2.66    0.026 *  
## time.Q        -1.123      0.765   -1.47    0.177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.53 on 9 degrees of freedom
## Multiple R-squared:  0.506,  Adjusted R-squared:  0.396 
## F-statistic:  4.6 on 2 and 9 DF,  p-value: 0.042


Em termos na ANOVA, mostramos a continuação como apresentar os resultados dos termos linear e quadrático:

summary.aov(mod3, split = list(time = list("Linear" = 1, "Quadrático" = 2)))
##                    Df Sum Sq Mean Sq F value Pr(>F)  
## time                2  21.57   10.79    4.60  0.042 *
##   time: Linear      1  16.53   16.53    7.05  0.026 *
##   time: Quadrático  1   5.04    5.04    2.15  0.177  
## Residuals           9  21.09    2.34                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Os resultados finais mostram a falta de comportamento quadrático e, portanto, o aumento na altura do crescimento do pão não é uma função linear do aumento do tempo de cocimento. Entre 35min e 40min de forno há um aumento no crescimento do pãp, depois de estabiliza o qual significa que não adianta deixar o pão mais do que 40min no forno, a massa do pão não cresce mais.

Os valores preditos podem ser obtidos utilizando o modelo estimado, como mostrado a seguir:

modelo = coef(mod3)[1]+coef(mod3)[2]*contrasts(bread$time)[,1]+coef(mod3)[3]*contrasts(bread$time)[,2]
modelo
##   35   40   45 
## 5.44 8.25 8.31


Uma outra maneira de enocntar os valores preditos é utilizando a função predict, como mostrado a seguir:

par(mar=c(4,4,1,1))
plot(height ~ time, bread, xlab = "Tempo", ylab = "Altura")
points(height ~ time, bread, pch =19 )
grid()
lines(predict(mod3) ~time, bread)


Se os níveis para o fator time foram criados com o comando ordered em vez do comando factor, acontece que o R cria automaticamente a matriz \(\pmb{X}\) usando os contrastes polinomiais ortogonais e a tabela de resumo pode ser obtida sem criar contrastes adicionais para time.


2.8.2 Comparações não planejadas


Quando um conjunto de comparações pré-planejadas pode ser expresso como um conjunto saturado de contrastes ortogonais, como os exemplos mostrados na última seção, essas comparações são independentes e equivalentes a particionar o teste \(F\) geral de \[ H_0 \, : \, \mu_1 = \cdots = \mu_t\cdot \] Contudo, se as comparações não forem planeadas antes da execução da experiência, o analista pode ficar tentado a escolher as comparações que gostaria de fazer com base nas médias dos dados.

Isto significa implicitamente que todas as comparações possíveis foram feitas. Ao testar todas as comparações possíveis, cada uma no nível de significância \(\alpha=0.05\), o nível de significância geral pode ser muito superior a 0.05, superior a 50% em alguns casos.

Isso significa que mesmo quando não há diferença nas médias \(\mu_1,\cdots,\mu_t\) pode haver uma alta probabilidade de encontrar uma ou mais comparações significativas quando cada uma é testada individualmente. Para reduzir a chance geral ou experimental de um erro tipo I, um ajuste deve ser feito.

Para comparações pareadas da forma \[ H_0 \, : \, \mu_i = \mu_j \qquad \mbox{para} i \neq j, \] teste HSD de Tukey ou diferença honestamente significativa (Honestly Significant Difference) ajusta a região crítica usando a estatística de rango padronizada em vez da distribuição \(t\)-Student. Usando o HSD rejeitamos \(H_0 \, : \, \mu_i = \mu_j\) em favor da alternativa \(H_1 \, : \, \mu_i\neq \mu_j\) se \[ | \, \widehat{\mu}_i-\widehat{\mu}_j \, | > q_{I,n-t,\alpha} s_{\widehat{\mu}_i-\widehat{\mu}_j}\sqrt{2} \] onde \(q_{I,n-t,\alpha}\) é o \(\alpha\) percentil superior da rango padronizado.

Isso só é aproximado quando os tamanhos das amostras são desiguais. Se \(X_1,\cdots,X_I\) são variáveis aleatórias independentes com distirbuição \(N(\mu,\sigma^2)\) e \[ R = \max_i\{X_i\} − \min_i\{X_i\} \] então \(R/\widehat{\sigma}\) segue a distribuição de rango padronizada (ver Tukey, 1949a).

A função R TukeyHSD fará comparações aos pares usando o método HSD de Tukey. O código abaixo ilustra como esta função é chamada para fazer comparações nos dados do experimento da beterraba sacarina.

mod4 = aov( yield ~ treat, data = sugarbeet )
mod4.tukey = TukeyHSD( mod4, ordered = T )
mod4.tukey
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = yield ~ treat, data = sugarbeet)
## 
## $treat
##     diff    lwr   upr p adj
## B-A  6.3  3.312  9.29 0.000
## D-A 10.0  7.166 12.83 0.000
## C-A 10.1  7.266 12.93 0.000
## D-B  3.7  0.866  6.53 0.009
## C-B  3.8  0.966  6.63 0.008
## C-D  0.1 -2.572  2.77 1.000


A primeira coluna da saída lista a comparação feita, a próxima coluna lista a diferença nas médias das células e as próximas duas colunas são limites para um intervalo de confiança de 95% na diferença das médias da forma \[ | \, {\mu}_i-{\mu}_j \, | \pm q_{I,n-t,0.05} s_{\widehat{\mu}_i-\widehat{\mu}_j}\sqrt{2}\cdot \]

A coluna final é um \(p\)-valor para o teste da hipótese nula de que as duas médias são iguais. Por exemplo, o intervalo de confiança para a última comparação, \(\mu_C - \mu_D\), inclui zero e o \(p\)-valor é > 0.05, indicando que o rendimento da beterraba sacarina para tratamento C - transmissão artificial aplicada em janeiro, não é significativo. ligeiramente diferente do rendimento para tratamento D - transmissão artificial aplicada em abril. Todas as outras comparações aos pares mostram uma diferença altamente significativa.

Um método menos conservador de comparar todas as médias celulares possíveis foi desenvolvido independentemente por Newman (1939) e Keuls (1952). Este método também se baseia na estatística de intervalo padronizada, mas é baseado no intervalo do par específico de médias que está sendo comparado, dentro de todo o conjunto de médias ordenadas, em vez do intervalo do maior para o menor como o método ou teste HSD de Tukey.

A comparação de médias utilizando o método de Newman-Keuls pode ser feita utilizando a função Snk.test do pacote R agricolae (de Mendiburu, 2012a). Os argumentos para a função Snk.test são semelhantes aos da função TukeyHSD e são ilustrados abaixo usando os dados do experimento com beterraba sacarina.

library(agricolae)
compare = SNK.test( mod4, "treat", alpha = 0.05 )
print(compare)
## $statistics
##   MSerror Df Mean   CV
##      2.11 14 45.7 3.18
## 
## $parameters
##   test name.t ntr alpha
##    SNK  treat   4  0.05
## 
## $snk
## NULL
## 
## $means
##   yield  std r    se  Min  Max  Q25  Q50  Q75
## A  38.7 1.65 4 0.727 36.9 40.1 37.5 38.9 40.1
## B  45.0 1.33 4 0.727 43.7 46.2 43.9 45.0 46.1
## C  48.8 1.10 5 0.650 47.7 50.6 48.2 48.6 48.9
## D  48.7 1.68 5 0.650 47.2 51.3 47.3 48.5 49.2
## 
## $comparison
## NULL
## 
## $groups
##   yield groups
## C  48.8      a
## D  48.7      a
## B  45.0      b
## A  38.7      c
## 
## attr(,"class")
## [1] "group"


A seção de faixa crítica da saída lista os valores críticos para diferenças de médias que variam em 2, 3 ou 4. Na última secção do resultado, as médias com o mesmo indicador de Grupo à esquerda não são significativamente diferentes. Isto mostra que o rendimento da beterraba sacarina para tratamento (difusão aplicada artificialmente em janeiro) não é significativamente diferente do rendimento para tratamento (difusão aplicada artificialmente em abril). Todas as outras comparações aos pares mostram uma diferença significativa, neste caso, os mesmos resultados do método HSD de Tukey.

A última seção da saída da função Snk.test ilustra uma maneira compacta de apresentar as diferenças significativas nas médias dos tratamentos que são encontradas por técnicas de comparação múltipla, como o método HSD de Tukey ou o método de Newman-Keuls.

Ao relatar resultados em texto escrito, este método de apresentação pode ser modificado listando as médias horizontalmente no texto, da menor para a maior, e usando um sublinhado no lugar do indicador do Grupo para mostrar quais médias não são significativamente diferentes. Médias que não são significativamente diferentes são sublinhadas com a mesma linha. O exemplo abaixo mostra as médias de um experimento para determinar o efeito do site de download no momento de baixar um arquivo.

Os resultados mostram que o tempo de download não é significativamente diferente entre os sites B e D, e não significativamente diferente entre os sites D e A, mas há uma diferença significativa no tempo de download entre os sites B e A.

Da mesma forma, não há diferença significativa nos tempos de download para os sites A e C, mas o tempo de download para o site C é significativamente maior do que para o site B ou D. Finalmente, o site E tem um tempo de download significativamente mais longo do que qualquer um dos outros sites.


2.8.3 Comparação de todas as médias com um controle ou o melhor


Em alguns experimentos, um dos níveis de tratamento é o nível atual ou padrão e os outros são níveis novos ou experimentais. Um dos principais objetivos neste tipo de experimento pode ser comparar a média de cada nível experimental com o padrão ou às vezes chamado de nível de controle. Dunnett (1955) desenvolveu um método para fazer isso e controlar a taxa de erro tipo I em termos de experimento.

Na experiência de produção de beterraba sacarina, o nível de tratamento (sem fertilizante artificial) pode ser considerado como o controle. Todos os outros níveis de tratamento podem ser comparados a este usando a função glht no pacote R multcomp (Hothorn et al., 2008).

Para carregar o pacote multcomp você também deve ter o pacote mvtnorm (Genz et al., 2012) e o pacote survival (Therneau, 2012) instalados. Ao usar o método Dunnett, a função glht, por padrão, usa o primeiro nível do fator de tratamento como controle.

summary(sugarbeet)
##  treat     yield     
##  A:4   Min.   :36.9  
##  B:4   1st Qu.:43.8  
##  C:5   Median :47.2  
##  D:5   Mean   :45.7  
##        3rd Qu.:48.6  
##        Max.   :51.3


Como pode ser visto acima, o nível de tratamento A - sem fertilizante artificial é o primeiro nível do fator de tratamento no quadro de dados para o experimento de produção de beterraba sacarina. O código para usar o método de Dunnett para comparar a média em cada nível do fator de tratamento com o controle (A) chamando a função glht é mostrado a seguir.

library(multcomp)
sugar.dun = glht(mod4,linfct = mcp(treat = "Dunnett"), alternative = "greater")
summary(sugar.dun)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = yield ~ treat, data = sugarbeet)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value  Pr(>t)    
## B - A <= 0    6.300      1.028    6.13 1.3e-05 ***
## C - A <= 0   10.100      0.975   10.36 < 1e-05 ***
## D - A <= 0   10.000      0.975   10.25 < 1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)


A saída acima é o resultado de testes unicaudais. Ao comparar todos níveis de tratamento para um controle, a direção desejada da diferença é frequentemente conhecida. Portanto, um teste unilateral, em vez de um teste bicaudal, pode ser necessário. Outras opções para alternative = no código acima são “less” (menos) ou “two.sided” (dois lados).

Quando não existe um nível de controle do fator de tratamento, ainda pode haver interesse em comparar todos os níveis de tratamento com o melhor nível. Por exemplo, no experimento para ver o efeito do site de download no tempo de download de um arquivo, descrito no final da Seção 2.8.2, pode ser interessante encontrar todos os sites cujos tempos médios de download não são significativamente significativos. mais longo que o site com o tempo médio mínimo de download observado. Nas médias mostradas na Seção 2.8.2, o site B teve o menor tempo de download observado.

Para comparar todas as médias de tratamento ao melhor nível e controlar a taxa de erro experimental, o procedimento MCB de Hsu (1984) pode ser usado. Este procedimento é equivalente ao método de Dunnett. Para utilizar este método, primeiro observe as médias observadas e decida qual é a melhor. A seguir, configure contrastes comparando cada nível com o melhor. Finalmente, chame a função glht para realizar o teste de Dunnett. Por exemplo, se os dados para o experimento de download de arquivos estivessem contidos em um quadro de dados chamado download, e o segundo nível (ou “B”) do fator site tivesse o tempo médio mínimo de download, o código abaixo configura os contrastes e chama a função glht para comparar as médias de tratamento dos locais “A”, “C”, “D” e “E” com a média do local “B” usando o método de Dunnett.

#aov.ex = aov(time ~ site, data = download)
K = rbind( c( 1, -1, 0, 0, 0), 
           c(0, -1, 1, 0, 0), 
           c(0, -1, 0, 1, 0), 
           c(0, -1, 0, 0, 1) )
rownames(K) = c( "A-B", "C-B", "D-B", "E-B" )
#colnames(K) = names(coef (aov.ex))
#dht = glht( aov.ex, linfct = mcp( site = "Dunnett" ), alternative = "two.sided")
#summary(dht)

2.9 Exercícios