3.1 Introdução


No Capítulo 2 examinamos projetos unifatoriais. Eles são úteis apenas quando um fator está em estudo. Quando vários fatores estão em estudo, uma abordagem clássica é estudar cada um separadamente, mantendo todos os outros constantes. R.A. Fisher (1935) destacou que esta abordagem é útil para demonstrar relações conhecidas aos alunos em cursos de laboratório quando a influência de outros fatores é conhecida. No entanto, esta abordagem é ineficiente e potencialmente enganosa quando se trata de descobrir novos conhecimentos através da experimentação.

Uma estratégia muito melhor para experimentar múltiplos fatores é usar um planejamento fatorial. Num planejamento fatorial as células consistem em todas as combinações possíveis dos níveis dos fatores em estudo. Os planejamentos fatoriais acentuam os efeitos dos fatores, permitem estimar a interdependência dos efeitos ou interações e são a primeira técnica na categoria do que é chamado de planejamento de tratamento.

Ao examinar todas as combinações possíveis de níveis de fator, o número de réplicas de um nível específico de um fator é aumentado pelo produto do número de níveis de todos os outros fatores no projeto e, assim, o mesmo poder ou precisão pode ser obtido com menos replicas. Além disso, se o efeito de um fator mudar dependendo do nível de outro fator, isso será visto em um plano fatorial. Este fenômeno passará despercebido na abordagem clássica, onde cada fator só varia em níveis constantes dos outros fatores. O exemplo da próxima seção ilustrará esses fatos.


3.2 Plano fatorial clássico, um de cada vez versus planos fatoriais


Nos exercícios do Capítulo 2, foi descrito um conjunto de experimentos com helicópteros de papel. Nessas experiências apenas um factor, o comprimento da asa, estava em estudo. Contudo, para maximizar o tempo de voo dos helicópteros de papel, seria aconselhável considerar mais de um fator. Por exemplo, considere variar o comprimento da asa em 4 níveis como antes e a largura do corpo em quatro níveis, como 4.25”, 4.0”, 3.75” e 3.5”.

O lado esquerdo da Figura 3.1 representa o plano clássico no qual um fator varia de cada vez. Os círculos no diagrama representam experimentos ou execuções. Usando esta abordagem, os experimentos na parte inferior da figura seriam concluídos variando o comprimento da asa (wing length), mantendo a largura do corpo (body width) constante em 3.5”. Em seguida, os três experimentos adicionais no lado esquerdo da figura seriam completados variando a largura do corpo, mantendo o comprimento da asa constante em seu nível baixo de 4.0”.

Se fossem feitas oito repetições para cada um desses experimentos, como sugerido no exercício 2 do Capítulo 2, seriam necessários um total de 56 experimentos.

Se o objetivo fosse encontrar a combinação com o maior tempo de voo, a abordagem clássica seria completar os experimentos primeiro com um fator. O próximo seria calcular as médias das células e então selecionar o nível com a média mais alta. Finalmente, o segundo factor seria variado, mantendo a primeira constante no seu nível óptimo e não no nível mais baixo, como mostra a Figura 3.1.

No entanto, a cautela de R.A. Fisher em randomizar diria que esta é uma má estratégia. Se alguma força desconhecida mudasse após o primeiro conjunto de experimentos, os resultados poderiam ser tendenciosos. Além disso, o nível ideal de um fator pode depender do nível do outro fator. Portanto, ao variar um fator de cada vez, o ótimo geral pode ser perdido.

Figura 3.1: Comparação entre experimentos individuais e fatoriais.

O diagrama do lado direito da Figura 3.1 representa um plano fatorial para os experimentos com helicópteros. Aqui pode ser visto que os experimentos são realizados em todas as combinações dos níveis dos dois fatores. Neste plano, se duas réplicas de cada célula forem concluídas, haverá oito réplicas de cada nível de comprimento de asa e oito réplicas de cada nível de largura do corpo, o que equivale ao plano um de cada vez. Portanto, o plano fatorial teria a mesma precisão ou poder para detectar efeitos de fator que o plano um de cada vez, mas é mais eficiente, pois requer apenas \(2\times 16 = 32\) execuções totais, em oposição às 56 exigidas pelo plano um de cada vez. O número de réplicas de cada nível de fator no experimento fatorial é igual ao número de réplicas por célula vezes o produto dos níveis de todos os outros fatores no experimento. Essa multiplicação é chamada de replicação oculta.

No caso mostrado na Figura 3.1, existem apenas dois fatores, cada um em quatro níveis; portanto, o número de repetições de cada nível de fator é \(2\times 4 = 8\). No plano fatorial, as 32 combinações de tratamento seriam randomizadas para unidades experimentais, evitando assim vieses de fontes desconhecidas.

Em problemas de pesquisa mais complicados, muitos fatores de tratamento podem estar em estudo. A eficiência dos planejamentos fatoriais pode ser demonstrada mesmo no caso mais simples, onde cada fator tem apenas dois níveis. Por exemplo, considere um experimento com quatro fatores. Um planejamento fatorial exigiria todas as combinações de quatro fatores em dois níveis ou \(2^4 = 16\) células. Se duas réplicas fossem executadas para cada célula, haveria um total de \(2\times 16 = 32\) experimentos ou execuções.

Para examinar o efeito de qualquer um dos quatro fatores, metade das execuções ou \(2\times 2^3 = 16\) devidas à replicação oculta, estariam em um nível do fator e a outra metade no outro nível. Assim, o efeito do tratamento consistiria em uma diferença de duas médias de 16.

Os resultados dos mesmos 32 experimentos podem ser usados para calcular o efeito do tratamento para cada um dos quatro fatores. Para ter igual precisão na comparação dos efeitos do tratamento usando um plano de cada vez, seriam necessárias 32 execuções para comparar os níveis de cada fator, mantendo os outros constantes. Isso resultaria em \(4\times 16 + 16 = 80\) experimentos ou 2.5 vezes o número necessário para um planejamento fatorial!


3.3 Interpretando interações


Se existir uma interacção ou efeito conjunto entre dois factores, o efeito de um factor sobre a resposta será diferente dependendo do nível do outro factor. Isto pode ser ilustrado graficamente na Figura 3.2. No lado esquerdo da figura há um gráfico de contorno que representa os resultados de uma série de experimentos com helicópteros de papel.

Figura 3.2: Gráfico de contorno do tempo de voo para experimento com helicópteros de papel.

Este gráfico pode ser interpretado como um mapa topológico com as linhas representando contornos de tempos de voo iguais. Você pode simular o que aconteceria em uma série de experimentos onde o comprimento da asa fosse mantido constante e a largura do corpo variasse, desenhando uma linha reta paralela ao eixo da largura do corpo através do gráfico de contorno.

O tempo de vôo para várias execuções pode ser lido como o rótulo das curvas de nível que a linha reta cruza. Por exemplo, se o comprimento da asa fosse mantido constante em um valor abaixo do seu ponto médio à esquerda do gráfico de contorno, os tempos de voo resultantes de cinco corridas com larguras de corpo variadas são representados como a linha preta traçada no gráfico à direita da Figura 3.2.

Se o comprimento da asa fosse mantido constante em um valor mais alto, a linha cinza indica como seria o resultado de uma série de experimentos com a largura do corpo. O fato de as duas linhas ou curvas do lado direito da figura não serem paralelas indica que há uma interação entre o comprimento da asa e a largura do corpo.

Elas mostram que o efeito da largura do corpo depende do comprimento da asa. As interações são comuns no mundo real, mas usar a estratégia clássica de experimentação uma de cada vez pressupõe tacitamente que as interações não existem. Para ver a falácia que resulta desta suposição, examine a Figura 3.3, que representa o que aconteceria se procurássemos a combinação ideal de comprimento de asa e largura de corpo.

Figura 3.3: Otimização uma de cada vez com helicópteros de papel.

O conjunto vertical de círculos é desenhado nas combinações de comprimento da asa e largura do corpo para uma série de experimentos que variam o comprimento da asa enquanto mantêm a largura do corpo constante. Os números dentro dos círculos representam os tempos de voo resultantes.

Depois de examinar o resultado desta série de experimentos, o comprimento ideal da asa seria escolhido e outra série de experimentos seria conduzida mantendo o comprimento da asa constante em seu valor ideal e variando a largura do corpo. Os resultados desses experimentos podem ser visualizados como uma série horizontal de círculos. O resultado máximo 2.8, não é o ideal geral, porque o comprimento ideal da asa depende da largura do corpo e vice-versa.

Quando o efeito dos fatores é próximo de linear a interação é mais fácil de explicar em palavras. A Tabela 3.1 mostra os resultados de um experimento fatorial conduzido por Derringer (1974) para determinar o efeito de compostos elastômeros na viscosidade da sílica a 100℃.

\[ \begin{array}{c|c}\hline \mbox{Óleo de nafteno (phr)} & \mbox{Enchimento (phr)} \\ \begin{array}{c} \\ 0 \\ 10 \\ 20 \\ 30 \end{array} & \begin{array}{cccccc} 0 & 12 & 24 & 36 & 48 & 60 \\ \hline 25 & 30 & 35 & 40 & 50 & 60 \\ 18 & 21 & 24 & 28 & 33 & 41 \\ 13 & 15 & 17 & 20 & 24 & 29 \\ 11 & 14 & 15 & 17 & 18 & 25 \\ \end{array} \\\hline \end{array} \] Tabela 3.1: Viscosidade Mooney da sílica B a 100℃. Um viscosímetro Mooney é um instrumento usado para medir a viscosidade Mooney de borrachas, foi inventado por Melvin Mooney.

Os compostos elastômeros foram óleo de Nafteno (Naphthene Oil), estudado em 4 níveis e conteúdo de carga (Filler), estudado em 6 níveis. A Figura 3.4 mostra uma representação gráfica dos dados da tabela. Essa figura é chamada de gráfico de interação. À medida que o conteúdo (Filler) aumenta de 0 a 60, a viscosidade aumenta ao longo de uma tendência bastante linear. No entanto, a inclinação da linha de tendência depende do nível do óleo de nafteno. Quando não há adição de óleo nafteno, aumentar o conteúdo de 0 para 60 faz com que a viscosidade aumente rapidamente de 25 para 60; mas quando há 30 phr de óleo de nafteno, aumentar o conteúdo de 0 para 60 causa um aumento mais gradual na viscosidade de 11 para 25.

Figura 3.4: Gráfico de interação de enchimento e óleo de nafteno.

Como as interações são comuns em experimentos fatoriais é importante aprender como explicar ou interpretar uma interação para apresentar claramente os resultados dos estudos de pesquisa. A melhor maneira de fazer isso é descrever o efeito de um fator sobre a resposta e, em seguida, contrastar ou comparar como esse efeito muda dependendo do nível do outro fator. Um gráfico de interação, como a Figura 3.4, pode orientar a descrição ou interpretação. Muitos outros exemplos de interpretações de interações serão dados ao longo deste capítulo e no restante do texto.


3.4 Criando um plano fatorial de dois fatores no R


Um experimento fatorial pode ser facilmente criado usando R de diversas maneiras. Por exemplo, laços aninhados podem ser usados, a função base R expand.grid ou várias funções disponíveis em pacotes desenvolvidos por usuários. Por exemplo, a função expand.grid, que cria um quadro ou conjunto de dados contendo todas as combinações possíveis de fatores fornecidos é ilustrada abaixo para criar um experimento fatorial para estudar helicópteros de papel.

D = expand.grid( BW = as.factor(c(3.25, 3.75, 4.25)), WL = as.factor(c(4, 5, 6)) )
D
##     BW WL
## 1 3.25  4
## 2 3.75  4
## 3 4.25  4
## 4 3.25  5
## 5 3.75  5
## 6 4.25  5
## 7 3.25  6
## 8 3.75  6
## 9 4.25  6


Como pode ser visto, este código cria um fatorial \(3^2\) não replicado nos fatores Largaura do corpo (Body width) = BW e Comprimento da asa (Wing length) = WL com os níveis fornecidos para esses fatores.

Este experimento é armazenado no conjunto de dados D. Para replicar cada execução no experimento, a função R rbind, que empilha uma cópia do experimento fatorial \(3^2\) sobre a outra é usada conforme mostrado abaixo.

D = rbind(D, D)


Para randomizar a ordem das execuções, a função sample pode ser usada para criar uma ordem aleatória dos números de execução de 1 a 18. Em seguida, as linhas no conjunto de dados D são ordenadas por esta lista de ordem aleatória.

Finalmente, as colunas de fator do conjunto de dados D podem ser copiadas para um novo quadro de dados Copterdes. Isso é ilustrado no código na próxima página.

set.seed(2591)
D = D[order(sample(1:18)), ]
CopterDes = D[ c( "BW", "WL" )]
CopterDes
##      BW WL
## 11 3.75  4
## 13 3.25  5
## 4  3.25  5
## 10 3.25  4
## 16 3.25  6
## 6  4.25  5
## 5  3.75  5
## 7  3.25  6
## 1  3.25  4
## 17 3.75  6
## 18 4.25  6
## 15 4.25  5
## 3  4.25  4
## 2  3.75  4
## 12 4.25  4
## 14 3.75  5
## 9  4.25  6
## 8  3.75  6
# Alternate way of creating the design with fac.des;
# using this function the factor levels are indicators
# BW: 1=3.25, 2=3.75, and 3=4.25. WL: 1=4, 2=5, and 3=6.
library(DoE.base)
D = fac.design(nlevels=c(3,3),factor.names=c("BW","WL"),replications=2,seed=5)
D
##    run.no run.no.std.rp BW WL Blocks
## 1       1           2.1  2  1     .1
## 2       2           3.1  3  1     .1
## 3       3           1.1  1  1     .1
## 4       4           8.1  2  3     .1
## 5       5           7.1  1  3     .1
## 6       6           5.1  2  2     .1
## 7       7           4.1  1  2     .1
## 8       8           9.1  3  3     .1
## 9       9           6.1  3  2     .1
## 10     10           7.2  1  3     .2
## 11     11           3.2  3  1     .2
## 12     12           8.2  2  3     .2
## 13     13           6.2  3  2     .2
## 14     14           2.2  2  1     .2
## 15     15           9.2  3  3     .2
## 16     16           1.2  1  1     .2
## 17     17           5.2  2  2     .2
## 18     18           4.2  1  2     .2
## class=design, type= full factorial 
## NOTE: columns run.no and run.no.std.rp  are annotation, 
##  not part of the data frame


A lista produzida por este código mostra que o primeiro experimento consistiria em construir um helicóptero com largura de corpo de 3.75” e comprimento de asa de 4”, largando-o de uma altura fixa e cronometrando seu voo. O segundo experimento consiste em construir um helicóptero com largura de corpo de 3.75” e comprimento de asa de 5”, largando-o da mesma altura fixa, cronometrando seu voo e assim por diante.

A randomização ajudará a evitar preconceitos de quaisquer variáveis ocultas, como mudanças nas correntes de ar, mudanças na temperatura ou curvas de aprendizagem no lançamento ou cronometragem de helicópteros. Ao remover a instrução set.seed(2591) no código acima, uma ordem aleatória diferente dos experimentos resultará cada vez que o código for executado.


3.5 Análise de um delineamento fatorial de dois fatores no R


O modelo matemático para um planejamento fatorial de dois fatores completamente aleatório pode ser escrito como: \[\begin{equation*} \tag{3.1} y_{ijk} = \mu_{ij}+\epsilon_{ijk}, \end{equation*}\] onde \(i\) representa o nível do primeiro fator, \(j\) representa o nível do segundo fator e \(k\) representa o número réplicas. Este modelo é chamado de modelo de média de célula e \(\mu_{ij}\) representa a resposta esperada na \(ij\)-ésima célula.

Outra forma de representar o modelo é o modelo de efeitos \[\begin{equation*} \tag{3.2} y_{ijk} = \mu+\alpha_i+\beta_j+\alpha\beta_{ij}+\epsilon_{ijk}\cdot \end{equation*}\]

Neste modelo, \(\alpha_i\), \(\beta_j\) são os efeitos principais e representam a diferença entre a média marginal de todos os experimentos no \(i\)-ésimo nível do primeiro fator e a média geral e a diferença entre a média marginal no \(j\)-ésimo nível do segundo fator e a média geral, respectivamente.

Os efeitos de interação, \(\alpha\beta_{ij}\), representam a diferença entre a média da célula, \(\mu_{ij}\) e \(\mu+\alpha_i + \beta_j\). Com essas definições, \(\sum_i\alpha_i = 0\), \(\sum_j \beta_j = 0\), \(\sum_i \alpha\beta_{ij} = 0\) e \(\sum_j\alpha\beta_{ij} = 0\).

As suposições usuais são que os erros experimentais são independentes e \(\epsilon_{ijk}\sim N(0,\sigma^2)\). A suposição de independência é garantida se as combinações de tratamento forem atribuídas aleatoriamente às unidades experimentais, as suposições de igual variância e normalidade podem ser verificadas com um gráfico dos resíduos versus previstos e um gráfico de probabilidade normal dos resíduos, conforme descrito na Seção 2.4.

As funções estimáveis, conforme descritas para o modelo de um fator na Seção 2.4.4, são combinações lineares das médias das células que podem ser expressas como uma combinação linear dos dados. Para o delineamento fatorial de dois fatores, as médias das células, \[ \mu_{ij} = \mu + \alpha_i + \beta_j + \alpha\beta_{ij}, \] são funções estimáveis, mas os efeitos individuais, \(\alpha_i\), \(\beta_j\) e \(\alpha\beta_{ij}\), não são funções estimáveis.

Contrastes entre efeitos como \(\sum_i c_i\alpha_i\) e \(\sum_j c_j\beta_j\), onde \(\sum_i c_i = 0\) e \(\sum_j c_j = 0\) são estimáveis apenas no modelo aditivo onde todos os \(\alpha\beta_{ij}\) são zero. Contrastes da forma \(\sum_i\sum_j b_{ij}\alpha\beta_{ij}\), onde \(\sum_i b_{ij} = 0\) e \(\sum_j b_{ij} = 0\) são estimáveis mesmo no modelo não aditivo.

As funções estimáveis e seus erros padrão podem ser calculados com a função estimable no pacote R gmodels. As médias marginais \(\mu+\alpha_i + \overline{\alpha\beta}_{i\bullet}\) e \(\mu+\beta_j + \overline{\alpha\beta}_{\bullet j}\) são funções estimáveis e elas e as médias das células podem ser calculadas usando a função R model.tables, que será ilustrado na Seção 5.4.1.


3.5.1 Representação matricial do modelo e análise


O modelo de efeitos pode ser representado em notação matricial como \[\begin{equation*} \tag{3.3} \pmb{y}=\pmb{X}\beta+\pmb{\epsilon}=\begin{pmatrix} 1 \; | \; X_A \; | \; X_B \; | \; X_{AB} \end{pmatrix}\begin{pmatrix} \mu \\ \beta_A \\ \beta_B \\ \beta_{AB} \end{pmatrix}+\pmb{\epsilon}\cdot \end{equation*}\]

Por exemplo, considere um caso em que o primeiro fator possui dois níveis, o segundo fator possui três níveis e há duas réplicas por célula. Então \[ \begin{pmatrix} y_{111} \\ y_{112} \\ y_{211} \\ y_{212} \\ y_{121} \\ y_{122} \\ y_{221} \\ y_{222} \\ y_{131} \\ y_{132} \\ y_{231} \\ y_{232} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \alpha\beta_{11} \\ \alpha\beta_{21} \\ \alpha\beta_{12} \\ \alpha\beta_{22} \\ \alpha\beta_{13} \\ \alpha\beta_{23} \end{pmatrix} + \begin{pmatrix} \epsilon_{111} \\ \epsilon_{112} \\ \epsilon_{211} \\ \epsilon_{212} \\ \epsilon_{121} \\ \epsilon_{122} \\ \epsilon_{221} \\ \epsilon_{222} \\ \epsilon_{131} \\ \epsilon_{132} \\ \epsilon_{231} \\ \epsilon_{232} \end{pmatrix}\cdot \]

O \(\pmb{X}^\top\pmb{X}\) é singular e para resolver as equações normais, usando a codificação de tratamento padrão, a função R lm elimina os indicadores para o primeiro nível de cada fator nas colunas de efeito principal e cria as colunas para a interação como todos os possíveis produtos das principais colunas de efeito.

Isto torna a matriz \(\pmb{X}^\top\pmb{X}\) de posto completo, como foi o caso do modelo de um factor nas Secções 2.3.1 e 2.3.2. Para o modelo de dois fatores onde o primeiro fator tem dois níveis, o segundo fator tem três níveis e há duas réplicas por célula, o vetor \(\pmb{y}\) e a matriz \(\pmb{X}\) recodificada seriam como mostrados abaixo. \[ \begin{pmatrix} y_{111} \\ y_{112} \\ y_{211} \\ y_{212} \\ y_{121} \\ y_{122} \\ y_{221} \\ y_{222} \\ y_{131} \\ y_{132} \\ y_{231} \\ y_{232} \end{pmatrix}, \qquad \pmb{X} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 0 & 1 \end{pmatrix}\cdot \]

Esta codificação de tratamento faz com que a célula \(\mu_{11}\) signifique o padrão, e as estimativas de efeito resultantes são mostradas abaixo. \[ \big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{X}^\top \pmb{y} = \widehat{\beta}=\begin{pmatrix} \widehat{\mu}+\widehat{\alpha}_1 +\widehat{\beta}_1 +\widehat{\alpha\beta}_{11} \\ \widehat{\alpha}_2-\widehat{\alpha}_1 +\widehat{\alpha\beta}_{21} - \widehat{\alpha\beta}_{11} \\ \widehat{\beta}_2-\widehat{\beta}_1 +\widehat{\alpha\beta}_{12} -\widehat{\alpha\beta}_{11} \\ \widehat{\beta}_3-\widehat{\beta}_1 +\widehat{\alpha\beta}_{13} -\widehat{\alpha\beta}_{11} \\ \widehat{\alpha\beta}_{11}+\widehat{\alpha\beta}_{22} -\widehat{\alpha\beta}_{12} -\widehat{\alpha\beta}_{21} \\ \widehat{\alpha\beta}_{11}+\widehat{\alpha\beta}_{23} -\widehat{\alpha\beta}_{13} -\widehat{\alpha\beta}_{21} \end{pmatrix}\cdot \]

A soma do erro dos quadrados \[ ssE = \pmb{y}^\top\pmb{y} − \widehat{\beta}^\top \pmb{X}^\top \pmb{y} = \pmb{y}^\top \big(\pmb{I} −\pmb{X}(\pmb{X}^\top \pmb{X})^{-1} X^\top\big)\pmb{y}, \] onde \(\widehat{\beta}=\big(\pmb{X}^\top\pmb{X}\big)^{-1}\pmb{X}^\top \pmb{y}\) são as estimativas produzidas pela função lm em R. Para testar as hipóteses \[ H_0 \, : \, \alpha_1 = \alpha_2 = 0, \qquad H_0 \, : \, \beta_1 =\beta_2 =\beta_3 = 0 \] e \[ H_0 \, : \, \alpha\beta_{11} = \alpha\beta_{21} = \alpha\beta_{12} = \alpha\beta_{22} = \alpha\beta_{13} = \alpha\beta_{23} = 0, \] os testes \(F\) de razão de verossimilhança são obtidos calculando proporções dos quadrados médios da ANOVA.

O que o lm designa como soma dos quadrados do fator A é \[ ssA = \widehat{\beta}^\top \pmb{X}^\top\pmb{y} - \big(\pmb{1}^\top \pmb{y}\big)^2/(\pmb{1}^\top\pmb{1}), \] onde o modelo é simplificado para incluir apenas os efeitos do primeiro fator, ou seja, \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_A \end{pmatrix}\). A soma dos quadrados dos erros para este modelo simplificado é \(ssE_A\).

A soma dos quadrados do fator A é denotada por \(R(\alpha \, | \, \mu)\). As somas dos quadrados para o fator B são denotadas como \(R(\beta \, | \, \mu )=ssE_A−ssE_B\) onde \(ssE_B\) é a soma dos erros dos quadrados do modelo reduzido onde \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_A \; | \; X_B \end{pmatrix}\). Finalmente, a soma dos quadrados para a interação AB é denotada \(R( \alpha\beta \, | \, \beta,\alpha,\mu )=ssE_B − ssE\).

Em geral, quando há níveis \(a\) do fator A, níveis \(b\) do fator B e \(r\) réplicas por célula, a tabela ANOVA para o experimento fatorial de dois fatores pode ser representada simbolicamente como mostrado na Tabela 3.2.

\[ \begin{array}{lcccc}\hline \mbox{Fonte} & \mbox{df} & \mbox{Soma de Quadrados} & \mbox{Quadrados Médios} & \mbox{Estatística } F \\\hline A & a-1 & R(\alpha \, | \, \mu) & ssA/(a-1) & F = msA/msE \\ B & b-1 & R(\beta \, | \, \mu ) & ssB/(b-1) & F = msB/msE \\ AB & (a-1)(b-1) & R( \alpha\beta \, | \, \beta,\alpha,\mu ) & ssAB/(a-1)(b-1) & F= msAB/msE\\ \hline \mbox{Erro} & ab(r-1) & ssE & ssE/ab(r-1) & \\ \hline \end{array} \] Tabela 3.2: Tabela de Análise de Variância.

As somas dos quadrados \(ssA\), \(ssB\) e \(ssAB\) também podem ser escritas na forma \[ \begin{array}{rcl} ssA & = & \big(\pmb{L}_\alpha\widehat{\beta} \big)^\top\big(\pmb{L}_\alpha(\pmb{X}^\top\pmb{X})^{-1} \pmb{L}_\alpha^\top\big)^{-1} \big(\pmb{L}_\alpha\widehat{\beta} \big)\\ ssB & = & \big(\pmb{L}_\beta\widehat{\beta} \big)^\top\big(\pmb{L}_\beta(\pmb{X}^\top\pmb{X})^{-1} \pmb{L}_\beta^\top\big)^{-1} \big(\pmb{L}_\beta\widehat{\beta} \big)\\ ssAB & = & \big(\pmb{L}_{\alpha\beta}\widehat{\beta} \big)^\top\big(\pmb{L}_{\alpha\beta}(\pmb{X}^\top\pmb{X})^{-1} \pmb{L}_{\alpha\beta}^\top\big)^{-1} \big(\pmb{L}_{\alpha\beta}\widehat{\beta} \big)\\ \end{array} \] onde \(\pmb{L}_\alpha\), \(\pmb{L}_\beta\) e \(\pmb{L}_{\alpha\beta}\) são matrizes de contraste calculadas internamente pela função lm, veja Goodnight (1980). Sob a hipótese nula, as razões \(F\) \(msA/msE\), \(msB/msE\) e \(msAB/msE\) seguem a distribuição \(F\) com os graus de liberdade mostrados na tabela e, na alternativa, seguem a distribuição \(F\) não central.

O parâmetro de não centralidade para \(F = msA/msE\) é dado pela expressão \[ \lambda_\alpha = \big(1/\sigma^2\big)\big(\pmb{L}_\alpha\widehat{\beta} \big)^\top\big(\pmb{L}_\alpha(\pmb{X}^\top\pmb{X})^{-1} \pmb{L}_\alpha^\top\big)^{-1} \big(\pmb{L}_\alpha\widehat{\beta} \big)\cdot \] Os parâmetros de não centralidade para as razões \(F\) \(msB/msE\) e \(msAB/msE\) são dados de forma semelhante. Quando há um número igual \(r\) de réplicas em cada célula, pode-se mostrar que os parâmetros de não centralidade são iguais a \[\begin{equation*} \tag{3.4} \lambda_\alpha = b \, r \sum_i \alpha_i^2/\sigma^2, \end{equation*}\] \[\begin{equation*} \tag{3.5} \lambda_\beta = a \, r \sum_j \beta_j^2/\sigma^2 \end{equation*}\] e \[\begin{equation*} \tag{3.6} \lambda_{\alpha\beta} = r \sum_i\sum_j \alpha\beta_{ij}^2/\sigma^2\cdot \end{equation*}\]

Para ilustrar a análise de um experimento fatorial de dois fatores usando a função R aov considere os dados da Tabela 3.3. Estes são os resultados de um experimento de dois fatores apresentado por Hunter (1983). Nestes dados, um experimento consistiu em queimar uma quantidade de combustível e determinar as emissões de CO liberadas.

A unidade experimental é a porção de um combustível padrão necessária para uma corrida e a resposta \(y\) é a concentração de emissões de monóxido de carbono (CO) em \(\mbox{gramas/metro}^3\) determinada a partir dessa corrida. Houve duas repetições para cada combinação de níveis de fator separados por vírgulas na Tabela 3.3. O fator A é a quantidade de etanol adicionada a uma unidade experimental ou porção do combustível padrão e o fator B é a relação combustível/ar usada durante a queima desse combustível.

\[ \begin{array}{ccc}\hline A = \mbox{adição de etanol} & B=\mbox{relação ar/combustível} & y=\mbox{emissões de CO} \\\hline 0.1 & 14 & 66, \, 62 \\ 0.1 & 15 & 72, \, 67 \\ 0.1 & 16 & 68, \, 66 \\ 0.2 & 14 & 78, \, 81 \\ 0.2 & 15 & 80, \, 81 \\ 0.2 & 16 & 66, \, 69 \\ 0.3 & 14 & 90, \, 94 \\ 0.3 & 15 & 75, \, 78 \\ 0.3 & 16 & 60, \, 58 \\\hline \end{array} \] Tabela 3.3: Dados do experimento do etanol combustível.

Os dados para este experimento são armazenados no conjunto ou quadro de dados COdata no pacote daewr onde os níveis de etanol e proporção são armazenados como os fatores Eth e Ratio.

Os comandos R para analisar os dados são mostrados a seguir.

library(daewr)
mod1 = aov( CO ~ Eth * Ratio, data = COdata )
summary(mod1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Eth          2  324.0   162.0   31.36 8.79e-05 ***
## Ratio        2  652.0   326.0   63.10 5.07e-06 ***
## Eth:Ratio    4  678.0   169.5   32.81 2.24e-05 ***
## Residuals    9   46.5     5.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Na tabela ANOVA resultante, mostrada acima, pode-se ver que a função aov produz uma tabela de somas de quadrados, conforme descrito anteriormente. Pode-se observar nas tabelas que os dois efeitos e a sua interação são significativos, conforme indicado pelos \(p\)-valores à direita dos valores \(F\).

A função model.tables produz os resultados mostrados a continuação. A linha superior é a estimativa da média geral \(\widehat{\mu}\). As próximas duas seções mostram as médias marginais de cada fator juntamente com o desvio padrão dos valores médios de cada média. Se a interação não fosse significativa, as médias marginais revelariam a direção dos efeitos dos fatores, mas outras comparações pré-planejadas ou outros procedimentos de comparação múltipla poderiam ser usados para tirar conclusões definitivas. A próxima seção mostra as médias das células e a seção final mostra os erros padrão das diferenças nas médias marginais e das células.

model.tables( mod1, type = "means", se = T )
## Tables of means
## Grand mean
##          
## 72.83333 
## 
##  Eth 
## Eth
##   0.1   0.2   0.3 
## 66.83 75.83 75.83 
## 
##  Ratio 
## Ratio
##   14   15   16 
## 78.5 75.5 64.5 
## 
##  Eth:Ratio 
##      Ratio
## Eth   14   15   16  
##   0.1 64.0 69.5 67.0
##   0.2 79.5 80.5 67.5
##   0.3 92.0 76.5 59.0
## 
## Standard errors for differences of means
##           Eth Ratio Eth:Ratio
##         1.312 1.312     2.273
## replic.     6     6         2


Duas estimativas de contrastes específicos dos efeitos principais, a função estimable no pacote R gmodels podem ser utilizadas. Para usá-lo, devemos primeiro construir contrastes para substituir os contrastes de tratamento padrão usados pela função R aov.

Por exemplo, na primeira afirmação abaixo construímos os coeficientes de contrastes para comparar o primeiro nível do fator com o terceiro em um fator de três níveis. Um segundo contraste ortogonal ao primeiro também é construído e a matriz de contraste cm é criada usando os dois contrastes como colunas.

c1 = c(-1/2, 0, 1/2)
c2 = c(0.5, -1, 0.5)
cm = cbind( c1, c2 )
cm
##        c1   c2
## [1,] -0.5  0.5
## [2,]  0.0 -1.0
## [3,]  0.5  0.5


Na chamada da função aov abaixo, a matriz de contrastes cm será usada para ambos os efeitos principais, em vez dos contrastes de tratamento padrão usados pela função aov. As próximas linhas carregam o pacote de funções gmodels e criam rótulos para os contrastes que comparam o primeiro nível de fator com o terceiro nível de fator.

O vetor após cada rótulo é um vetor indicador para o qual o coeficiente do modelo é exibido. Seleciona o primeiro coeficiente para etanol (Eth) e proporção (Ratio). Finalmente, a função estimable é chamada com as entradas sendo o modelo em mod2, que foi criado pela função aov e os rótulos e definições de contraste em c. Os resultados são mostrados abaixo.

mod2 = aov( CO ~ Eth * Ratio, contrasts = list( Eth = cm, Ratio = cm ), data = COdata)
library(gmodels)
c = rbind( 'Ethanol 0.3 vs 0.1' = c(0,1,0,0,0,0,0,0,0), 
           'Ratio 16 vs 14' = c(0,0,0,1,0,0,0,0,0) )
estimable(mod2,c)
##                    Estimate Std. Error    t value DF     Pr(>|t|)
## Ethanol 0.3 vs 0.1        9   1.312335   6.858007  9 7.406588e-05
## Ratio 16 vs 14          -14   1.312335 -10.668011  9 2.083651e-06


Ambas são funções estimáveis e as estimativas junto com seus respectivos erros padrão e razões ou valor da estatística \(t\) para testar as hipóteses, \[ H_0 \, : \, \sum_i c_i\alpha_i = 0 \qquad \mbox{e} \qquad H_0 \, : \, \sum_j c_j \beta_j =0, \] também foram mostrados.

Essas estimativas seriam estimáveis e significativas se não houvesse interação significativa entre o nível de adição de etanol e a relação ar/combustível, mas neste caso há uma interação significativa e a diferença nas emissões de CO causada pela alteração da quantidade de adição de etanol dependerá na relação ar/combustível e a diferença na emissão de CO causada pela alteração da relação ar/combustível dependerá da quantidade de etanol adicionado.

Um gráfico de interação é a melhor maneira de interpretar esses resultados. Um gráfico de interação pode ser gerado usando a função R interaction.plot conforme mostrado abaixo. Este código usa a instrução with para chamar a função interaction.plot usando nomes de variáveis no quadro de dados COdata para produzir a Figura 3.5.

par(mar = c(4,4,1,1))
with(COdata, interaction.plot(Eth, Ratio, CO, type = "b", pch = c(18,24,22), leg.bty = "o", 
                               main = "Interaction Plot of Ethanol and air/fuel ratio", 
                               xlab = "Ethanol",ylab = "CO emissions"))
grid()

Figura 3.5: Gráfico da interação entre o etanol e a relação ar/combustível.

Neste gráfico podemos ver mais claramente a dependência dos efeitos. Aumentar a quantidade de etanol adicionado ao combustível de 0.1 para 0.3 faz com que as emissões de CO aumentem linearmente de 64 gramas/litro para 92 gramas/litro quando a relação ar/combustível está em seu nível baixo de 14. Isso é mostrado pela linha pontilhada com losangos pretos representando as médias das células.

No entanto, quando a relação ar/combustível está em seu alto nível de 16 (ilustrado pela linha sólida com quadrados representando as médias das células), aumentar o etanol adicionado ao combustível de 0.1 para 0.3 na verdade causa uma diminuição nas emissões de CO de 67 gramas/litro para 59 gramas/litro ao longo de uma tendência quase linear. Finalmente, quando a relação ar/combustível é mantida constante em seu nível médio de 15 (ilustrado pela linha tracejada com triângulos representando as médias das células), aumentar o etanol de 0.1 para 0.2 faz com que as emissões de CO aumentem em 11 gramas/litro; mas um aumento adicional no etanol para 0.3 provoca uma diminuição nas emissões de CO de 4 gramas/litro para 76.5.

A interpretação acima ilustra novamente o princípio de comparação do efeito de um factor através dos níveis do outro factor, a fim de descrever uma interacção. Isto foi feito comparando o efeito da alteração da adição de etanol entre os níveis de relação ar/combustível. Também poderia ser feito de maneira oposta. Por exemplo, o código R abaixo inverte o gráfico de interação conforme mostrado na Figura 3.6.

par(mar = c(4,4,1,1))
with(COdata, interaction.plot(Ratio, Eth, CO, type = "b", pch = c(18,24,22), leg.bty = "o", 
                              main="Interaction Plot of Ethanol and air/fuel ratio", 
                              xlab = "Ratio", ylab = "CO emissions"))
grid()

Figura 3.6: Gráfico da interação o etanol e a relação ar/combustível.

Neste gráfico, a linha sólida, com quadrados representando as médias das células, mostra o efeito do aumento da relação ar/combustível quando o etanol é adicionado na alta taxa de 0.3. As emissões de monóxido de carbono diminuem linearmente de 92 gramas/litro para 59 gramas/litro.

No entanto, quando o etanol é adicionado a uma taxa baixa de 0.1, as emissões de CO na verdade aumentam ligeiramente de 64 gramas/litro para 67 gramas/litro como resultado do aumento da relação ar/combustível de 14 para 16. Isto pode ser visto no gráfico pontilhado. linha com losangos pretos representando as médias das células.

Quando o etanol é adicionado na taxa média de 0.2, há pouca mudança nas emissões de CO quando a relação ar/combustível é aumentada de 14 para 15, mas há uma diminuição nas emissões de CO de 13 gramas/litro causada pelo aumento da relação ar/combustível. proporção de 15 a 16. O último resultado pode ser visualizado na linha tracejada com triângulos representando as médias das células.

Qualquer forma de apresentar e interpretar a interação é válida, desde que se discuta como o efeito de um fator muda dependendo do nível do outro. Os efeitos fatoriais que devem ser comparados dependem de qual deles é de maior interesse em um determinado problema de pesquisa. Outra coisa a notar sobre as duas interpretações é que são assumidas relações de causa e efeito. Dizemos que a mudança na resposta é causada pela mudança no fator ou que a mudança na resposta é o resultado da mudança no fator. Esta afirmação não pôde ser feita ao discutir os resultados de um estudo observacional.


3.5.2 Determinando o número de réplicas


Um dos dois métodos possíveis pode ser seguido para determinar o número de repetições para um experimento fatorial que resultará em um poder entre 0.80 e 0.90, para detectar diferenças que tenham significância prática. O primeiro método é considerar a detecção de diferenças entre as médias das células. O segundo método é considerar a detecção de uma diferença prática de tamanho nas médias marginais dos fatores do experimento.

Ao procurar diferenças entre as médias das células, as células no fatorial são consideradas um grupo não estruturado, como em um experimento de um fator. Usando o modelo de médias celulares \(y_{ijk} = \mu_{ij} +\epsilon_{ijk}\) o procedimento é o mesmo descrito para o modelo de um fator \(y_{ij} = \mu_i + \epsilon_{ij}\) na Seção 2.6.

O parâmetro de não centralidade para o teste \(F\) é: \[ \lambda=\left(\dfrac{r}{\sigma^2}\right)\sum_{i=1}^a\sum_{j=1}^b \big(\overline{\mu}_{ij}-\overline{\mu}_{\bullet\bullet}\big)^2\cdot \] Ao procurar diferenças nas médias marginais dos fatores, o fator de não centralidade para o primeiro efeito principal é: \[ \lambda_a = b \, r\sum_i \dfrac{\alpha_i^2}{\sigma^2} = b \, r \dfrac{1}{\sigma^2}\sum_i \big(\overline{\mu}_{i\bullet}-\overline{\mu}_{\bullet\bullet}\big)^2, \] e para o segundo efeito principal o fator de não centralidade é: \[ \lambda_b = a \, r\sum_j \dfrac{\beta_j^2}{\sigma^2} = a \, r \dfrac{1}{\sigma^2} \sum_j \big(\overline{\mu}_{\bullet j}-\overline{\mu}_{\bullet\bullet}\big)^2\cdot \]

Se \(\Delta\) for considerado o tamanho de uma diferença prática nas médias das células, então o menor \(\lambda\) poderia ser com duas células diferindo em pelo menos \(\Delta\) é \(r\Delta^2/2\sigma^2\). Da mesma forma, se \(\Delta\) for considerado o tamanho de uma diferença prática nas médias marginais para o fator A, o menor \(\lambda_a\) poderia ser com duas médias marginais diferindo em pelo menos \(\Delta\) é \(b\, r \Delta^2/2\sigma^2\).

Aqui, novamente, podemos ver a eficiência dos planos ou planejamentos fatoriais porque o fator de não centralidade para detectar diferenças nas médias marginais do fator A é maior do que o fator de não centralidade para detectar diferenças nas médias das células por um fator de \(b\), o número de níveis do fator B.

Considere o seguinte exemplo. Um experimento com helicópteros de papel está planejado para investigar os efeitos de quatro níveis do fator A = comprimento da asa (wing length) e quatro níveis do fator B = largura do corpo (body width), sobre o tempo de vôo. Se os experimentos piloto com nove réplicas de um projeto resultaram em tempos de voo de 2.8, 2.6, 3.5, 3.0, 3.1, 3.5, 3.2, 3.4 e 3.4 segundos. Quantas réplicas seriam necessárias para detectar uma diferença nos tempos de voo de 1 segundo com uma poder de 0.90?

A partir dos testes piloto o desvio padrão do erro experimental pode ser estimado como \(s = 0.32\). Se \(\Delta = 1.0\) for considerado uma diferença prática de tamanho nas médias das células, o código R na Seção 2.7 pode ser modificado para fornecer a resposta. Em um plano fatorial 4 por 4 há 16 células, então o número de níveis do fator é considerado 16. O código R modificado é mostrado abaixo.

library(daewr)
rmin = 2 # smallest number of replicates
rmax = 8 # largest number of replicates
sigma = .32
alpha = .05
Delta = 1
nlev = 16
nreps = c(rmin:rmax)
power = Fpower1(alpha, nlev, nreps, Delta, sigma)
options(digits = 5)
power
##      alpha nlev nreps Delta sigma   power
## [1,]  0.05   16     2     1  0.32 0.24173
## [2,]  0.05   16     3     1  0.32 0.48174
## [3,]  0.05   16     4     1  0.32 0.69246
## [4,]  0.05   16     5     1  0.32 0.83829
## [5,]  0.05   16     6     1  0.32 0.92326
## [6,]  0.05   16     7     1  0.32 0.96664
## [7,]  0.05   16     8     1  0.32 0.98655


Os resultados da execução deste código mostram que seriam necessárias 6 réplicas por célula para obter um poder de pelo menos 0.90.

Se \(\Delta = 1.0\) for considerada uma diferença prática de tamanho nas médias marginais para um dos fatores, os resultados serão diferentes. Os graus de liberdade para o numerador seriam \(\nu_1 = 4-1\), os graus de liberdade para o denominador seriam \(\nu_2 = 16(r-1)\), o fator de não centralidade para um efeito principal A seria \(\lambda_a= b \, r \Delta^2/2\sigma^2\) e o fator de não-entralidade para o efeito principal B seria \(\lambda_b = a\, r \Delta^2/2\sigma^2\).

O código R abaixo demonstra o uso da função Fpower2 no pacote daewr para determinar o número de réplicas necessárias para detectar uma diferença nas médias marginais dos fatores em um fatorial de dois fatores.

Os argumentos que devem ser fornecidos ao Fpower2 são: alfa, nlev (um vetor de comprimento 2 contendo o número de níveis do primeiro fator (A) e do segundo fator (B)), nreps=\(r\), Delta=\(\Delta\) e sigma = \(\sigma\).

library(daewr)
rmin = 2 # smallest number of replicates
rmax = 4 # largest number of replicates
alpha = .05
sigma = .32
Delta = 1.0
nlev = c(4,4)
nreps = c(rmin:rmax)
result = Fpower2(alpha, nlev, nreps, Delta, sigma)
options(digits = 5)
result
##      alpha a b nreps Delta sigma  powera  powerb
## [1,]  0.05 4 4     2     1  0.32 0.99838 0.99838
## [2,]  0.05 4 4     3     1  0.32 1.00000 1.00000
## [3,]  0.05 4 4     4     1  0.32 1.00000 1.00000


Pode ser visto que com apenas \(r\) = duas réplicas por célula o poder para detectar uma diferença de \(\Delta= 1.0\) nas médias marginais para o fator A ou B é maior do que o poder para detectar diferenças de \(\Delta= 1.0\) nas médias das células com \(r = 8\) réplicas por célula. Novamente, isso demonstra a eficiência dos experimentos fatoriais por meio da replicação oculta.

Com a capacidade de calcular o poder rapidamente é possível explorar muitos projetos potenciais antes de realmente executar os experimentos. O número de fatores, o número de níveis de cada fator e o número de réplicas em cada célula afetam o poder de detectar diferenças. Os cálculos do poder ajudam o experimentador a determinar o uso eficiente de seus recursos.


3.5.3 Análise com número desigual de réplicas por célula


Embora seja incomum planejar um experimento fatorial com um número desigual de réplicas por célula, os dados de um experimento fatorial podem acabar com um número desigual de réplicas devido a experimentos que não puderam ser concluídos, respostas que não puderam ser medidas ou simplesmente dados perdidos.

Contanto que a chance de perder uma observação não esteja relacionada aos níveis dos fatores de tratamento, os dados de um experimento fatorial com um número desigual de réplicas por célula ainda podem ser analisados e interpretados de maneira semelhante à forma como seria feito para o caso replicado igual. Porém, as fórmulas computacionais para análise dos dados diferem para o caso com número desigual de réplicas.

Para ilustrar por que a análise mostrada na Seção 3.5.1 é inadequada, considere novamente os dados do experimento com etanol combustível descrito na Seção 3.5.1. Desta vez, assuma uma observação na célula onde faltava a relação ar/combustível = 16 e o nível de etanol = 0.3. Então a Tabela 3.4 mostra os dados com cada valor de resposta escrito acima do seu valor simbólico esperado.

\[ \begin{array}{c|ccc}\hline & \mbox{ar/combustível} & \mbox{ar/combustível} & \mbox{ar/combustível} \\ \mbox{Etanol} & 14 & 15 & 16 \\\hline & 66 & 72 & 68 \\ 0.1 & 62 & 67 & 66 \\ & \mu+\alpha_1+\beta_1+\alpha\beta_{11} & \mu+\alpha_1+\beta_2+\alpha\beta_{12} & \mu+\alpha_1+\beta_3+\alpha\beta_{13}\\ \hline & 78 & 80 & 66 \\ 0.2 & 81 & 81 & 69 \\ & \mu+\alpha_2+\beta_1+\alpha\beta_{21} & \mu+\alpha_2+\beta_2+\alpha\beta_{22} & \mu+\alpha_2+\beta_3+\alpha\beta_{23} \\\hline & 90 & 75 & 60 \\ 0.3 & 94 & 78 & \\ & \mu+\alpha_3+\beta_1+\alpha\beta_{31} & \mu+\alpha_3+\beta_2+\alpha\beta_{32} & \mu+\alpha_3+\beta_3+\alpha\beta_{33} \\\hline \end{array} \] Tabela 3.4: Experimento de combustível com repetições desiguais.

O código R abaixo cria um quadro ou conjunto de dados contendo os dados da Tabela 3.4, inserindo um valor ausente na 18ª linha e na terceira coluna.

COdatam = COdata
COdatam[18, 3] = NA
tail(COdatam)
##    Eth Ratio CO
## 13 0.2    14 81
## 14 0.2    15 81
## 15 0.2    16 69
## 16 0.3    14 94
## 17 0.3    15 78
## 18 0.3    16 NA


A média da coluna marginal para os níveis do fator de relação ar/combustível calculado usando a instrução model.tables conforme mostrado na Seção 3.5.1 e o quadro de dados modificado COdatam seria 78.5, 75.5 e 65.8, respectivamente.

options(digits = 5)
(66+78+90+62+81+94)/6
## [1] 78.5
(72+80+75+67+81+78)/6
## [1] 75.5
(68+66+60+66+69)/5
## [1] 65.8


O valor esperado das médias marginais para as duas primeiras colunas seria: \(\mu+\beta_1\), \(\mu+\beta_2\), já que \((\alpha_1 +\alpha_2 + \alpha_3)/3 = 0\) e \((\alpha\beta_{1i}+\alpha\beta_{2i}+\alpha\beta_{3i})/3 = 0\) para \(i = 1,2\). No entanto, o valor esperado da média da última coluna marginal seria \[ \mu+\beta_3+(2\alpha_1+2\alpha_2+\alpha_3)/5+(2\alpha\beta_{13}+2\alpha\beta_{23}+\alpha\beta_{33})/5 \] e não é uma estimativa nao viciada de \(\mu+\beta_3\).

A comparação entre as médias da primeira e da terceira coluna não seria uma estimativa não viciada de \(\beta_1-\beta_3\). Da mesma forma, a média da última linha marginal não seria uma estimativa não viciada de \(\mu+\beta_3\).

Se a tabela ANOVA dos dados no arquivo COdatam for produzida com a função R aov, os testes \(F\) não testarão as mesmas hipóteses que testam no caso de igual número de réplicas nas células. Quando há um número desigual de réplicas nas células, o parâmetro de não centralidade para o teste \(F\) da hipótese \[ H_0 \, : \, \alpha_1 = \cdots = \alpha_a, \] que é baseado em \(R( \alpha \, | \, \mu )\) não será \(\lambda_a = r \, b \sum_i \alpha_i^2\) mas sim uma forma quadrática envolvendo os elementos de \(\alpha\), \(\beta\) assim como \(\alpha\beta\).

A não centralidade para o teste \(F\) de \[ H_0 \, : \, \beta_1 = \cdots =\beta_b \] com base em \(R(\beta \, | \, \mu,\alpha)\) será uma forma quadrática envolvendo os elementos de \(\beta\) e \(\alpha\beta\).

Para calcular somas de quadrados ajustadas para a hipótese nula para os efeitos principais, use a opção contr.sum na função R lm e a função Anova do pacote R car (Fox and Weisberg, 2011).

A opção type = “II” na função Anova calcula as somas dos quadrados do tipo II e a opção type = “III” produz as somas dos quadrados do tipo III. A soma dos quadrados do tipo II para os fatores A e B pode ser representada como \(ssA_{II} = R( \alpha \, | \, \mu,\beta)\) e \(ssB_{II} = R(\beta \, | \, \mu,\alpha)\), observando que \(R( \alpha \, | \, \mu,\beta)\) é a diferença nas somas dos quadrados dos erros para o modelo reduzido onde \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_B \end{pmatrix}\) e o modelo completo onde \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_A \; | \; X_B \; | \; X_{AB} \end{pmatrix}\).

O fator de não centralidade correspondente para o teste \(F\) correspondente será uma forma quadrática que envolve apenas \(\alpha^\top = (\alpha_1,\alpha_2,\alpha_3)\) e \(\alpha\beta\). Quando há um número igual de replicações por célula, as somas dos quadrados calculadas pela função aov são idênticas às somas dos quadrados do tipo II.

A soma dos quadrados do tipo III para os fatores A e B podem ser representadas como \(ssA_{III} = R(\alpha \, | \, \mu,\beta,\alpha\beta)\) e \(ssB_{III} = R(\beta \, | \, \mu,\alpha,\alpha\beta)\). \(R(\alpha \, | \, \mu,\beta,\alpha\beta)\) é a diferença nas somas dos quadrados dos erros para o modelo reduzido onde \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_B \; | \; X_{AB} \end{pmatrix}\) e o modelo completo onde \(\pmb{X} = \begin{pmatrix} 1 \; | \; X_A \; | \; X_B \; | \; X_{AB} \end{pmatrix}\).

O fator de não centralidade correspondente para o teste \(F\) correspondente será uma forma quadrática que envolve apenas \(\alpha^* = c(\alpha_1,\alpha_2,\alpha_3)\). Quando há um número igual de replicações por célula, as somas dos quadrados calculadas pela função aov são idênticas às somas dos quadrados do tipo III.

Alguns analistas preferem usar somas de quadrados do tipo II e outros preferem somas de quadrados do tipo III quando há um número desigual de réplicas por célula. Aqui ilustramos as somas dos quadrados do tipo III e os testes de hipóteses, embora as somas dos quadrados do tipo II possam ser obtidas alterando a opção de type = “III” para type = “II” na chamada à função Anova.

O código para produzir o experimento de etanol combustível com tabela ANOVA tipo III após a retirada da observação com o valor 58, da célula com relação ar/combustível = 16 e nível de etanol = 0.3 é mostrado abaixo.

library(car)
mod2 = lm( CO ~ Eth*Ratio, data = COdatam, 
           contrasts = list( Eth = contr.sum, Ratio = contr.sum ))
Anova( mod2, type="III" )
## Anova Table (Type III tests)
## 
## Response: CO
##             Sum Sq Df F value  Pr(>F)    
## (Intercept)  86198  1 15496.4 1.9e-14 ***
## Eth            319  2    28.7 0.00022 ***
## Ratio          511  2    46.0 4.1e-05 ***
## Eth:Ratio      555  4    24.9 0.00014 ***
## Residuals       44  8                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Para obter médias que tenham esperanças \(\mu+\beta_1\), \(\mu+\beta_2\) e \(\mu+\beta_3\) quando houver um número desigual de réplicas por célula, as médias ajustadas devem ser calculadas. As médias ajustadas, às vezes chamadas de médias de mínimos quadrados ou lsmeans para abreviar, são calculadas calculando as estimativas das médias marginais de células \((\widehat{\mu}+\widehat{\alpha}_1+\widehat{\beta}_j+\widehat{\alpha\beta}_{ij})\) obtidas a partir das estimativas de mínimos quadrados dos parámetros do modelo. Lembre-se de que as médias das células são funções estimáveis.

O código R abaixo calcula as estimativas das médias das células usando as estimativas de efeito do modelo mod2 criado pela função lm mostrada acima e, em seguida, calcula as médias marginais ajustadas ou médias marginais de mínimos quadrados para a relação ar/combustível usando a função R tapply.

p = data.frame( expand.grid( Eth = c(.1, .2, .3), Ratio = c(14,15,16) ) )
p[] = lapply(p, factor)
p = cbind( yhat = predict( mod2, p), p) 
with(p, tapply(yhat, Ratio, mean) )
##     14     15     16 
## 78.500 75.500 64.833


Nestes resultados pode-se observar que as médias das duas primeiras colunas são iguais à média aritmética simples das respostas nas duas primeiras colunas, conforme mostrado no começo desta seção, mas a média da terceira coluna é diferente e é uma estimativa mais precisa de \(\mu+\beta_3\). O pacote R lsmeans calcula automaticamente estes resultados e, além disso, calcula seus erros padrão e limites de confiança.

O código R abaixo ilustra o uso deste pacote para calcular as médias marginais ajustadas para etanol e relação ar/combustível.

library(lsmeans)
lsmeans(mod2, ~ Eth)
##  Eth lsmean    SE df lower.CL upper.CL
##  0.1   66.8 0.963  8     64.6     69.1
##  0.2   75.8 0.963  8     73.6     78.1
##  0.3   76.2 1.112  8     73.6     78.7
## 
## Results are averaged over the levels of: Ratio 
## Confidence level used: 0.95
lsmeans(mod2, ~ Ratio)
##  Ratio lsmean    SE df lower.CL upper.CL
##  14      78.5 0.963  8     76.3     80.7
##  15      75.5 0.963  8     73.3     77.7
##  16      64.8 1.112  8     62.3     67.4
## 
## Results are averaged over the levels of: Eth 
## Confidence level used: 0.95


Em geral, as somas de quadrados do tipo II ou III e lsmeans devem ser usadas, porque elas testarão as hipóteses corretas e fornecerão médias não viciadas dos níveis dos fatores se há um número igual ou desigual de réplicas por célula.


3.5.4 Teste de interação com uma réplica por célula


Quando há poder adequado para detectar efeitos principais com \(r = 1\) réplica por célula, faria sentido executar um planejamento fatorial com apenas uma observação por célula e um total de observações \(a\times b\). Adicionar uma réplica adicional a cada célula duplicaria o esforço e normalmente não seria necessário.

No entanto, com apenas uma réplica por célula em um planejamento fatorial, não há como calcular a ANOVA \(ssE\) e, portanto, não há como fazer testes \(F\) sobre os principais efeitos e interações da maneira tradicional. Se os termos de interações forem assumidos como zero, então os testes \(F\) sobre os efeitos principais podem ser feitos usando o modelo aditivo \(y_{ij} = \mu+\alpha_i +\beta_j + \epsilon_{ij}\). Mesmo assim, isso pode ser perigoso se a interação realmente existir. Existem maneiras de testar se a interação é zero neste caso.

Se os níveis de ambos os fatores forem quantitativos, como nos experimentos propostos com helicóptero de papel ou no experimento com etanol combustível, as somas dos quadrados para o termo de interação podem ser particionadas em graus de liberdade polinomiais ortogonais únicos. Por exemplo, se houverem três níveis quantitativos igualmente espaçados do fator A e três níveis quantitativos igualmente espaçados para o fator B, então as somas dos quadrados para a interação podem ser particionadas em quatro graus únicos de liberdade, a saber: linear\(\times\)linear, linear\(\times\)quadrático, quadrático\(\times\)linear e quadrático\(\times\)quadrático.

Usando a filosofia da Série de Taylor de que polinômios de ordem inferior podem aproximar a maioria das relações funcionais, os três termos de ordem superior podem ser assumidos como insignificantes e agrupados para estimar o \(ssE\), que poderia então ser usado como um termo de erro para testar a porção linear\(\times\)linear da interação.

Isto será ilustrado com os dados da experiência do etanol combustível apresentados na Tabela 3.4. Primeiro, considere as médias das duas réplicas em cada célula da Tabela 3.3 como o resultado de um único experimento. O código R mostrado abaixo calcula a média dos dados em cada célula para produzir o conjunto de dados cells com uma observação por célula. Ajustar o modelo 3.2 a estes dados com a função R lm resulta em uma ANOVA com zero graus de liberdade para \(ssE\) e sem testes \(F\) (NaN).

library(daewr)
data(COdata)
Cellmeans = tapply( COdata$CO, list(COdata$Eth, COdata$Ratio), mean )
dim(Cellmeans) = NULL
Eth = factor(rep(c(.1, .2, .3), 3))
Ratio = factor(rep(c(14,15,16), each=3))
cells = data.frame( Eth, Ratio, Cellmeans )
cells
##   Eth Ratio Cellmeans
## 1 0.1    14      64.0
## 2 0.2    14      79.5
## 3 0.3    14      92.0
## 4 0.1    15      69.5
## 5 0.2    15      80.5
## 6 0.3    15      76.5
## 7 0.1    16      67.0
## 8 0.2    16      67.5
## 9 0.3    16      59.0
modnr = lm(Cellmeans ~ Eth*Ratio, data=cells )
anova(modnr)
## Analysis of Variance Table
## 
## Response: Cellmeans
##           Df Sum Sq Mean Sq F value Pr(>F)
## Eth        2    162    81.0     NaN    NaN
## Ratio      2    326   163.0     NaN    NaN
## Eth:Ratio  4    339    84.7     NaN    NaN
## Residuals  0      0     NaN


Para obter as somas dos quadrados para a porção linear\(\times\)linear da interação, os fatores Eth e Ratio são primeiro convertidos em fatores ordenados, conforme mostrado.

Ethc = as.ordered(cells$Eth)
Ratioc = as.ordered(cells$Ratio)


Quando fatores ordenados são usados, a função R lm usa contrastes polinomiais ortogonais para colunas na matriz \(\pmb{X}\) em vez das codificações de tratamento padrão. No código R abaixo, o modelo mbo é ajustado usando apenas os contrastes polinomiais ortogonais lineares\(\times\)lineares para a interação de Ethc e Ratioc.

EthLin = contr.poly(Ethc)[Ethc,".L"]
RatioLin = contr.poly(Ratioc)[Ratioc,".L"]
mbo = lm(Cellmeans~Ethc + Ratioc + EthLin:RatioLin, data = cells)
anova(mbo)
## Analysis of Variance Table
## 
## Response: Cellmeans
##                 Df Sum Sq Mean Sq F value Pr(>F)   
## Ethc             2    162      81    16.2 0.0247 * 
## Ratioc           2    326     163    32.6 0.0092 **
## EthLin:RatioLin  1    324     324    64.8 0.0040 **
## Residuals        3     15       5                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


O erro ou soma residual dos quadrados nesta tabela ANOVA é a diferença nas somas dos quadrados de interação de 4 graus de liberdade mostradas na tabela anterior a esta e o único grau de liberdade linear por somas de quadrados de interação linear.

Essa diferença é usada para construir o denominador dos testes \(F\) da tabela acima. Os resultados mostram que a porção linear por linear da interação é significativa e é responsável pela maior parte das somas dos quadrados da interação.

Como a interação é significativa, o modelo aditivo \(y_{ij} = \mu+\alpha_i +\beta_j +\epsilon_{ij}\) é inadequado e os efeitos dos efeitos principais serão diferentes dependendo do nível do outro fator. Os resultados podem ser melhor interpretados examinando o gráfico de interação.

O gráfico de interação que inclui apenas a parte linear\(\times\)linear da interação pode ser construído traçando as previsões do modelo mbo. No código abaixo, o comando R predict é usado para obter as previsões do modelo e o comando aggregate é usado para criar o conjunto de dados pred.means que combina as previsões do modelo com os níveis dos fatores. A seguir, o comando interact.plot é usado como anteriormente para criar o gráfico.

O resultado é mostrado na Figura 3.7, que deve ser comparado com a Figura 3.6. A Figura 3.7 é bastante semelhante à Figura 3.6, confirmando o que foi visto na tabela ANOVA, ou seja, a maior parte da variação causada pela interação é capturada na parte linear por linear.

Pred = predict(mbo, newdata=data.frame(Ethc, Ratioc, EthLin, RatioLin))
pred.means = aggregate(Pred, by=list(Ethc = Ethc, Ratioc = Ratioc), "mean")
Ethanol = pred.means$Ethc
par(mar = c(4,4,1,1))
interaction.plot(pred.means$Ratioc, Ethanol, pred.means$x, type="b", 
                 pch = c(18,24,22), leg.bty ="o", xlab = "Ratio", ylab = "predicted CO emissions")
grid()

Figura 3.7: Gráfico de interação linear por linear Etanol e relação ar/combustível.

Quando o etanol está em seu nível alto (0.3), o aumento da relação ar/combustível de 14 para 16 causa uma diminuição acentuada nas emissões de CO. Quando o etanol está no seu nível médio (0.2), o aumento da relação ar/combustível de 14 para 16 causa uma ligeira diminuição nas emissões de CO representada pela suave linha inclinada negativa. No entanto, quando o etanol está no seu nível baixo (0.1), o aumento da relação ar/combustível de 14 para 16 causa, na verdade, um aumento nas emissões de CO ilustrado pela linha positivamente inclinada.

Quando há apenas uma réplica por célula em um experimento fatorial e os fatores não têm níveis quantitativos, particionar as somas de interação dos quadrados em contrastes polinomiais ortogonais e combinar os termos de ordem superior como somas de quadrados de erro pode não ser apropriado. No entanto, Tukey (1949b) desenvolveu um método alternativo para testar um único grau de liberdade particionado de somas de quadrados de interação. Este método equivale a restringir o \(\alpha\beta_{ij}\) no modelo 3.2 da Seção 3.5 a ser uma função polinomial de segundo grau dos efeitos principais \(\alpha_i\) e \(\beta_j\), ver Scheffé (1959).

Fazendo isso, as somas dos quadrados \[\begin{equation*} \tag{3.7} ssAB = \dfrac{a\, b \Big( \sum_i\sum_j y_{ij}\overline{y}_{i\bullet}\overline{y}_{\bullet j}-\big(ssA+ssB+a \, b \, \overline{y}_{\bullet\bullet}^2 \big)\overline{y}_{\bullet\bullet}\Big)^2}{(ssA)(ssB)}, \end{equation*}\] para testar a hipótese restrita \(H_0 \, : \, \alpha\beta_{ij} = 0\) para todos \(i\) e \(j\) terá um grau de liberdade e a diferença entre ele e o termo de erro para o modelo aditivo formarão as somas dos quadrados dos erros semelhantes ao exemplo acima com níveis de fator quantitativos.

Para ilustrar o uso do teste de grau de liberdade único de Tukey para interação, considere os dados da Tabela 3.5, que é uma parte dos dados de um estudo para validar um ensaio de contaminação viral relatado por Lin and Stephenson (1998). Ensaios de contaminação viral são usados para determinar a presença e quantidade de um vírus específico em produtos biológicos, como o Fator Oito de coagulação sanguínea.

\[ \begin{array}{cc|cccccc}\hline & & & & \mbox{Amostras} & & & \\ & & 1 & 2 & 3 & 4 & 5 & 6 \\\hline & 3 & 1.87506 & 1.74036 & 1.79934 & 2.02119 & 1.79934 & 1.59106 \\ \mbox{Diluição} & 4 & 1.38021 & 1.36173 & 1.25527 & 1.39794 & 1.20412 & 1.25527 \\ & 5 & 0.60206 & 0.90309 & 0.95424 & 1.00000 & 0.60206 & 0.60206 \\\hline \end{array} \] Tabela 3.5: Ensaio \(\log_{10}\) (PFU/mL) de contaminação viral.

Um experimento ou corrida, consiste em fazer uma solução com uma contaminação viral conhecida, permitir que o vírus cresça em uma solução contaminada e depois medir o resultado. A unidade experimental é a amostra viral específica em combinação com o local e a época onde ela pode crescer. O fator A representa o número da amostra ou solução com a qual a amostra viral é misturada ou enriquecida. O fator B representa diferentes diluições da amostra enriquecida. A resposta medida é o \(\log_{10}\) das unidades formadoras de placa por mL de solução.

Como o fator A (amostras) não é um fator quantitativo, seria inapropriado usar contrastes polinomiais ortogonais para particionar suas somas de quadrados ou as somas de quadrados de sua interação com o fator B (diluição). Para determinar se o modelo aditivo \(y_{ij} = \mu+\alpha_i +\beta_j +\epsilon_{ij}\) é apropriado para esses dados, teste para ver se há uma interação significativa usando o método de Tukey.

A função Tukey1df no pacote R daewr calcula as somas de quadrados não aditivas ou de interação, mostradas na Equação (3.7) e imprime um relatório. O código para abrir os dados da Tabela 3.5 e chamar a função é mostrado abaixo. A primeira coluna no conjunto de dados usado por esta função é uma resposta numérica, a segunda coluna é o indicador do fator A e a terceira coluna é o indicador do fator B.

O número de linhas no quadro de dados deve ser exatamente igual ao número de níveis do fator A vezes o número de níveis do fator B, uma vez que o experimento não possui réplicas.

library(daewr)
data("virus")
virus
##          y Sample Dilution
## 1  1.87506      1        3
## 2  1.38021      1        4
## 3  0.60206      1        5
## 4  1.74036      2        3
## 5  1.36173      2        4
## 6  0.90309      2        5
## 7  1.79934      3        3
## 8  1.25527      3        4
## 9  0.95424      3        5
## 10 2.02119      4        3
## 11 1.39794      4        4
## 12 1.00000      4        5
## 13 1.79934      5        3
## 14 1.20412      5        4
## 15 0.60206      5        5
## 16 1.59106      6        3
## 17 1.25527      6        4
## 18 0.60206      6        5
Tukey1df(virus)
## Source           df     SS        MS        F     Pr>F 
## A                 5   0.1948    0.039 
## B                 2   3.1664    1.5832 
## Error            10   0.1283    0.0513 
## NonAdditivity     1   0.0069    0.0069    0.51    0.4932 
## Residual          9   0.1214    0.0135


Nos resultados percebe-se que a interação ou não aditividade não é significativa. Portanto, para estes dados, seria apropriado ajustar o modelo aditivo, \(y_{ij} =\mu +\alpha_i +\beta_j +\epsilon_{ij}\), com a função R lm ou aov.


3.6 Experimentos fatoriais com múltiplos fatores


Os experimentos fatoriais de dois fatores são mais eficientes do que estudar cada fator separadamente em experimentos de um fator. Da mesma forma, quando muitos fatores estão em estudo, é mais eficiente estudá-los juntos em um planejamento fatorial multifatorial do que estudá-los separadamente em grupos de dois usando planejamentos fatoriais de dois fatores. Quando múltiplos fatores são estudados simultaneamente, o poder de detecção dos efeitos principais aumenta em relação ao que seria em experimentos fatoriais de dois fatores separados.

Além disso, é possível detectar interações entre qualquer um dos fatores. Se os fatores fossem estudados separadamente em fatoriais de dois fatores, as interações de dois fatores só poderiam ser detectadas entre fatores estudados em conjunto no mesmo delineamento. Em um fatorial multifatorial não só é possível detectar interações de dois fatores entre qualquer par de fatores, mas também é possível detectar interações de ordem superior entre grupos de fatores. Uma interação de três fatores entre os fatores A, B e C, por exemplo, significa que o efeito do fator A difere dependendo da combinação dos níveis dos fatores B e C. Exemplos de interações de ordem superior serão apresentados nos exemplos a seguir.

As combinações de tratamento em um fatorial multifatorial consistem em todas as combinações possíveis dos níveis de todos os fatores. Um experimento pode ser produzido usando a função expand.grid em R, semelhante ao plano completamente aleatório criado na Seção 3.4, usando a função gen.factorial no pacote AlgDesign, ou usando funções de outros pacotes que serão descritos posteriormente.

O modelo para análise é uma extensão da Equação (3.2) e a análise pode ser feita utilizando a função R lm semelhante aos exemplos mostrados anteriormente. Considere um exemplo de projeto ou planejamento fatorial multifatorial em pesquisa de marketing. Uma empresa cujas vendas são feitas on-line através de uma página da Web gostaria de aumentar a proporção de visitantes do seu site que se inscrevem para receber seus serviços, configurando sua página da Web de maneira otimizada. Para comprar da empresa, o cliente deve se cadastrar e preencher um formulário informando seu endereço de e-mail e outros campos obrigatórios.

Depois que um cliente se inscreve, a empresa tem informações de contato em seu banco de dados e pode enviar anúncios por e-mail, ofertas especiais e assim por diante. A empresa gostaria de experimentar testando diferentes configurações de sua página da Web para ver se consegue aumentar o número de visitantes de seu site que realmente se inscrevem.

As unidades experimentais neste estudo serão indivíduos que visitam o site da empresa. A resposta é binária; o cliente se inscreve ou não. Os fatores em estudo foram características que alteram a aparência da página web. Por exemplo, o fator A eram as alternativas de plano de fundo para a página com três opções. O fator B foi o tamanho da fonte do banner principal, com três níveis; o fator C foi a cor do texto com duas alternativas; e o fator D foi uma escolha entre um botão ou link de inscrição. Com base nesses fatores havia \(3\times 3\times 2\times 2 = 36\) configurações possíveis da página Web ao considerar todas as combinações possíveis dos níveis de cada fator.

Um experimento fatorial de quatro fatores consistiria em atribuir aleatoriamente os visitantes do site a uma das configurações possíveis e registrar sua resposta binária. Existem variáveis ocultas que podem afetar a chance de um visitante do site se inscrever. Por exemplo, a ordem de posição em que o link, para o site da empresa, aparece em uma pesquisa na Web sobre os produtos que ela vende, as promoções oferecidas pelos concorrentes e a atratividade dos sites dos concorrentes.

A atribuição aleatória de cada visitante sequencial do site a uma das configurações alternativas em estudo deve minimizar a chance de viés de mudanças nas variáveis ocultas ao longo do tempo.

A probabilidade de um visitante do site se inscrever pode ser expressa pelo modelo: \[\begin{equation*} \tag{3.8} \begin{array}{rcl} p_{ijkl} & = & \mu+\alpha_i+\beta_j +\alpha\beta_{ij} +\gamma_k+\alpha\gamma_{ik}+\beta\gamma_{jk}+\alpha\beta\gamma_{ijk}+\delta_l + \alpha\delta_{il} \\ & & \qquad +\beta\delta_{jl}+\alpha\beta\delta_{ijl}+\gamma\delta_{kl}+\alpha\gamma\delta_{ikl}+\beta\gamma\delta_{jkl}+\alpha\beta\gamma\delta_{ijkl}, \end{array} \end{equation*}\] onde \(\alpha_i\) representa o efeito da escolha do plano de fundo, \(\beta_j\) representa o efeito do tamanho da fonte no banner principal, \(\gamma_k\) representa o efeito da cor do texto e \(\delta_l\) representa o efeito do link de inscrição versus botão.

O experimento foi conduzido através da construção de 36 sites consistindo de todas as combinações possíveis dos quatro fatores descritos acima. Cada cliente potencial que visitou o site da empresa durante o período de teste foi redirecionado aleatoriamente para uma das 36 configurações.

O número de visitantes \(n_{ijkl}\) para a configuração \(ijkl\) e o número de inscritos \(x_{ijkl}\) foram registrados. Então \(x_{ijkl}\) têm distribuição binomial \[\begin{equation*} \tag{3.9} \mbox{Binomial}(x_{ijkl};n_{ijkl},p_{ijkl})=\binom{n_{ijkl}}{x_{ijkl}}p_{ijkl}^{x_{ijkl}}(1-p_{ijkl})^{(n_{ijkl}-x_{ijkl})}\cdot \end{equation*}\]

Abaixo está o código R para abrir os dados brutos e imprimir as primeiras seis linhas do conjunto ou quadro de dados.

library(daewr)
data(web)
head(web)
##   A B C D visitors signup
## 1 1 1 1 1     1016     22
## 2 1 1 1 2     1145     16
## 3 1 1 2 1     1145     17
## 4 1 1 2 2     1082     19
## 5 1 2 1 1     1121     28
## 6 1 2 1 2     1103     28


O procedimento correto deve ser utilizado para analisar os dados, determinar se algum dos efeitos dos fatores é significativo e prever a configuração ideal da página da Web. Como as respostas dos visitantes individuais para cada configuração do site são Bernoulli, os dados de resposta agregados são binomiais com tamanhos de amostra grandes e aproximadamente iguais, ou seja, número de visitantes para cada configuração possível da Web.

A transformação arco seno da raiz quadrada mostrada na Tabela 2.4 da Seção 2.6.2 poderia ser aplicada e a função R lm poderia ser usada para análise. No entanto, o problema com a utilização deste procedimento é que as respostas individuais foram somadas para obter as respostas agregadas e ao usar essas respostas binomiais agregadas, não há observações replicadas em nenhuma das células e, portanto, não há maneira de calcular \(ssE\).

Isso seria semelhante a somar ou calcular a média das respostas replicadas em cada célula se os dados fossem distribuídos normalmente, deixando apenas uma observação por célula e nenhuma maneira de calcular \(ssE\).

A alternativa é utilizar o método de máxima verossimilhança para ajustar o modelo (3.8). Isso pode ser feito usando a função R glm. Ele definirá automaticamente \(\sigma^2 = 1.0\) e as somas de quadrados do tipo III da forma \[ \big(\pmb{L}\widehat{\beta} \big)^\top\big(\pmb{L}(\pmb{X}^\top\pmb{X})^{-1} \pmb{L}^\top\big)^{-1} \big(\pmb{L}\widehat{\beta} \big) \] serão distribuídas assintoticamente como qui-quadrados sob a hipótese nula.

Os comandos para analisar os dados usando glm são mostrados a continuação.

modb = glm( cbind(signup, visitors-signup) ~ A * B * C * D, data = web, family = binomial )
summary( modb )
## 
## Call:
## glm(formula = cbind(signup, visitors - signup) ~ A * B * C * 
##     D, family = binomial, data = web)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.8107     0.2155  -17.68   <2e-16 ***
## A2           -0.0927     0.3083   -0.30    0.764    
## A3           -0.2608     0.3260   -0.80    0.424    
## B2            0.1462     0.2883    0.51    0.612    
## B3           -0.3262     0.3208   -1.02    0.309    
## C2           -0.3843     0.3258   -1.18    0.238    
## D2           -0.4458     0.3314   -1.35    0.179    
## A2:B2        -0.2756     0.4304   -0.64    0.522    
## A3:B2        -0.2942     0.4589   -0.64    0.521    
## A2:B3         0.7389     0.4286    1.72    0.085 .  
## A3:B3         0.3492     0.4651    0.75    0.453    
## A2:C2         0.7202     0.4376    1.65    0.100 .  
## A3:C2         0.4806     0.4686    1.03    0.305    
## B2:C2        -0.2064     0.4587   -0.45    0.653    
## B3:C2         0.6378     0.4572    1.40    0.163    
## A2:D2         0.6734     0.4432    1.52    0.129    
## A3:D2         0.5781     0.4671    1.24    0.216    
## B2:D2         0.4624     0.4279    1.08    0.280    
## B3:D2         0.5601     0.4660    1.20    0.229    
## C2:D2         0.6164     0.4724    1.30    0.192    
## A2:B2:C2     -0.1910     0.6396   -0.30    0.765    
## A3:B2:C2      0.6113     0.6564    0.93    0.352    
## A2:B3:C2     -0.7430     0.5948   -1.25    0.212    
## A3:B3:C2     -0.3681     0.6423   -0.57    0.567    
## A2:B2:D2     -0.1329     0.5946   -0.22    0.823    
## A3:B2:D2     -0.1188     0.6339   -0.19    0.851    
## A2:B3:D2     -0.9304     0.6127   -1.52    0.129    
## A3:B3:D2     -0.4225     0.6481   -0.65    0.514    
## A2:C2:D2     -0.8702     0.6200   -1.40    0.160    
## A3:C2:D2     -0.6864     0.6586   -1.04    0.297    
## B2:C2:D2      0.1489     0.6287    0.24    0.813    
## B3:C2:D2     -0.3522     0.6396   -0.55    0.582    
## A2:B2:C2:D2   0.0342     0.8585    0.04    0.968    
## A3:B2:C2:D2  -0.7181     0.8955   -0.80    0.423    
## A2:B3:C2:D2   0.4874     0.8389    0.58    0.561    
## A3:B3:C2:D2   0.3755     0.8833    0.43    0.671    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5.6641e+01  on 35  degrees of freedom
## Residual deviance: 6.2195e-13  on  0  degrees of freedom
## AIC: 251.1
## 
## Number of Fisher Scoring iterations: 3


No arquivo de dados, signup/visitors é a proporção observada de inscrições. A opção family = binomial declara a resposta como distribuída binomialmente e a opção test = “Chisq” na chamada à função anova solicita uma tabela de somas de quadrados tipo III e testes qui-quadrado.

O comando summary(modb) imprime a tabela de estimativas dos parâmetros produzidas pelo método de máxima verossimilhança que, como pode ser verificado na saída acima, é semelhante ao resumo de um objeto criado pela função lm.

anova(update(modb, .~ A+B + A:B + C + A:C + B:C + A:B:C + D + A:D+B:D + 
               A:B:D + C:D + A:C:D + B:C:D + A:B:C:D ), test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cbind(signup, visitors - signup)
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                       35       56.6            
## A        2    10.25        33       46.4   0.0059 **
## B        2     6.93        31       39.5   0.0313 * 
## C        1     2.96        30       36.5   0.0852 . 
## D        1     4.62        29       31.9   0.0317 * 
## A:B      4     6.31        25       25.6   0.1772   
## A:C      2     0.68        23       24.9   0.7110   
## B:C      2     4.03        21       20.9   0.1333   
## A:D      2     0.28        19       20.6   0.8691   
## B:D      2     3.77        17       16.8   0.1518   
## C:D      1     0.08        16       16.7   0.7796   
## A:B:C    4     2.21        12       14.5   0.6967   
## A:B:D    4     6.33         8        8.2   0.1757   
## A:C:D    2     6.15         6        2.0   0.0461 * 
## B:C:D    2     0.02         4        2.0   0.9891   
## A:B:C:D  4     2.01         0        0.0   0.7346   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Nesta saída, podemos ver que, no nível de significância \(\alpha = 0.05\), os fatores A (estilo de fundo) e o fator D (botão ou link de inscrição) foram significativos junto com a interação tripartida A:C:D, onde o fator C representa a cor do texto.

Como existe uma interação significativa, os efeitos principais A e D não podem ser interpretados separadamente. Se é melhor usar um botão ou link de inscrição depende se a cor do texto é preto ou branco e qual estilo de fundo é escolhido.

Para interpretar a interação tripartida, é necessário fazer uma tabela da proporção de inscritos em cada combinação dos fatores A, C e D e uma série de gráficos de interação.

O código R para fazer isso é mostrado abaixo.

prop = web$signup / web$visitors
webp = data.frame(web,prop)
par ( mfrow = c(1,3), mar = c(4,4,2,1) )
webp1 = subset(webp, A == 1)
interaction.plot(webp1$C, webp1$D, webp1$prop, type = "l", legend=FALSE, 
                 ylim = c(.015,.0275), main = "Background = 1", 
                 xlab = "Text Color", ylab = "Proportion Signing-up")
grid()
webp2 = subset(webp, A == 2 )
interaction.plot( webp2$C, webp2$D, webp2$prop, type = "l", legend = FALSE, 
                  ylim = c(.015,.0275), main = "Background = 2", 
                  xlab = "Text Color", ylab = " ")
grid()
lines( c(1.7,1.85), c(.016,.016), lty = 2)
lines( c(1.7,1.85), c(.017,.017) ,lty = 1)
text(1.3, .017, "Sign-up link ")
text(1.3, .016, "Sign-up Button" )
text(1.4, .018, "LEGEND" )
webp3 = subset(webp, A == 3)
interaction.plot(webp3$C, webp3$D, webp3$prop, type = "l", legend=FALSE, 
                 ylim = c(.015,.0275), main="Background = 3", 
                 xlab = "Text Color", ylab = " ")
grid()

Figura 3.8: Efeito de cor do texto por tipo de fundo e botão ou link.

O resultado na figura acima mostra o efeito do fator D (cor do texto) para cada combinação dos níveis dos fatores A (tipo de fundo) e C (link ou botão de inscrição). A forma comum de interpretar a interação é comparar o efeito da variável representada no eixo horizontal entre combinações de níveis dos demais fatores.

Para este exemplo, uma interpretação pode ser feita da seguinte forma. Ao usar o plano de fundo (background) tipo 2 ou tipo 3, pode-se observar que alterar a cor do texto de preto = 1 para branco = 2 causa um aumento na proporção de visitantes do site que se inscrevem. O aumento, representado pela inclinação das linhas, é maior quando um botão de inscrição é usado em vez de um link de inscrição porque a taxa geral de inscrição é maior quando um link é usado, independentemente da cor do texto, e não há muito espaço para melhorias.

Entretanto, quando o tipo de plano de fundo 1 é usado, o efeito da cor do texto é totalmente diferente. Neste caso, alterar a fonte de preto = 1 para branco = 2, na verdade, causa uma diminuição na proporção de inscrições quando um botão de inscrição é usado e há um grande aumento na proporção de inscrições ao mudar da fonte preta para branca ao usar um link de inscrição.

Isto é o oposto dos efeitos observados para os tipos de fundo 2 e 3. Qualquer um dos três factores poderia ser colocado no eixo horizontal e uma interpretação equivalente poderia ser feita. Às vezes, a interpretação resultante da colocação do fator com maior efeito principal no eixo horizontal é mais fácil de explicar.

Ao interpretar uma interação de dois fatores, apenas um gráfico foi necessário para ilustrar o fato de que o efeito de um fator dependia do nível de outro fator. Contudo, neste caso é necessário mais do que um gráfico para ilustrar como o efeito de um factor depende da combinação dos níveis dos outros dois factores. As duas linhas em cada gráfico mostram como o efeito da cor do texto muda quando há um link de inscrição versus um botão, e os diferentes gráficos mostram como o efeito muda quando o plano de fundo é alterado.

A partir da inspeção dos três gráficos ou da tabela de médias que poderiam ser produzidas com a função tapply, pode-se observar que a maior proporção de inscrições seria para a página da Web com link de inscrição, texto em branco e tipo de plano de fundo 2. Aqui prevê-se que, em média, um pouco mais de 2.7% se inscrevam. O tamanho da fonte é insignificante, portanto não importa o tamanho da fonte usado.


3.7 Fatoriais de dois níveis


À medida que fatores adicionais são adicionados a um experimento fatorial, o número de combinações de tratamento ou execuções no experimento aumenta exponencialmente. O exemplo da última seção continha quatro fatores e 36 combinações de tratamento.

Se houvesse cinco fatores em um projeto, cada um com quatro níveis, o número de combinações de tratamento seria \(4\times 4\times 4\times 4\times 4 = 4^5 = 1.024\) execuções no projeto. Pode-se ver que não seriam necessários muitos fatores para tornar o projeto impraticável.

Em outras palavras, haveria muitas combinações de tratamento para serem executadas em um período de tempo razoável. No entanto, é melhor reduzir o número de níveis de cada fator e permanecer com o planejamento fatorial usando todos os fatores do que reverter para experimentos um de cada vez ou dois de cada vez e perder a eficiência dos experimentos fatoriais.

Com experimentos separados, a capacidade de detectar interações de ordem superior e a capacidade de detectar interações entre qualquer par de fatores é perdida. Se cinco fatores em um planejamento fatorial fossem estudados com apenas dois níveis cada, o número de combinações de tratamento seria reduzido para \(2^5 = 32\). Por esse motivo, os planejamentos fatoriais com dois níveis para cada fator ou fatoriais de dois níveis, são populares. Uma abreviação para um fatorial de dois níveis com \(k\) fatores é um planejamento \(2^k\).

Em fatoriais de dois níveis, se um fator tiver níveis quantitativos, os dois níveis são denotados simbolicamente por (-) e (+), onde (-) representa o nível mais baixo que o experimentador consideraria e (+) representa o nível mais alto que o experimentador consideraria. Os níveis altos e baixos são geralmente espalhados tanto quanto possível a fim de acentuar o sinal ou a diferença na resposta entre os dois níveis. Se um fator tem níveis qualitativos, as designações (-) e (+) são arbitrárias, mas os dois níveis escolhidos normalmente seriam aqueles que o experimentador acredita que deveriam resultar na diferença máxima na resposta.


3.7.1 Efeitos principais e inclinações da regressão


O modelo para um experimento fatorial com três fatores pode ser escrito como: \[\begin{equation*} \tag{3.10} y_{ijkl}=\mu+\alpha_i+\beta_j+\alpha\beta_{ij}+\gamma_k+\alpha\gamma_{ik}+\beta\gamma_{jk}+\alpha\beta\gamma_{ijk}+ \epsilon_{ijkl} \end{equation*}\] onde \(\alpha_i\), \(\beta_j\) e assim por diante são os efeitos definidos anteriormente.

No entanto, no caso em que cada fator tem apenas dois níveis representados por (−) e (+), \(i\), \(j\), \(k\) e \(l\) podem ser substituídos por (−) ou (+) e \(\alpha_- = −\alpha_+\), já que \[ \alpha_-=\overline{y}_{-\cdot\cdot\cdot}-\overline{y}_{\cdot\cdot\cdot\cdot}, \quad \alpha_+=\overline{y}_{+\cdot\cdot\cdot}-\overline{y}_{\cdot\cdot\cdot\cdot} \quad \mbox{e} \quad \overline{y}_{\cdot\cdot\cdot\cdot}-=(\overline{y}_{-\cdot\cdot\cdot}+\overline{y}_{+\cdot\cdot\cdot})/2\cdot \]

Uma igualdade semelhante será verdadeira para todos os efeitos e interações. Como os dois efeitos para cada fator são o mesmo valor com sinais diferentes, uma maneira mais compacta de definir os principais efeitos para um fatorial de dois níveis é \(E_A=\overline{y}_{+\cdot\cdot\cdot}-\overline{y}_{-\cdot\cdot\cdot}\). Isto pode ser visualizado no lado esquerdo da Figura 3.9 e representa a mudança na resposta média causada por uma mudança no fator do seu nível baixo (-) para o seu nível alto (+). Este efeito pode então ser representado pela diferença entre duas médias \(\overline{y}_{+\cdot\cdot\cdot}\) e \(\overline{y}_{-\cdot\cdot\cdot}\).

Figura 3.9: Efeito e coeficiente de regressão para o fatorial de dois níveis.

A inclinação da regressão \(\beta_A\) mostrada no lado direito da Figura 3.9 é a mudança vertical na resposta média para uma mudança de uma unidade, ou seja, de 0 a +1 no nível do fator em unidades simbólicas. Portanto, a inclinação \(\beta_A\), é apenas metade do efeito \(E_A\) ou a diferença entre duas médias dividida por 2.

As combinações de tratamento em um fatorial de dois níveis também podem ser representadas geometricamente como os cantos de um cubo, conforme mostrado na Figura 3.10. No lado esquerdo desta figura há uma lista das combinações de tratamento ou execuções listadas na ordem padrão ou de Yates, com a primeira coluna mudando mais rapidamente com sinais alternados - e +, a segunda coluna mudando em pares de sinais - e +, e o terceira coluna mudando mais lentamente em grupos de quatro sinais - e +.

As combinações de tratamento em experimentos fatoriais de dois níveis têm sido tradicionalmente escritas em ordem padrão para facilitar o cálculo manual dos efeitos principais e dos efeitos de interação usando o algoritmo de Yates, ver Daniel (1976). O efeito principal para o fator A pode ser visualizado na figura como a diferença entre a média das respostas do lado direito do cubo nos círculos sombreados em cinza e a média das respostas do lado esquerdo do cubo nos círculos brancos.

Com programas de computador modernos, como a função R lm, metade dos efeitos principais ou coeficientes de regressão, mostrados no lado direito da Figura 3.9, podem ser calculados por regressão e não precisamos mais do algoritmo de Yates.

Figura 3.10: Representação geométrica de \(2^3\) projetos e cálculo do efeito principal.

Uma das propriedades desejáveis de um plano fatorial \(2^k\) é que os efeitos dos fatores não sejam obscurecidos por mudanças planejadas em outros fatores. Na lista de experimentos para planejamento \(2^k\), mostrada na Figura 3.10, isso é evidente pelo fato de que no nível alto de cada fator, há um número igual de níveis altos e baixos de todos os outros fatores.

Também no nível baixo de cada fator, há um número igual de níveis altos e baixos de todos os outros fatores. Assim, o efeito de um fator ou a diferença na resposta média entre o nível alto e baixo desse fator, representa o efeito daquele fator sozinho, porque a influência de todos os outros fatores foi calculada em média. Matematicamente esta propriedade é chamada de ortogonalidade.


3.7.2 Interações


Quando todos os fatores têm apenas dois níveis, o efeito da interação AB é definido como metade da diferença no efeito simples do fator A, \((\overline{y}_{++\cdot\cdot}-\overline{y}_{-+\cdot\cdot})\) quando o fator B é mantido constante em seu nível alto (+), e o efeito simples do fator A \((\overline{y}_{+-\cdot\cdot}-\overline{y}_{--\cdot\cdot})\), quando o fator B é mantido constante em seu nível baixo (-), ou seja, \[ \dfrac{1}{2}\big((\overline{y}_{++\cdot\cdot}-\overline{y}_{-+\cdot\cdot})-(\overline{y}_{+-\cdot\cdot}-\overline{y}_{--\cdot\cdot})\big) \cdot \] Isso é ilustrado no lado esquerdo da Figura 3.11.

O efeito de interação também pode ser definido como metade da diferença no efeito simples do fator B, calculado como \((\overline{y}_{++\cdot\cdot}-\overline{y}_{+-\cdot\cdot})\), quando o fator A é mantido constante em seu nível alto (+) e o efeito simples do fator B, \((\overline{y}_{-+\cdot\cdot}-\overline{y}_{--\cdot\cdot})\), quando o fator A é mantido constante em seu nível baixo (-).

Isso é ilustrado no lado direito da Figura 3.11. De qualquer forma, o efeito de interação é \[ E_{AB} = \dfrac{1}{2}\big(\overline{y}_{++\cdot\cdot}+\overline{y}_{--\cdot\cdot}\big)-\dfrac{1}{2}\big(\overline{y}_{+-\cdot\cdot}+\overline{y}_{-+\cdot\cdot}\big), \] que é a diferença de duas médias.

Figura 3.11: Definição de um efeito de interação para fatorial de dois níveis.

Pode-se determinar quais respostas devem ser calculadas e quais médias devem ser subtraídas das outras para calcular um efeito de interação conforme ilustrado na Figura 3.12. Para calcular a interação AB, adicionamos uma coluna de sinais \(X_A\times X_B\), à lista de combinações de tratamento no lado esquerdo da figura.

Figura 3.12: Representação geométrica de planejamentos \(2^3\) e efeito de interação.

Os elementos nesta nova coluna são apenas os produtos elementares dos sinais na coluna para \(X_A\) e \(X_B\), ou seja, (−)(−) = +, (−)(+) = − etc. Agora o efeito de interação pode ser visualizado na figura como a diferença na resposta média em uma diagonal representada por círculos cinza e a resposta média na outra diagonal representada por círculos brancos.

A partir desta representação, também pode ser visto que os efeitos da interacção não são obscurecidos por mudanças planeadas noutros fatores ou, por outras palavras, são ortogonais aos efeitos principais.

Metade deste efeito de interação ou o coeficiente de regressão pode ser calculado usando um programa de regressão como a função R lm adicionando um termo \(X_A\times X_B\) ao modelo. Os efeitos de interação de ordem superior podem ser definidos de forma semelhante. Portanto, uma maneira mais simples de escrever o modelo para um fatorial de dois níveis é usar a equação de regressão familiar, \[\begin{equation*} \tag{3.11} y = \beta_0+\beta_A X_A + \beta_B X_B +\beta_{AB}X_AX_b +\beta_C X_C + \beta_{AC}X_{A}X_C \\ \qquad \qquad \qquad \qquad + \beta_{BC}X_B X_C +\beta_{ABC}X_A X_B X_C+\epsilon, \end{equation*}\] onde os \(\beta\)’s são metade dos efeitos e \(X_A = −1\) se o fator A estiver em seu nível baixo e \(X_A = +1\) se o fator A estiver em seu nível alto.

Se escrevermos este modelo em termos matriciais, \(y = X\beta+\epsilon\), a propriedade de ortogonalidade do projeto é expressa pelo fato de que as colunas da matriz \(X\) são ortogonais e a matriz \(X^\top X\) é diagonal com elementos diagonais \(r2^k\), onde \(r\) é o número de réplicas de cada célula.


3.7.3 Exemplo de um fatorial \(2^3\)


Para ilustrar o projeto e a análise de um experimento fatorial \(2^3\), considere o seguinte exemplo em Lawson and Erjavec (2001). Estudantes de um laboratório de eletrônica universitária frequentemente reclamavam que as medições de tensão feitas em um circuito que construíram em sala de aula eram inconsistentes.

O auxiliar de ensino do laboratório (TA - teaching assistant) decidiu realizar um experimento para tentar identificar a origem da variação. Os três fatores que ele variou foram A = temperatura ambiente onde a medição de tensão foi feita, B = tempo de aquecimento do voltímetro e C = tempo em que a energia foi conectada ao circuito antes da medição ser feita.

A resposta foi a tensão medida em milivolts. Os dois níveis para o fator A foram − = 22℃ (temperatura ambiente) e + = 32℃ (próximo à temperatura em alguns ambientes industriais). Foi utilizado um forno e o circuito foi deixado estabilizar por pelo menos cinco minutos antes das medições. As configurações para os fatores B e C foram − =30 segundos ou menos e + =5 minutos.

O mesmo circuito foi medido para cada combinação de fatores de tratamento, de modo que a unidade experimental nada mais era do que o ensaio ou momento em que a combinação específica de níveis de fator de tratamento foi aplicada para fazer a medição.

Duas réplicas de cada uma das oito combinações experimentais foram executadas em ordem aleatória para ajudar a evitar vieses. Os resultados do experimento são mostrados na Tabela 3.6.

\[ \begin{array}{c}\hline \begin{array}{cccccc} & \mbox{Nǘveis dos} & \mbox{Códigos dos} & & & & & & \\ & \mbox{fatores} & \mbox{fatores} & & & & & & \end{array} \\ \begin{array}{cccccccccc} \mbox{Corrida} & \mbox{A} & \mbox{B} & \mbox{C} & X_A & X_B & X_C & \mbox{Rep} & \mbox{Ordem} & y \\ \hline 1 & 22 & 0.5 & 0.5 & - & - & - & 1 & 5 & 705 \\ 2 & 32 & 0.5 & 0.5 & + & - & - & 1 & 14 & 620 \\ 3 & 22 & 5.0 & 0.5 & - & + & - & 1 & 15 & 700\\ 4 & 32 & 5.0 & 0.5 & + & + & - & 1 & 1 & 629 \\ 5 & 22 & 0.5 & 5.0 & - & - & + & 1 & 8 & 672 \\ 6 & 32 & 0.5 & 5.0 & + & - & + & 1 & 12 & 668 \\ 7 & 22 & 5.0 & 5.0 & - & + & + & 1 & 10 & 715 \\ 8 & 32 & 5.0 & 5.0 & + & + & + & 1 & 9 & 647 \\ 1 & 22 & 0.5 & 0.5 & - & - & - & 1 & 4 & 680 \\ 2 & 32 & 0.5 & 0.5 & + & - & - & 1 & 7 & 651 \\ 3 & 22 & 5.0 & 0.5 & - & + & - & 1 & 2 & 685 \\ 4 & 32 & 5.0 & 0.5 & + & + & - & 1 & 3 & 635 \\ 5 & 22 & 0.5 & 5.0 & - & - & + & 1 & 11 & 654 \\ 6 & 32 & 0.5 & 5.0 & + & - & + & 1 & 16 & 691 \\ 7 & 22 & 5.0 & 5.0 & - & + & + & 1 & 6 & 672 \\ 8 & 32 & 5.0 & 5.0 & + & + & + & 1 & 13 & 673 \\ \end{array} \end{array} \] Tabela 3.6: Configurações de fator e resposta para experimento com voltímetro.

Nesta tabela, as configurações reais dos fatores são mostradas à esquerda e os níveis codificados - e + são mostrados à direita. As configurações reais à esquerda formam uma lista de receitas ou instruções para realizar cada experimento.

O número de ordem na extrema direita ao lado da resposta foi criado com um gerador de números aleatórios e representa a ordem em que os experimentos devem ser executados. Os níveis de fator codificados são usados como variáveis independentes em um programa de regressão para calcular os coeficientes de regressão ou meios efeitos.

Os níveis de fator codificados podem ser calculados a partir das configurações de fator reais usando a fórmula de codificação e escala. Nesta fórmula, subtraímos o ponto médio das duas configurações de fator e depois dividimos pela metade do intervalo.

Por exemplo, para o fator A, o ponto médio entre 22 e 32 é 27 e metade do intervalo é 5, portanto \[ X_A = \left( \dfrac{\mbox{Configuração do fator real}-27}{5} \right)\cdot \]

A função R contr.FrF2 executa essa codificação e escalonamento em fatores R. Um conjunto de dados volt, no pacote daewr, contém fatores R com os níveis reais dos fatores e resposta da Tabela 3.6. O código para abrir o conjunto de dados, codificar e dimensionar os fatores e ajustar o modelo de regressão com a função lm, junto com a tabela resultante de coeficientes de regressão, é mostrado abaixo.

A função contr.FrF2 rotula os contrastes lineares codificados e escalonados A1, B1 e C1 na saída, em vez de \(X_A\), \(X_B\) e \(X_C\) como na equação 3.11.

library(daewr)
modv = lm( y ~ A*B*C, data = volt, contrast = list(A=contr.FrF2, B=contr.FrF2, C=contr.FrF2))
summary(modv)
## 
## Call:
## lm.default(formula = y ~ A * B * C, data = volt, contrasts = list(A = contr.FrF2, 
##     B = contr.FrF2, C = contr.FrF2))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -21.5  -11.8    0.0   11.8   21.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  668.562      4.518  147.99  4.9e-15 ***
## A1           -16.813      4.518   -3.72   0.0059 ** 
## B1             0.937      4.518    0.21   0.8408    
## C1             5.437      4.518    1.20   0.2632    
## A1:B1         -6.687      4.518   -1.48   0.1771    
## A1:C1         12.563      4.518    2.78   0.0239 *  
## B1:C1          1.813      4.518    0.40   0.6988    
## A1:B1:C1      -5.813      4.518   -1.29   0.2342    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.1 on 8 degrees of freedom
## Multiple R-squared:  0.772,  Adjusted R-squared:  0.572 
## F-statistic: 3.87 on 7 and 8 DF,  p-value: 0.0385


Na saída, pode-se observar que o fator A (temperatura ambiente) e a interação A×C, ou interação entre a temperatura ambiente e o tempo de aquecimento do circuito, são significativas. O efeito principal tem interpretação direta.

O efeito do fator A é duas vezes o coeficiente de regressão mostrado acima ou \[ E_A = 2\times \widehat{\beta}_A = 2(-16.8125)=-33.625\cdot \] Isso significa que, em média, quando a temperatura ambiente aumenta de 22℃ para 32℃, a medição de tensão diminuirá em 33.6 milivolts. Contudo, como a interação é significativa neste exemplo, não faz sentido falar sobre o efeito principal médio porque o efeito da temperatura ambiente depende do tempo de aquecimento do circuito.

Descrever ou interpretar a interação é melhor feito observando o gráfico de interação mostrado na Figura 3.13. Aqui pode ser visto que quando o tempo de aquecimento do circuito é curto, 0.5 minutos ou 30 segundos, a mudança da temperatura ambiente de 22℃ para 32℃ causa uma grande diminuição, 58.7 milivolts na leitura da tensão.

Entretanto, quando o tempo de aquecimento do circuito é longo, 5 minutos, a mudança da temperatura ambiente de 22℃ para 32℃ causa apenas uma pequena diminuição, 8.5 milivolts, na leitura da tensão. Portanto, para tornar as leituras de tensão mais consistentes, o TA do laboratório recomendou que seus alunos deixassem seus circuitos aquecerem 5 minutos antes de fazer medições de tensão.

A regressão foi realizada nos níveis dos fatores codificados para que os coeficientes de regressão produzidos pela função lm sejam exatamente metade dos efeitos. No entanto, os nomes e níveis reais dos fatores devem ser usados para maior clareza ao apresentar os resultados graficamente para inclusão num relatório ou apresentação, conforme mostrado na Figura 3.13. A maioria dos leitores ou ouvintes não se lembrará do que os níveis - e + representam.

O código para produzir a Figura 3.13 é mostrado a continuação. Neste código não houve necessidade de produzir uma tabela de valores previstos a partir do modelo reduzido que contenha apenas os fatores significativos A e AC, como foi feito na criação da Figura 3.7, pois o desenho ortogonal garante que a resposta média nos quatro A por As combinações C serão iguais aos valores previstos do modelo reduzido. Um gráfico de interação ainda mais simples pode ser feito rapidamente com o comando (IAPlot(modv, select = c(1,3)) usando a função IAPlot no pacote FrF2.

C_Warmup = volt$C
par(mar=c(4,4,1,1))
with(volt, (interaction.plot(A, C_Warmup, y, type = "b", pch = c(24,22), leg.bty = "o",
                             xlab = "Temperature",ylab = "Voltage")))
## NULL
grid()


A ortogonalidade do projeto também permite que uma equação de predição reduzida seja escrita a partir dos resultados da regressão, simplesmente eliminando os termos insignificantes.

Esta equação pode ser usada para prever a leitura de tensão em milivolts para qualquer temperatura ambiente entre 22℃ e 32℃ e qualquer tempo de aquecimento do circuito entre 30 segundos e 5 minutos.

\[ y = 668.563-16.813\left(\dfrac{\mbox{Temp}-27}{5}\right)-6.688\left(\dfrac{\mbox{CWarm}-2.75}{2.25}\right)\left(\dfrac{\mbox{Temp}-27}{5}\right)\cdot \]


3.7.4 Fórmula para determinação do número de réplicas


Wheeler (1974) desenvolveu um atalho de aproximação à fórmula para calcular o número de execuções necessárias para atingir um poder igual a 0.95 quando o nível de significância para um fatorial de dois níveis é \(\alpha= 0.05\). Esta fórmula é \[\begin{equation*} \tag{3.12} N=\left( \dfrac{8\sigma}{\Delta}\right)^2, \end{equation*}\] onde \(\sigma\) é o desvio padrão do erro experimental, \(\Delta\) é o tamanho prático de um efeito. Neste caso, a diferença na resposta média entre o nível baixo e alto de um fator e \(N = r\times 2^k\) é o número total de experimentos no fatorial \(2^k\).

Como esta fórmula é tão compacta, é fácil de usar no local em reuniões onde experimentos estão sendo planejados. Como ilustração de seu uso, considere um exemplo baseado no experimento do medidor de tensão apresentado na última seção.

O instrutor do laboratório sentiu que o desvio padrão do erro experimental era cerca de \(\sigma= 15.0\) e o tamanho prático de um efeito era cerca de \(\Delta= 30.0\). \(\sigma\) seria conhecido pela experiência dos instrutores de laboratório em fazer medições repetidas de tensão do mesmo circuito exatamente sob as mesmas condições, ou seja, níveis de fator e \(\Delta\) seria conhecido pela quantidade de inconsistência nas medições reivindicadas pelos alunos que estavam ficando com leituras inconsistentes. \(\Delta\) é o tamanho do efeito que o TA gostaria de detectar em seus experimentos.

Usando a fórmula de atalho, isso diz que \[ N = ((8\times 15.0)/30.0)^2 = 16 \] ou que \(r = 2\) réplicas de cada uma das \(2^3 = 8\) execuções devem resultar em um poder de 0.95 para detectar efeitos de tamanho 30.0 em nível de significância \(\alpha= 0.05\).

Usando uma fórmula mais exata, como mostrada na Seção 3.5.2, o poder real para \(r = 2\) réplicas está mais próximo de 0.94 do que de 0.95. No entanto, esta fórmula aproximada é suficientemente precisa para a maioria dos propósitos de planeamento.

A fórmula simples também pode ser usada de trás para frente resolvendo \(\Delta\) como uma função de \(N\), ou seja, \(\Delta = 8\times \sigma /sqrt{N}\). Dessa forma, se um experimentador conhece seu orçamento para experimentação, que determina o maior \(N\) que pode ser, ele pode calcular o tamanho do efeito \(\Delta\) que provavelmente será capaz de detectar. Se o experimentador não tiver uma estimativa precisa de \(\sigma\), a fórmula ainda pode ser usada falando sobre o tamanho do efeito prático em unidades do desconhecido \(\sigma\).

Por exemplo, se o orçamento de um experimentador lhe permite fazer no máximo \(N = 64\) experimentos, ele pode esperar detectar efeitos que não sejam mais do que um desvio padrão do erro experimental, ou seja, \(\Delta= 8\times \sigma/\sqrt{64} = \sigma\).

Este resultado será verdadeiro independentemente do número de fatores no experimento de dois níveis. Consequentemente, com 64 execuções ele pode ter um fator com \(r = 32\) réplicas de cada nível, ou seis fatores com \(r = 1\) réplica de cada uma das \(2^6 = 64\) combinações de tratamento.


3.7.5 Análise com uma réplica por célula


Projetos fatoriais com uma réplica por célula são frequentemente chamados de experimentos não replicados. Quando há poder adequado para detectar efeitos com \(r = 1\) replicação por célula ou combinação de tratamentos, não há necessidade de duplicar o trabalho experimental replicando cada experimento.

Contudo, num factorial não replicado, surge o mesmo problema que foi discutido na Secção 3.5.4. Haverá zero graus de liberdade para calcular o \(ssE\) e, portanto, nenhum teste \(F\) para os efeitos. Porém, quando existem múltiplos fatores em um fatorial de dois níveis, existem ferramentas gráficas simples que permitem a detecção dos efeitos significativos.

Como não se espera que todos os efeitos e interações principais em um experimento \(2^k\) sejam significativos, os níveis de fatores não significantes e as combinações de níveis definidos pelas interações não significantes são equivalentes a ter réplicas no projeto. As ferramentas gráficas permitem reconhecer os efeitos significativos ou coeficientes de regressão equivalentes.

A ferramenta gráfica mais comum usada para detectar efeitos significativos são os gráficos normais ou semi-normais que foram sugeridos pela primeira vez por Daniel (1959). Estes são fáceis de produzir usando a função DanielPlot no pacote R FrF2 (Groemping, 2011a) ou a função LGB no pacote daewr. Ferramentas gráficas adicionais, como Lenth Plots e Bayes Plots, também são úteis para detectar efeitos e interações significativas e podem ser geradas usando funções do pacote BsMD (Barrios, 2009).

Para ilustrar a análise de um fatorial de dois níveis não replicado, considere um exemplo da indústria química. Os princípios de projeto experimental foram desenvolvidos por R.A. Fisher no início do século XX e foram originalmente usados em experimentos agrícolas.

Dentro de 40 anos houve amplo uso de técnicas de planejamento experimental na indústria química. A Figura 3.14 é um diagrama de um processo químico contínuo. Neste processo, fluxos contínuos de dois reagentes, A e B, são combinados em uma junção chamada mistura-T, onde eles começam a reagir.

Figura 3.14: Diagrama de um processo químico.

A mistura então flui para um reator e é combinada com solvente e um catalisador e a reação é completada. O resultado da reação flui para um tanque separador onde o produto final flutua até o topo em uma fase solvente enquanto o catalisador e a água vão para o fundo do tanque. O catalisador é concentrado e enviado de volta ao reator, enquanto o produto, subprodutos e solvente são levados para uma coluna de destilação onde o produto é removido e o solvente é reciclado para o reator.

Um dos problemas vivenciados nesse processo foi a produção de subproduto (alcatrão). Com o tempo, esses alcatrões obstruiriam o reator e forçariam o encerramento do processo de limpeza. Também exigiu uma etapa adicional do processo para purificar o produto final.

Os engenheiros decidiram realizar experimentos para ver se conseguiriam aumentar a conversão percentual, o que reduziria o valor de subprodutos. Os fatores que eles consideraram que poderiam afetar a conversão percentual são mostrados na tabela abaixo.

\[ \begin{array}{cl} \mbox{Símbolo} & \mbox{Nome do fator} \\\hline \mbox{A} &\mbox{Excesso de Reagente A (acima da quantidade molar)} \\ \mbox{B} & \mbox{Concentração catalisadora} \\ \mbox{C} & \mbox{Pressão no rReator} \\ \mbox{D} & \mbox{Temperatura da Mistura-T Revestida} \end{array} \]

Foram escolhidos dois níveis de cada fator que foram distribuídos tão amplamente quanto os engenheiros consideraram viável, a fim de maximizar a chance de detectar efeitos de fator com apenas dois níveis. Durante a experimentação, os níveis dos fatores seriam alterados após um intervalo fixo de tempo.

A unidade experimental para isso seriam os reagentes, o catalisador e o solvente específico que entram na zona de reação durante uma determinada corrida e a resposta \(Y\), seria a porcentagem de conversão calculada a partir do produto produzido durante uma corrida.

Considerou-se que se a percentagem de conversão pudesse ser aumentada em \(\Delta = 12%\) ou mais, reduziria substancialmente a manutenção e o processamento extra actualmente necessários e valeria a pena detectar. A partir de experiências anteriores com o processo, o desvio padrão na conversão percentual neste processo para produtos produzidos nos mesmos intervalos de comprimento que as execuções no projeto experimental proposto, sem alterações nos níveis dos fatores, foi \(\sigma = 6%\).

Usando a fórmula de atalho, o número de execuções necessárias para ter um poder de 0.95 para detectar efeitos de fator de \(\Delta = 12%\) ou mais foi \[ N= (8\sigma/\Delta)^2 = \big((8)(6)/12\big)^2 = 16\cdot \]

Dezesseis execuções com quatro fatores significam que o projeto não seria replicado. As condições experimentais, em unidades codificadas, uma lista de ordens de execução aleatórias e os resultados dos experimentos são mostrados na Tabela 3.7.

\[ \begin{array}{cccccc}\hline \mbox{No. Corrida aleatória} & \mbox{A} & \mbox{B} & \mbox{C} & \mbox{D} & Y \\\hline 15 & - & - & - & - & 45 \\ 13 & + & - & - & - & 41 \\ 11 & - & + & - & - & 90 \\ 1 & + & + & - & - & 67 \\ 10 & - & - & + & - & 50 \\ 2 & + & - & + & - & 39 \\ 3 & - & + & + & - & 95 \\ 12 & + & + & + & - & 66 \\ 16 & - & - & - & + & 47 \\ 8 & + & - & - & + & 43 \\ 9 & - & + & - & + & 95 \\ 14 & + & + & - & + & 69 \\ 6 & - & - & + & + & 40 \\ 5 & + & - & + & + & 51 \\ 7 & - & + & + & + & 87 \\ 4 & - & + & + & + & 72 \\ \end{array} \] Tabela 3.7: Lista de experimentos e resultados para processo químico.

O código R para abrir o conjunto de dados e realizar uma análise de regressão para estimar os efeitos meios ou coeficientes de regressão é mostrado abaixo. O conjunto de dados chem no pacote daewr contém vetores numéricos de valores codificados e escalonados para os níveis dos fatores e a função contr.FrF2 não era necessária.

library(daewr)
modf = lm( y ~ A*B*C*D, data = chem)
summary(modf)
## 
## Call:
## lm.default(formula = y ~ A * B * C * D, data = chem)
## 
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  62.3125        NaN     NaN      NaN
## A            -6.3125        NaN     NaN      NaN
## B            17.8125        NaN     NaN      NaN
## C             0.1875        NaN     NaN      NaN
## D             0.6875        NaN     NaN      NaN
## A:B          -5.3125        NaN     NaN      NaN
## A:C           0.8125        NaN     NaN      NaN
## B:C          -0.3125        NaN     NaN      NaN
## A:D           2.0625        NaN     NaN      NaN
## B:D          -0.0625        NaN     NaN      NaN
## C:D          -0.6875        NaN     NaN      NaN
## A:B:C        -0.1875        NaN     NaN      NaN
## A:B:D        -0.6875        NaN     NaN      NaN
## A:C:D         2.4375        NaN     NaN      NaN
## B:C:D        -0.4375        NaN     NaN      NaN
## A:B:C:D      -0.3125        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:   NaN 
## F-statistic:  NaN on 15 and 0 DF,  p-value: NA


Nos resultados acima, ao contrário do exemplo da Seção 3.7.3, não haverá estimativa de \(ssE\) e, portanto, não haverá testes \(t\) nos coeficientes de regressão no resumo do lm.

Os coeficientes de regressão para os efeitos principais A e B, juntamente com a interação A×B, são os maiores efeitos, mas um gráfico deve ser usado para determinar quais são significativos.

Os efeitos em um fatorial de dois níveis são a diferença de duas médias. Se as mudanças nos níveis dos fatores não causarem uma mudança na resposta, o efeito será apenas a diferença nas médias dos dados aleatórios, devido a flutuações aleatórias no erro experimental.

Se nenhum dos fatores ou interações causar alterações na resposta, todo o conjunto de efeitos ou coeficientes de regressão, deverá aparecer como uma amostra da distribuição normal com média zero devido ao Teorema do Limite Central. Portanto, se fizermos um gráfico de probabilidade normal dos efeitos, os efeitos não significativos deverão estar ao longo de uma linha reta e quaisquer efeitos ou interações significativas deverão aparecer como valores discrepantes no gráfico.

O código R a seguir ilustra como criar um gráfico normal dos coeficientes de regressão no objeto modf que foi criado pela função lm. Este código chama a função fullnormal do pacote daewr. Um gráfico semelhante, com o eixo invertido, pode ser feito com a função DanielPlot no pacote FrF2.

Neste gráfico normal, Figura 3.15, a maioria dos pontos está ao longo de uma linha reta traçada através da origem em (0, 0). No entanto, os pontos que representam os efeitos principais A, B e a interação A × B tendem a estar abaixo e à esquerda ou acima e à direita da linha reta. Isso indica que esses três efeitos são significativos. Os pontos ao longo da reta são não significativos e a inclinação desta reta é uma estimativa do erro padrão de um efeito \(\widehat{\sigma}_\beta\).

library(daewr)
par(mar=c(4,4,1,1), pch =19)
fullnormal(coef(modf)[-1],alpha=.025)
grid()

Figura 3.15: Gráfico de probabilidade normal de coeficientes de regressão.

Outra maneira de observar os efeitos graficamente é criar um gráfico meio normal ou semi-normal. O gráfico semi-normal é criado essencialmente dobrando o gráfico normal em 45&deg ao longo da linha \(y = -x\). Isso coloca a origem do gráfico no canto inferior esquerdo e faz com que a linha reta dos efeitos não significativos pareça sair da origem.

Por esta razão é fácil adicionar uma linha de referência ao gráfico e os efeitos significativos são identificados como aqueles acima da linha de referência na extrema direita. Para fazer um gráfico semi-normal, represente graficamente o valor absoluto dos efeitos ou coeficientes de regressão no eixo vertical e suas pontuações semi-normais correspondentes no eixo horizontal.

Uma maneira de criar um gráfico meio normal dos coeficientes no objeto R modf é usar a função LGB no pacote daewr. Por padrão, esta função desenha um limite de significância e produz um relatório indicando quais efeitos são significativos. A opção rpt = FALSE suprime o relatório impresso. O código abaixo produz a Figura 3.16.

library(daewr)
par(mar=c(4,4,1,1), pch =19)
LGB( coef(modf)[-1], rpt = FALSE)
grid()

Figura 3.16: Gráfico meio Normal do módulo absoluto dos coeficientes de regressão.

O gráfico semi-normal mostra que os efeitos A, B e a interação A × B caem acima da linha limite superior de predição e são claramente significativos. A linha de referência e o limite superior da Figura 3.16 foram adicionados automaticamente pela função LGB utilizando o método descrito por Lawson et al. (1998).

Ao desenhar a linha de referência manualmente em um gráfico é mais fácil usar o gráfico semi-normal do que o gráfico normal, porque a primeira metade a dois terços dos pontos com tendência a partir do canto inferior esquerdo quase sempre formarão uma linha reta. Porém, no gráfico semi-normal os sinais dos coeficientes são perdidos.

Por exemplo, no resumo do lm e no gráfico normal, pode-se ver claramente que o efeito principal A, o excesso do reagente A, tem um efeito negativo e que aumentar o excesso do reagente A causará, em média, uma diminuição na percentagem de conversão. No gráfico semi-normal, pode-se observar que o efeito principal A é significativo, mas é necessário consultar a tabela de coeficientes para ver se tem um efeito positivo ou negativo.

É claro que neste exemplo há uma interação significativa e, portanto, os efeitos principais não podem ser interpretados separadamente. Para auxiliar na interpretação da interação, um gráfico de interação deve ser criado.

O código R para criar o gráfico de interação entre o fator A, excesso do reagente A, e o fator B, concentração do catalisador, está abaixo e o gráfico resultante é mostrado na Figura 3.17.

library(daewr)
par(mar=c(4,4,1,1))
with(chem, (interaction.plot( A, B, y, type = "b", pch = c(18,24), 
                              main = "Interaction Plot of Catalyst by Excess A",
                              xlab = "Excess Reactant A", ylab = "Percent Conversion")))
## NULL
grid()

Figura 3.17: Gráfico de interação catalisador por excesso A.

Esta figura mostra que aumentar o nível de catalisador aumenta a conversão. Aumentar o excesso do reagente A tem pouco efeito na conversão quando é utilizado um nível baixo de catalisador. Contudo, se for utilizado um nível elevado de catalisador, o aumento do excesso do reagente A diminui a conversão em mais de 20%.

Portanto, para atingir o nível mais alto de conversão, um nível elevado de catalisador e um nível baixo de reagente A em excesso devem ser utilizados. Duas outras ferramentas gráficas usadas para identificar efeitos e interações significativas em fatoriais não replicados são o gráfico de Lenth (Lenth, 1989) e o gráfico de Bayes (Box and Meyer, 1986a).

O gráfico Lenth é semelhante ao gráfico de análise de médias (Ott, 1967) com limites adicionais fornecidos. É pouco provável que as estimativas que se enquadram nos limites da margem de erro (ME) sejam significativas; as estimativas que estão dentro dos limites da margem de erro simultânea (SME), mas fora dos limites da margem de erro (ME), são possivelmente significativas, enquanto aquelas que estão fora dos limites da margem de erro são provavelmente significativas.

O gráfico de Bayes representa graficamente as probabilidades a posteriori de Bayes de que os efeitos estejam ativos. As funções LenthPlot e BsProb do pacote BsMD facilitam a confecção desses gráficos. O código R para fazer esses gráficos com os dados do experimento de processo químico é mostrado abaixo.

par( mfrow = c(2,1), mar = c(1,1,1,1) )
library(BsMD)
par( mar = c(4,4,1,1))
LenthPlot(modf, main = "Lenth Plot of Effects")
##  alpha    PSE     ME    SME 
## 0.0500 1.6875 4.3379 8.8065
grid()
X = model.matrix(modf)[ , 2:16]
y = chem$y
Chem.BsProb = BsProb( X = X, y = y, blk = 0, mFac = 15, 
                      mInt = 1, p = 0.2, g = 2.49, ng = 1, nMod = 10)
par(mar = c(4,4,1,1))
plot( Chem.BsProb, main = "Bayes Plot of Effects" )
grid()

Figura 3.18: Gráficos de efeitos de Lenth e Bayes do experimento de processo químico.

O resultado é mostrado na Figura 3.18 acima, onde pode ser visto que os efeitos A, B e a interação AxB são identificados como provavelmente significativos no gráfico de Lenth. O gráfico de Bayes rotula o efeito pela ordem de Yates, ou seja, x1=A, x2=B, x3=C, x4=D, x5=A:B, x6=A:C, … etc. em vez de pelo nome.

No entanto, os mesmos efeitos A, B e AxB são identificados como tendo grandes probabilidades posteriores. Esses resultados são consistentes com o que foi observado no gráfico de efeitos normal e meio normal.


3.8 Verificando as suposições do modelo


Quando houver unidades experimentais replicadas em cada célula de um modelo fatorial ou quando um termo de interação puder ser considerado desprezível e removido do modelo, como resultado de um teste preliminar como aqueles descritos na Seção 3.5.4 ou Seção 3.7.5, as suposições de normalidade e variância constante do modelo fatorial podem ser verificadas com gráficos de resíduos conforme descrito na Seção 2.4.

Porém, no caso do experimento \(2^k\) com apenas uma réplica por célula é um pouco mais difícil verificar a suposição de normalidade. A suposição de normalidade é mais frequentemente violada por ter um valor discrepante ou atípico. Os principais efeitos e interações calculados em um fatorial de dois níveis sempre podem ser representados como a diferença de duas médias, \(\overline{y}_+ =\overline{y}_−\).

Quando os erros experimentais seguem uma distribuição normal, os efeitos calculados para fatores e interações que têm uma influência insignificativa na resposta devem ser normalmente distribuídos com média zero. A significância dos potenciais factores influentes é julgada pela sua relação com uma linha de pontos de referência num gráfico normal ou semi-normal de efeitos formado pelos factores negligenciáveis. No entanto, se acontecer um valor atípico este irá distorcer cada efeito calculado positiva ou negativamente longe de zero.

A variabilidade dos efeitos calculados para os fatores e interações não influentes será muito maior e será mais difícil julgar a significância dos efeitos, muito menos verificar a suposição de normalidade dos resíduos.

Daniel (1960) propôs um método manual para detectar e corrigir um valor discrepante ou atípico em um planejamneto \(2^k\) não replicado. Este método consiste em três etapas:

Primeiro, a presença de um valor discrepante é detectada por uma lacuna no centro de um gráfico normal de efeitos.

Segundo, o outlier é identificado combinando os sinais dos efeitos insignificantes com os sinais dos níveis dos fatores codificados e das interações de cada observação.

Terceiro, estimar a magnitude da discrepância e corrigir o valor atípico.

Como exemplo, considere o gráfico normal de efeitos (Figura 3.19) de um experimento \(2^4\) não replicado descrito por Box (1991). Neste gráfico parece que os efeitos principais B e C podem ser significativos, mas há uma lacuna vertical na linha de efeitos insignificantes que indica que um valor discrepante pode estar presente.

Lawson and Gatlin (2006) automatizaram o procedimento de Daniel identificando e corrigindo um valor atípico fazendo duas passagens pelos dados. Se a estatística do gap (a razão do gap vertical na Figura 3.19 dividida pela estatística PSE de Lenth) estiver acima do percentil 50 de sua distribuição de referência na primeira passagem pelos dados, o PSE é recalculado após a correção do valor discrepante e a estatística do gap é testada. novamente em uma segunda passagem pelos dados.

Se a estatística de lacuna estiver acima do percentil 95 na segunda passagem pelos dados, a função imprime uma tabela mostrando qual observação é o valor atípico potencial e quais efeitos são significativos após a correção do valor atípico. Um gráfico semi-normal é usado para identificar os efeitos significativos após a correção do outlier.

Uma função R chamada Gaptest para realizar este procedimento está disponível no pacote daewr. A chamada de função é mostrada abaixo e os resultados estão abaixo dela. O conjunto de dados BoxM na chamada da instrução ou comando contém os dados de Box (1991). As primeiras quatro colunas no conjunto de dados são os fatores A, B, C e D no experimento \(2^4\) não replicado e a última coluna é uma resposta numérica \(y\).

# Figure 3.19 (without Gap annotation)
modB = lm(y ~ A*B*C*D, data = BoxM) 
fullnormal(coef(modB)[-1],alpha=.2)

Figura 3.19: Gráfico normal de efeitos com lacuna vertical.

O método detectou um outlier na 13ª execução e corrigiu o resultado da resposta alterando 59.15 para 52.15. A reanálise dos dados corrigidos mostra que os efeitos principais B, C e a interação AxC são significativos.

A interação AxC não pôde ser detectada na Figura 3.19 porque o valor discrepante inflou as estimativas dos efeitos não significativos e o efeito AxC foi enterrado no ruído.

library(daewr)
data(BoxM)
par ( mar = c(4,4,2,1), pch = 19 )
Gaptest(BoxM)

## Effect Report 
##    
## Label     Half Effect    Sig(.05) 
## A          -0.400        no         
## B          -2.110        no         
## C           1.855        no         
## D           0.505        no         
## AB          0.455        no         
## AC         -1.245        no         
## AD         -0.290        no         
## BC         -0.400        no         
## BD         -0.590        no         
## CD          0.745        no         
## ABC         0.600        no         
## ABD         0.360        no         
## ACD         0.200        no         
## BCD        -0.790        no         
## ABCD        0.760        no         
##    
## Lawson, Grimshaw & Burt Rn Statistic =  1 
## 95th percentile of Rn =  1.201 
## Initial Outlier Report 
## Standardized-Gap =  3.3532 Significant at 50th percentile 
## Final Outlier Report 
## Standardized-Gap =  13.189 Significant at 99th percentile 
##     
##     Corrrected Data Report   
## Response  Corrected Response   Detect Outlier 
##    47.46         47.46             no         
##    49.62         49.62             no         
##    43.13         43.13             no         
##    46.31         46.31             no         
##    51.47         51.47             no         
##    48.49         48.49             no         
##    49.34         49.34             no         
##    46.10         46.10             no         
##    46.76         46.76             no         
##    48.56         48.56             no         
##    44.83         44.83             no         
##    44.45         44.45             no         
##    59.15         52.75             yes        
##    51.33         51.33             no         
##    47.02         47.02             no         
##    47.90         47.90             no

## Effect Report 
##    
## Label     Half Effect    Sig(.05) 
## A        -4.5143e-15        no         
## B        -1.7100e+00        yes        
## C         1.4550e+00        yes        
## D         1.0500e-01        no         
## AB        5.5000e-02        no         
## AC       -8.4500e-01        yes        
## AD        1.1000e-01        no         
## BC        2.1316e-15        no         
## BD       -1.9000e-01        no         
## CD        3.4500e-01        no         
## ABC       2.0000e-01        no         
## ABD      -4.0000e-02        no         
## ACD       6.0000e-01        no         
## BCD      -3.9000e-01        no         
## ABCD      3.6000e-01        no         
##    
## Lawson, Grimshaw & Burt Rn Statistic =  1.6261 
## 95th percentile of Rn =  1.201

Figura 3.20: Gráfico meio normal de coeficientes calculados com dados corrigidos.

A Figura 3.20 mostra o gráfico semi-normal dos coeficientes calculados com os dados corrigidos. Neste gráfico pode-se observar que a interação AxC é claramente significativa além dos efeitos principais B e C que foram identificados na Figura 3.19.

Sempre que um valor discrepante for descoberto, usando este método ou gráficos de resíduos e as conclusões da análise mudarem quando o valor discrepante for removido ou corrigido, o experimentador deve proceder com cautela. Quando há mais de duas réplicas nas mesmas configurações de fator onde o valor discrepante foi encontrado, pode ficar claro que algo está errado. No entanto, se houver duas ou menos observações em configurações de fatores onde o valor discrepante foi encontrado, pode ser aconselhável reexecutar o experimento questionável.


3.9 Exercícios


1- Um consultor foi chamado para auxiliar o departamento de polícia de uma grande cidade metropolitana na avaliação do curso de relações humanas para novos policiais. Ele planejou um experimento fatorial de dois fatores em que os tratamentos eram A — o tipo de ronda para o qual os policiais eram designados, e B — a duração do curso de relações humanas. Foi escolhida uma amostra de 45 novos oficiais, e 5 foram designados aleatoriamente para cada uma das 9 combinações de tratamento. Foi desenvolvido um teste para medir a atitude dos agentes em relação aos grupos minoritários e foi aplicado aos agentes participantes após a sua formação ter terminado e terem passado duas semanas na sua ronda. Melhores atitudes resultam em pontuações mais altas neste teste. A análise dos dados revelou um efeito significativo de interação A×B entre o tipo de batida e a duração do curso de relações humanas. A tabela abaixo mostra as pontuações médias dos testes para as 9 combinações de níveis de tratamento.

\[ \begin{array}{c|c} & \mbox{Duração do Curso de Relações Humanas} \\ \begin{array}{l} \mbox{Tipo de batida}\\\hline \mbox{batida da classe alta} \\ \mbox{batida de classe média} \\ \mbox{batida do centro da cidade} \\\hline \end{array} & \begin{array}{ccc} \mbox{5 horas} & \mbox{10 horas} & \mbox{15 horas} \\\hline 34.4 & 35.5 & 39.2 \\ 30.2 & 32.4& 34.7 \\ 20.1 & 39.4 & 54.3 \\\hline \end{array} \end{array} \]

  1. Construa um gráfico de interação.

  2. Escreva uma interpretação da interação em algumas frases completas.

2- Uma catapulta de madeira pode ser usada para arremessar uma bola de espuma. A catapulta possui três fatores que podem ser ajustados: o ângulo inicial (start angle), o ângulo de parada (stop angle) e a altura do pivô (pivot height). A distância percorrida pela bola pode ser medida com uma fita métrica.

  1. Se fossem realizados experimentos com a catapulta lançando a bola e medindo a distância, qual seria a unidade experimental?

  2. Usando os números 1, 2 e 3 para representar os níveis do ângulo inicial e do ângulo final, e mantendo a altura do pivô constante em seu nível alto, faça uma lista aleatória de experimentos para um experimento fatorial 3×3 com \(r = 2\) réplicas por célula.

  3. Se a variância do erro experimental na distância medida fosse \(\sigma^2 = 12\) polegadas, calcule o número de repetições que você precisaria para ter um poder de 0.90 para detectar uma diferença de 10 polegadas nas médias das células.

  4. Calcule o número de réplicas necessárias para ter um poder de 0.90 para detectar uma diferença de 24 polegadas nas médias marginais para qualquer um dos fatores.

  5. Se uma catapulta estiver disponível, realize a lista de experimentos que você escreveu na parte (b).

  6. Calcule a ANOVA com os dados resultantes e teste os principais efeitos e interação.

  7. Explique ou interprete quaisquer efeitos significativos, use gráficos se necessário para ajudar na sua explicação.

3- Em um experimento para maximizar a resolução \(Y\) = de um pico em um cromatógrafo a gás, foi encontrada uma interação significativa entre \(A\) = temperatura da coluna e \(C\) = taxa de fluxo de gás. A tabela abaixo mostra a resolução média em cada combinação de temperatura da coluna e vazão de gás.

\[ \begin{array}{c|cc} \mbox{Temperatura} & \mbox{Taxa de} & \mbox {fluxo de gás} \\ \mbox{da coluna} & \mbox{Baixo} & \mbox{Alto} \\\hline 120 & 10 & 13 \\ 180 & 12 & 18 \\\hline \end{array} \] (a) Construa um gráfico de interação.

  1. Escreva uma ou duas frases para interpretar essa interação.

4- Considere realizar experimentos para determinar o efeito da marca da pipoca, do nível de potência e do tempo na porcentagem de grãos de pipoca comestíveis feitos em um forno de micro-ondas, observando que àqueles não estourados ou queimados são não comestíveis. O objetivo é maximizar a proporção de grãos comestíveis. Comece com 13 xícaras de grãos e faça um experimento piloto para determinar a gama de fatores que você gostaria de estudar e forneça uma estimativa do desvio padrão das réplicas feitas nas mesmas condições.

  1. Qual é a unidade experimental deste experimento?

  2. Determine quantas réplicas serão necessárias para detectar uma diferença máxima nas médias marginais, efeito principal, de 0.25.

  3. Determine o número de níveis dos fatores que você gostaria de estudar e crie uma lista aleatória dos experimentos para um planejamento fatorial.

  4. Realize realmente os experimentos e colete os dados.

  5. Analisar os dados. Divida quaisquer efeitos significativos de nível de potência e tempo em contrastes polinomiais ortogonais. Interprete os resultados e determine a combinação ideal de níveis de fatores dentre aqueles que você testou.

  6. Quantos experimentos seriam necessários para obter o mesmo poder de detecção dos efeitos principais usando um planejamento de experimentos variando um fator por vez? Você detectaria o mesmo ótimo usando o plano de variação de um fator por vez?

5- Modifique o código R na seção 3.5.2 para verificar que o poder seria 0.80 para detectar diferenças nas médias marginais de \(\Delta = 1.0\) polegada nos experimentos de helicóptero de papel com quatro níveis de comprimento de asa, quatro níveis de largura do corpo e \(r = 2\) replicações por célula, conforme mostrado na Seção 3.5.2.

6- Numa continuação das experiências para investigar o efeito dos reguladores de crescimento das plantas e das escamas dos botões de lança no alongamento da lança em espargos, descritas no exercício 5 do Capítulo 2, Yang-Gyu and Woolley (2006) conduziram uma experiência factorial 4×3 variando a concentração do regulador de crescimento vegetal CPPU na solução e o tempo em que as lanças de aspargos foram mergulhadas na solução. Os resultados são mostrados na tabela abaixo.

\[ \begin{array}{c|c} & \mbox{Tempo de imersão (seg)} \\ \mbox{CPPU Conc.} & \begin{array}{ccc} 30 \quad & 60 & \; 90 \end{array} \\\hline 0 \mbox{ control} & \begin{array}{ccc} 92.5 & 92.9 & 91.3 \end{array} \\ 0.5 \mbox{ ppm} & \begin{array}{ccc} 97.8 & 94.9 & 101.3 \end{array} \\ 1.0 \mbox{ ppm} & \begin{array}{ccc} 97.0 & 98.5 & 101.6 \end{array} \\ 10.0 \mbox{ ppm} & \begin{array}{ccc} 103.4 & 102.9 & 98.6 \end{array} \\ \end{array} \]

  1. Divida as somas dos quadrados para a concentração de CPPU e tempo de imersão em contrastes polinomiais ortogonais. Particione a interação do CPPU conc. mergulhando o tempo em linear×linear, quadrático×linear, cúbico×linear, linear×quadrático, quadrático×quadrático e cúbico×quadrático. Agrupe (soma) todas as somas dos quadrados, exceto a porção linear×linear, para formar uma soma errada dos quadrados e use-a para testar a parte linear×linear da interação e os contrastes polinomiais encontrados para os efeitos principais. Interprete seus resultados.

  2. Use a função Tukey1df para testar a significância da interação.

  3. Com base nos resultados obtidos em (a) e (b), você recomendaria o modelo aditivo para esses dados?

7- Kenett and Steinberg (1987) descreveram um experimento fatorial de dois níveis conduzido por estudantes para estudar o tempo necessário para ferver 1 litro de água. Os fatores foram A = nível de chama (baixo ou alto), B = tamanho da panela (pequena ou grande), C = tampa da panela (nenhuma ou tampa de vidro) e D = sal adicionado à água (não ou sim).

  1. Se o desvio padrão no tempo de ebulição, testado nas mesmas condições, for \(\widehat{\sigma}\) = 0.236 minutos, use a fórmula de atalho para determinar quantos experimentos você precisará realizar para ter um poder de 0.95 para detectar efeitos de tamanho \(\Delta\)=0.50 minutos. Essa resposta mudaria se você decidisse realizar um experimento apenas com 3 dos 4 fatores?

  2. Crie uma lista de experimentos em ordem aleatória para realizá-los.

  3. Realize realmente os experimentos fervendo a água e colete os dados.

  4. Analisar os dados para determinar quais efeitos e interações são significativos.

  5. Interprete e explique quaisquer efeitos significativos que você encontrar.

8- Considere os dados da Tabela 3.7 com a terceira observação (90) substituída por um valor atípico (78) para essas configurações de fator.

  1. Calcule os efeitos para esses dados.

  2. Faça um gráfico normal ou semi-normal dos efeitos.

  3. Que efeitos você considera significativos?

  4. Execute a função Gaptest para detectar e corrigir um outlier com esses dados.

9- Nyberg (1999) mostrou que o nitreto de silício (SiNx) cultivado por Deposição Química de Vapor Aprimorada por Plasma (PECVD) é um candidato promissor para um revestimento anti-reflexo (ARC) em células solares comerciais de silício cristalino. O nitreto de silício foi cultivado em pastilhas de silício 4A polidas (100) orientadas usando um reator PECVD da Plasma Technology de placa paralela. O diâmetro dos eletrodos do PECVD é de 24 cm e o diâmetro do chuveiro, e entram os gases, é de 2A. A frequência de RF foi de 13.56 MHz. A espessura do nitreto de silício era um quarto do comprimento de onda da luz no nitreto, sendo o comprimento de onda de 640 nm. Espera-se que este comprimento de onda esteja próximo do ideal para fins de células solares de silício. Os gases do processo foram amônia e uma mistura de 3% de silano em argônio. Os experimentos foram realizados de acordo com um planejamento fatorial \(2^5\). Os resultados são mostrados na tabela abaixo.

  1. Ajuste o modelo fatorial à resposta \(y_1\) incluindo todas as interações até ordem 5.

  2. Faça um gráfico normal dos efeitos ou coeficientes de regressão para determinar quais efeitos principais e interações são significativos.

  3. Elimine termos não significativos do modelo e faça gráficos de resíduos para verificar as premissas de ajuste do modelo.

  4. Use seu modelo para determinar as condições, ou seja, níveis de fator que minimizarão o índice de refração.

\[ \begin{array}{cccccccc} \mbox{No.} & A & B & C & D & E & y_1 & y_2 \\\hline 1 & 0.1 & 40 & 300 & 300 & 10 & 1.92 & 1.79 \\ 2 & 0.9 & 40 & 300 & 300 & 10 & 3.06 & 10.1 \\ 3 & 0.1 & 220 & 300 & 300 & 10 & 1.96 & 3.02 \\ 4 & 0.9 & 220 & 300 & 300 & 10 & 3.33 & 15 \\ 5 & 0.1 & 40 & 1200 & 300 & 10 & 1.87 & 19.7 \\ 6 & 0.9 & 40 & 1200 & 300 & 10 & 2.62 & 11.2 \\ 7 & 0.1 & 220 & 1200 & 300 & 10 & 1.92 & 35.7 \\ 8 & 0.9 & 220 & 1200 & 300 & 10 & 2.96 & 36.2 \\ 9 & 0.1 & 40 & 300 & 460 & 10 & 1.94 & 2.31 \\ 10 & 0.9 & 40 & 300 & 460 & 10 & 3.53 & 5.58 \\ 11 & 0.1 & 220 & 300 & 460 & 10 & 2.06 & 2.75 \\ 12 & 0.9 & 220 & 300 & 460 & 10 & 3.75 & 14.5 \\ 13 & 0.1 & 40 & 1200 & 460 & 10 & 1.96 & 20.7 \\ 14 & 0.9 & 40 & 1200 & 460 & 10 & 3.14 & 11.7 \\ 15 & 0.1 & 220 & 1200 & 460 & 10 & 2.15 & 31 \\ 16 & 0.9 & 220 & 1200 & 460 & 10 & 3.43 & 39 \\ 17 & 0.1 & 40 & 300 & 300 & 60 & 1.95 & 3.93 \\ 18 & 0.9 & 40 & 300 & 300 & 60 & 3.16 & 12.4 \\ 19 & 0.1 & 220 & 300 & 300 & 60 & 2.01 & 6.33 \\ 20 & 0.9 & 220 & 300 & 300 & 60 & 3.43 & 23.7 \\ 21 & 0.1 & 40 & 1200 & 300 & 60 & 1.88 & 35.3 \\ 22 & 0.9 & 40 & 1200 & 300 & 60 & 2.14 & 15.1 \\ 23 & 0.1 & 220 & 1200 & 300 & 60 & 1.98 & 57.1 \\ 24 & 0.9 & 220 & 1200 & 300 & 60 & 2.81 & 45.9 \\ 25 & 0.1 & 40 & 300 & 460 & 60 & 1.97 & 5.27 \\ 26 & 0.9 & 40 & 300 & 460 & 60 & 3.67 & 12.3 \\ 27 & 0.1 & 220 & 300 & 460 & 60 & 2.09 & 6.39 \\ 28 & 0.9 & 220 & 300 & 460 & 60 & 3.73 & 30.5 \\ 29 & 0.1 & 40 & 1200 & 460 & 60 & 1.98 & 30.1 \\ 30 & 0.9 & 40 & 1200 & 460 & 60 & 2.99 & 14.5 \\ 31 & 0.1 & 220 & 1200 & 460 & 60 & 2.19 & 50.3 \\ 32 & 0.9 & 220 & 1200 & 460 & 60 & 3.39 & 47.1 \\ \end{array} \]

A = Taxa de fluxo de silano para amônia

B = Taxa total de fluxo de gás (sccm)

C = Pressão (mtorr)

D = Temperatura (℃)

E = Potência (W)

\(y_1\) = Índice de refração

\(y_2\) = Taxa de crescimento (nm/min)