CE 074 - Controle de Processos Industriais

http://www.leg.ufpr.br/ce074


Universidade Federal do Paraná Prof. Fernando de Pol Mayer

Curso de Graduação em Estatística Lab. de Estatística e Geoinformação - LEG

Departamento de Estatística - UFPR


Experimentos fatoriais fracionados

Nos experimentos fatoriais \(2^k\), o número de corridas cresce rapidamente com o número de fatores.

k 2 3  4  5  6   7   8   9   10
n 4 8 16 32 64 128 256 512 1024

Com isso, se existe necessidade de estudar muitos fatores, esse tipo de planejamento oferece limitação devido ao elevado número de corridas exigidas na sua execução. Ao correr um experimento \(2^k\), é possível estimar efeitos principais, interações duplas e interações até de grau \(k\). Porém, essas interações de ordem elevada podem ser desconsideradas, já que normalmente não são de interesse prático e muitas vezes possuem efeitos desprezíveis. O que se considera mais importante são os efeitos principais e interações baixas (duplas, triplas). Assim, correr o experimento \(2^k\) iria permitir estimar interações altas que de fato não seriam aproveitadas. Então o ponto é: não seria possível reduzir o número de corridas do experimento já qua há pouco interesse nas interações altas? Sim. É exatamente esse o elemento central dos experimentos fatoriais fracionados.

Os fatoriais fracionários tem um uso importante nos experimentos exploratórios (screening experiments). Estes são experimentos em que muitos fatores são considerados, com a finalidade de se identificar aqueles fatores que possuem efeitos grandes. Experimentos exploratórios geralmente são realizados nos estágios iniciais de um projeto, quando é provavel que muitos dos fatores considerados tenham pouco ou nenhum efeito na resposta. Os fatores que forem identificados como importantes serão então investigados mais a fundo em experimentos subsequentes.

A idéia básica é dividir as corridas do fatorial \(2^k\) em frações, exatamente da mesma forma que fazemos para aplicar blocagem. Na blocagem escolhemos um efeito para propositalmente confundirmos com o efeito dos blocos. Nos fatoriais fracionados é similar: define-se efeitos que serão propositalmente confundidos. O confundimento entre efeitos é o preço que se paga para reduzir o número de corridas. Perceba que nosso objetivo é estudar \(k\) fatores onde seria possível estudar com todas as corridas um fatorial com k* < k. O ponto chave ao delinear o experimento é saber escolher os efeitos a serem confundidos para tirar a maior vantagem possível do experimento que é inferir sobre os efeitos principais e interações baixas.

Fatoriais fracionados \(2^{k-1}\), ou meia-fração de um fatorial \(2^k\)

Nesses experimentos avaliamos k fatores fazendo metade das corridas previstas para o fatorial \(2^k\), ou seja \(2^{k-1} = 2^k 2^{-1} = \frac{2^k}{2}\). Para saber que níveis de cada fator que serão corridos, temos que montar a matriz de delineamento do experimento com \(k-1\) fatores. Os níveis do \(k\)-ésimo fator são determinados a partir dos precedentes e do efeito escolhido para o confundimento.

Considere como exemplo, um experimento fatorial \(2^{3-1}\), ou seja, uma meia-fração de um \(2^3\). Esse planejamento possui apenas quatro corridas, em contraste com as oito corridas do planejamento original.

A tabela de siais para o planejamento \(2^3\) é mostrada abaixo.

    I  A  B  C AB AC BC ABC
a   1  1 -1 -1 -1 -1  1   1
b   1 -1  1 -1 -1  1 -1   1
c   1 -1 -1  1  1 -1 -1   1
abc 1  1  1  1  1  1  1   1
(1) 1 -1 -1 -1  1  1  1  -1
ab  1  1  1 -1  1 -1 -1  -1
ac  1  1 -1  1 -1  1 -1  -1
bc  1 -1  1  1 -1 -1  1  -1

Repare que esta tabela está propositalmente ordenada em ordem decrescente pela coluna ABC. Se selecionarmos as quatro combinações de tratamento onde a coluna ABC possui sinal positivo, ou seja, \(a\), \(b\), \(c\), \(abc\), estamos selecionando uma meia-fração desse planejamento. Sendo assim, dizemos que \(ABC\) é o gerador dessa fração particular. Além disso, o elemento identidade \(I\) possui também o mesmo sinal (positivo) para as quatro corridas. Assim, chamamos

\[ I = ABC \]

de relação de definção para esse planejamento fracionário. Essa fração com sinal positivo na relação de definição é chamada de fração principal, enquanto que a fração com sinal negativo em ABC, \(I = -ABC\) é chamada de fração alternada.

Ainda analisando a parte superior da tabela acima, note que os sinais que estimam os contrastes para o fator \(A\) são os mesmos que estimam os contrastes para a interação \(BC\). Da mesma forma, os contrastes de \(B\) são os mesmos de \(AC\), e assim por diante. Isso ocorre pois a relação de definição gera uma série de confundimentos no fatorial fracionário. De fato, todos os efeitos e interações estão confundidos com algum outro efeito ou interação. Se multiplicarmos cada fator ou interação pela relação de definição, e usarmos as propriedades da tabela de sinais, podemos ver que:

\[ \begin{align} A \cdot I = A \cdot ABC = A^2BC = BC \\ B \cdot I = B \cdot ABC = AB^2C = AC \\ C \cdot I = C \cdot ABC = ABC^2 = AB \\ \end{align} \]

Portanto, só precisamos estimar os efeitos de A, B, e C, pois BC, AC, e AB estão confundidos com estes efeitos. À estes efeitos que estão confundidos chamamos de pares associados, ou aliases. Em experimentos fatoriais fracionados \(2^{k-1}\) todo efeito possui um par associado. Nesse caso específico de um fracionário \(2^{3-1}\), cada efeito principal tem um par associado de segunda ordem.

Repare então que ao estimar o afeito de A, por exemplo, também estamos estimando o efeito de BC. Se este efeito for grande, não temos como saber se ele é causado por A ou por BC isoladamente, pois eles estão confundidos. No entanto, como vimos anteriormente, é muito mais comum que efeitos principais sejam mais importantes do que efeitos de ordens mais altas. Assim, se assumirmos que efeitos de ordens mais altas tem pouco ou nenhum efeito, então podemos concluir que um efeito grande de A seja exclusivamente desse efeito principal. O mesmo vale para os outros efeitos.

Novamente aqui cabe uma frase já escrita anteriormente: o confundimento entre efeitos é o preço que se paga para reduzir o número de corridas. A ideia do experimento fatorial fracionário é de justamente correr experimentos com muitos fatores em frações menores, verificar os efeitos que realmente importam, e daí sim planejar fatoriais completos, mas com menos fatores, aqueles que realmente importam. Dessa forma, podemos correr uma sequência de experimentos pequenos e eficientes, combinar informações por meio de vários experimentos, e tirar vantagem de aprender sobre o processo que estamos experimentando à medida que continuamos. Esse é um bom exemplo do conceito de experiência sequencial.

Um planejamento \(2^{k-1}\) pode ser construído escrevendo as combinações dos tratamentos para um fatorial completo com \(k-1\) fatores, chamado de planejamento básico. Depois, adiciona-se o \(k\)-ésimo fator, identificando seus níveis alto e baixo com os sinais mais e menos da interação de mais alta ordem.

Por exemplo, umm fatorial fracionário \(2^{4-1}\) é construído escrevendo o planejamento básico \(2^3\), e então igualando o fator \(D\) com a interação \(ABC\), pois a relação de definição é \(I = ABCD\) e \(D \cdot I = D \cdot ABCD = ABC\).

##======================================================================
## Exemplo 2^(4-1), tabela de sinais do 2^3
da <- expand.grid(A = c(-1, 1), B = c(-1, 1), C = c(-1, 1))
da
   A  B  C
1 -1 -1 -1
2  1 -1 -1
3 -1  1 -1
4  1  1 -1
5 -1 -1  1
6  1 -1  1
7 -1  1  1
8  1  1  1
## Níveis de D gerados pelo efeito ABC
da$D <- with(da, A*B*C)
da
   A  B  C  D
1 -1 -1 -1 -1
2  1 -1 -1  1
3 -1  1 -1  1
4  1  1 -1 -1
5 -1 -1  1  1
6  1 -1  1 -1
7 -1  1  1 -1
8  1  1  1  1

Veja que propositalmente atribuímos os níveis de D sobre o efeito ABC, ou seja, confundimos estes dois efeitos. Mas sabemos que este não é o único par de efeitos que está confundido. Se obtermos a tabela de sinais para todos os efeitos veremos mais efeitos confundidos.

(tab <- model.matrix(~ A*B*C*D, da))
  (Intercept)  A  B  C  D A:B A:C B:C A:D B:D C:D A:B:C A:B:D A:C:D B:C:D
1           1 -1 -1 -1 -1   1   1   1   1   1   1    -1    -1    -1    -1
2           1  1 -1 -1  1  -1  -1   1   1  -1  -1     1    -1    -1     1
3           1 -1  1 -1  1  -1   1  -1  -1   1  -1     1    -1     1    -1
4           1  1  1 -1 -1   1  -1  -1  -1  -1   1    -1    -1     1     1
5           1 -1 -1  1  1   1  -1  -1  -1  -1   1     1     1    -1    -1
6           1  1 -1  1 -1  -1   1  -1  -1   1  -1    -1     1    -1     1
7           1 -1  1  1 -1  -1  -1   1   1  -1  -1    -1     1     1    -1
8           1  1  1  1  1   1   1   1   1   1   1     1     1     1     1
  A:B:C:D
1       1
2       1
3       1
4       1
5       1
6       1
7       1
8       1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

Os pares de efeitos \((A, BCD), (B, ACD), \ldots, (A:B, C:D)\) estão confundidos, pois têm os mesmos sinais. Verificar aos pares as colunas para encontrar os pares associados é demorado, por isso é mais fácil fazer operações com o contraste de definição, \(I = ABCD\). Assim, para sabermos quem está confundido com quem, podemos fazer

\[ \begin{align} A \cdot I = A \cdot ABCD = A^2BCD = BCD \\ B \cdot I = B \cdot ABCD = AB^2CD = ACD \\ C \cdot I = C \cdot ABCD = ABC^2D = ABD \\ D \cdot I = D \cdot ABCD = ABCD^2 = ABC \\ \end{align} \]

para os efeitos principais, e

\[ \begin{align} AB \cdot I = AB \cdot ABCD = A^2B^2CD = CD \\ AC \cdot I = AC \cdot ABCD = A^2BC^2D = BD \\ AD \cdot I = AD \cdot ABCD = A^2BCD^2 = BC \\ \end{align} \]

para os efeitos de segunda ordem. Veja que aqui, os efeitos principais estão confundidos com efeitos de terceira ordem, e efeitos de segunda ordem estão confundidos entre si.

No R, podemos obter a tabela de sinais para os efeitos únicos através de

(tab <- unique(tab, MARGIN = 2))
  (Intercept)  A  B  C  D A:B A:C B:C
1           1 -1 -1 -1 -1   1   1   1
2           1  1 -1 -1  1  -1  -1   1
3           1 -1  1 -1  1  -1   1  -1
4           1  1  1 -1 -1   1  -1  -1
5           1 -1 -1  1  1   1  -1  -1
6           1  1 -1  1 -1  -1   1  -1
7           1 -1  1  1 -1  -1  -1   1
8           1  1  1  1  1   1   1   1

Veja que ao usar esta tabela para estimar os efeitos de I, A, B, C, D, AB, AC, e BC, também estamos estiamdo BCD, ACD, ABD, ABC, CD, BD, e AD, respectivamente (efeitos confundidos).

Podemos escolher qualquer contraste para confundir, no entanto, usar as interações mais altas fazem com que os efeitos principais fiquem confundidos com estas, que podem ser assumidas como nulas. Então que o valor observado do efeito pode ser considerado como exclusivamente do efeito principal.

Uma outra forma de se construir um planejamento \(2^{k-1}\) é usar a própria relação de definição para “fracionar” o experimento original. Nesse mesmo exemplo do fatorial \(2^{4-1}\), podemos então obter o experimento completo \(2^4\)

db <- expand.grid(A = c(-1, 1), B = c(-1, 1),
                  C = c(-1, 1), D = c(-1, 1))
db
    A  B  C  D
1  -1 -1 -1 -1
2   1 -1 -1 -1
3  -1  1 -1 -1
4   1  1 -1 -1
5  -1 -1  1 -1
6   1 -1  1 -1
7  -1  1  1 -1
8   1  1  1 -1
9  -1 -1 -1  1
10  1 -1 -1  1
11 -1  1 -1  1
12  1  1 -1  1
13 -1 -1  1  1
14  1 -1  1  1
15 -1  1  1  1
16  1  1  1  1

Como a relação de definição é \(I = ABCD\), e ela é a geradora da meia-fração do experimento, então so precisamos calcular \(ABCD\)

db$ABCD <- with(db, A*B*C*D)
db
    A  B  C  D ABCD
1  -1 -1 -1 -1    1
2   1 -1 -1 -1   -1
3  -1  1 -1 -1   -1
4   1  1 -1 -1    1
5  -1 -1  1 -1   -1
6   1 -1  1 -1    1
7  -1  1  1 -1    1
8   1  1  1 -1   -1
9  -1 -1 -1  1   -1
10  1 -1 -1  1    1
11 -1  1 -1  1    1
12  1  1 -1  1   -1
13 -1 -1  1  1    1
14  1 -1  1  1   -1
15 -1  1  1  1   -1
16  1  1  1  1    1
db[order(db$ABCD, decreasing = TRUE), ]
    A  B  C  D ABCD
1  -1 -1 -1 -1    1
4   1  1 -1 -1    1
6   1 -1  1 -1    1
7  -1  1  1 -1    1
10  1 -1 -1  1    1
11 -1  1 -1  1    1
13 -1 -1  1  1    1
16  1  1  1  1    1
2   1 -1 -1 -1   -1
3  -1  1 -1 -1   -1
5  -1 -1  1 -1   -1
8   1  1  1 -1   -1
9  -1 -1 -1  1   -1
12  1  1 -1  1   -1
14  1 -1  1  1   -1
15 -1  1  1  1   -1

e selecionar a fração principal de \(ABCD\), ou seja, aquela com sinal positivo em \(ABCD\),

db <- subset(db, ABCD == 1)
db
    A  B  C  D ABCD
1  -1 -1 -1 -1    1
4   1  1 -1 -1    1
6   1 -1  1 -1    1
7  -1  1  1 -1    1
10  1 -1 -1  1    1
11 -1  1 -1  1    1
13 -1 -1  1  1    1
16  1  1  1  1    1

Veja que automaticamente ficamos com a relação anterior de \(D = ABC\)

db$D == with(db, A*B*C)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Portanto as duas formas mostradas acima podem ser utilizadas para construir qualquer experimento fatorial fracionado \(2^{k-1}\).

Exemplo de um experimento \(2^4\), convertido para um experimento \(2^{4-1}\)

Exemplo 14-5, Montgomery (EAPE), pg. 385: processo de ataque químico localizado sobre nitreto, por meio de uma sonda de plasma de pastilha única. Variável resposta: taxa de ataque do nitreto de silício. Variáveis explicativas:

                Fatores do experimento
 Nível     Espaçamento   Pressão    Vazão de C2F6   Potência
              (cm)       (m Torr)   (cm³ pad/min)     (W)
 Baixo (-)    0.80         450           125          275
 Alto  (+)    1.20         550           200          325

Este exemplo já foi analsado na aula 05. Este é um experimento fatorial \(2^4\) completo, ou seja, foram corridas todas as 16 combinações de tratamentos possíveis. Vamos re-analisar estes dados completos aqui, e depois realizar o fracionamento desse experimento, com aintenção de verificar se chegamos à mesma concusão.

##----------------------------------------------------------------------
## Importação
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-5.txt"
dados <- read.table(url, header = TRUE)
row.names(dados) <- apply(dados, 1,
                          function(i) paste(letters[1:4][i==1], collapse = ""))
row.names(dados)[1] <- "(1)"
dados
      A  B  C  D    y
(1)  -1 -1 -1 -1  550
a     1 -1 -1 -1  669
b    -1  1 -1 -1  604
ab    1  1 -1 -1  650
c    -1 -1  1 -1  633
ac    1 -1  1 -1  642
bc   -1  1  1 -1  601
abc   1  1  1 -1  635
d    -1 -1 -1  1 1037
ad    1 -1 -1  1  749
bd   -1  1 -1  1 1052
abd   1  1 -1  1  868
cd   -1 -1  1  1 1075
acd   1 -1  1  1  860
bcd  -1  1  1  1 1063
abcd  1  1  1  1  729
## Definições do experimento
k <- 4
n <- 1

##----------------------------------------------------------------------
## Montando a tabela de sinais
(tab <- model.matrix(~A * B * C * D, data = dados))
     (Intercept)  A  B  C  D A:B A:C B:C A:D B:D C:D A:B:C A:B:D A:C:D
(1)            1 -1 -1 -1 -1   1   1   1   1   1   1    -1    -1    -1
a              1  1 -1 -1 -1  -1  -1   1  -1   1   1     1     1     1
b              1 -1  1 -1 -1  -1   1  -1   1  -1   1     1     1    -1
ab             1  1  1 -1 -1   1  -1  -1  -1  -1   1    -1    -1     1
c              1 -1 -1  1 -1   1  -1  -1   1   1  -1     1    -1     1
ac             1  1 -1  1 -1  -1   1  -1  -1   1  -1    -1     1    -1
bc             1 -1  1  1 -1  -1  -1   1   1  -1  -1    -1     1     1
abc            1  1  1  1 -1   1   1   1  -1  -1  -1     1    -1    -1
d              1 -1 -1 -1  1   1   1   1  -1  -1  -1    -1     1     1
ad             1  1 -1 -1  1  -1  -1   1   1  -1  -1     1    -1    -1
bd             1 -1  1 -1  1  -1   1  -1  -1   1  -1     1    -1     1
abd            1  1  1 -1  1   1  -1  -1   1   1  -1    -1     1    -1
cd             1 -1 -1  1  1   1  -1  -1  -1  -1   1     1     1    -1
acd            1  1 -1  1  1  -1   1  -1   1  -1   1    -1    -1     1
bcd            1 -1  1  1  1  -1  -1   1  -1   1   1    -1    -1    -1
abcd           1  1  1  1  1   1   1   1   1   1   1     1     1     1
     B:C:D A:B:C:D
(1)     -1       1
a       -1      -1
b        1      -1
ab       1       1
c        1      -1
ac       1       1
bc      -1       1
abc     -1      -1
d        1      -1
ad       1       1
bd      -1       1
abd     -1      -1
cd      -1       1
acd     -1      -1
bcd      1      -1
abcd     1       1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
##----------------------------------------------------------------------
## Estimativa dos efeitos
## Podemos usar as colunas da tabela de sinais para calcular os
## contrastes
y <- dados$y
## Contrastes
contr <- t(tab[, -1]) %*% y
## Efeitos = contraste/(n2^{k-1})
(ef <- contr/(n * 2^(k - 1)))
            [,1]
A       -101.625
B         -1.625
C          7.375
D        306.125
A:B       -7.875
A:C      -24.875
B:C      -43.875
A:D     -153.625
B:D       -0.625
C:D       -2.125
A:B:C    -15.625
A:B:D      4.125
A:C:D      5.625
B:C:D    -25.375
A:B:C:D  -40.125
##----------------------------------------------------------------------
## Gráfico de probabilidade normal dos efeitos estimados
qqaux <- qqnorm(ef, col = 2, pch = 19); qqline(ef)
text(qqaux$x, qqaux$y, rownames(qqaux$y), cex = 0.8, pos = 3)

##----------------------------------------------------------------------
## Ajuste do modelo
## Como vimos no gráfico anterior, as interações de ordem maior do que 2
## são desnecessárias, portanto vamos ajustar o modelo com interações de
## até segunda ordem apenas
m0 <- lm(y ~ (A + B + C + D)^2, data = dados)
anova(m0)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq  F value    Pr(>F)    
A          1  41311   41311  20.2765  0.006382 ** 
B          1     11      11   0.0052  0.945391    
C          1    218     218   0.1068  0.757069    
D          1 374850  374850 183.9879 3.903e-05 ***
A:B        1    248     248   0.1218  0.741351    
A:C        1   2475    2475   1.2148  0.320582    
A:D        1  94403   94403  46.3357  0.001042 ** 
B:C        1   7700    7700   3.7794  0.109498    
B:D        1      2       2   0.0008  0.978978    
C:D        1     18      18   0.0089  0.928641    
Residuals  5  10187    2037                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Modelo com efeitos significativos apenas
m1 <- update(m0, . ~ A * D, data = dados)
anova(m1)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
A          1  41311   41311  23.767 0.0003816 ***
D          1 374850  374850 215.661 4.951e-09 ***
A:D        1  94403   94403  54.312 8.621e-06 ***
Residuals 12  20858    1738                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m1)
## Quantis normais
qqnorm(res); qqline(res)

Para ilustrar o uso da meia-fração, suponha que decidamos usar o planejamento \(2^{4-1}\), com a relação de definição \(I = ABCD\), para investigar os quatro fatores deste experimento. Esse planejamento será construído como o planejamento básico \(2^3\) nos fatores A, B, C, e então estabelecendo os níveis do quarto fator \(D = ABC\), pois, através da relação de definição, temos que \(D \cdot I = D \cdot ABCD = ABC\).

O primeiro passo é então criar um planejamento \(2^3\) completo, e incluir a coluna \(D\) como sendo \(D = ABC\)

dados2 <- expand.grid(A = c(-1, 1), B = c(-1, 1), C = c(-1, 1))
dados2$D <- with(dados2, A*B*C)
row.names(dados2) <- apply(dados2, 1,
                          function(i) paste(letters[1:4][i==1], collapse = ""))
row.names(dados2)[1] <- "(1)"
dados2
      A  B  C  D
(1)  -1 -1 -1 -1
ad    1 -1 -1  1
bd   -1  1 -1  1
ab    1  1 -1 -1
cd   -1 -1  1  1
ac    1 -1  1 -1
bc   -1  1  1 -1
abcd  1  1  1  1

Outra forma de montar esse experimento é a seguinte. Se a relação de definição é \(I = ABCD\), então podemos usar os sinais dessa interação para separar as frações do experimento.

dados$ABCD <- with(dados, A*B*C*D)
dados[order(dados$ABCD, decreasing = TRUE), ]
      A  B  C  D    y ABCD
(1)  -1 -1 -1 -1  550    1
ab    1  1 -1 -1  650    1
ac    1 -1  1 -1  642    1
bc   -1  1  1 -1  601    1
ad    1 -1 -1  1  749    1
bd   -1  1 -1  1 1052    1
cd   -1 -1  1  1 1075    1
abcd  1  1  1  1  729    1
a     1 -1 -1 -1  669   -1
b    -1  1 -1 -1  604   -1
c    -1 -1  1 -1  633   -1
abc   1  1  1 -1  635   -1
d    -1 -1 -1  1 1037   -1
abd   1  1 -1  1  868   -1
acd   1 -1  1  1  860   -1
bcd  -1  1  1  1 1063   -1

As combinações de tratamento com sinal positivo em ABCD serão da fração principal, e os demais serão da fração alternada. Dessa forma, se ficarmos apenas com o sinal positivo de ABCD

dados <- subset(dados, ABCD == 1)
dados
      A  B  C  D    y ABCD
(1)  -1 -1 -1 -1  550    1
ab    1  1 -1 -1  650    1
ac    1 -1  1 -1  642    1
bc   -1  1  1 -1  601    1
ad    1 -1 -1  1  749    1
bd   -1  1 -1  1 1052    1
cd   -1 -1  1  1 1075    1
abcd  1  1  1  1  729    1

Ficamos essencialmente com o mesmo planejamento fracionário. Note que dessa forma, automaticamente temos \(D = ABC\)

dados$D == (dados$A*dados$B*dados$C)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Agora sabemos que as combinações de tratamento que serão corridas na meia-fração serão

row.names(dados)
[1] "(1)"  "ab"   "ac"   "bc"   "ad"   "bd"   "cd"   "abcd"

O próximo passo então é criar a tabela de sinais

(tab <- model.matrix(~ A*B*C*D, data = dados))
     (Intercept)  A  B  C  D A:B A:C B:C A:D B:D C:D A:B:C A:B:D A:C:D
(1)            1 -1 -1 -1 -1   1   1   1   1   1   1    -1    -1    -1
ab             1  1  1 -1 -1   1  -1  -1  -1  -1   1    -1    -1     1
ac             1  1 -1  1 -1  -1   1  -1  -1   1  -1    -1     1    -1
bc             1 -1  1  1 -1  -1  -1   1   1  -1  -1    -1     1     1
ad             1  1 -1 -1  1  -1  -1   1   1  -1  -1     1    -1    -1
bd             1 -1  1 -1  1  -1   1  -1  -1   1  -1     1    -1     1
cd             1 -1 -1  1  1   1  -1  -1  -1  -1   1     1     1    -1
abcd           1  1  1  1  1   1   1   1   1   1   1     1     1     1
     B:C:D A:B:C:D
(1)     -1       1
ab       1       1
ac       1       1
bc      -1       1
ad       1       1
bd      -1       1
cd      -1       1
abcd     1       1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

Note que os tratamentos possuem pares combinados (aliases), pois muitas colunas são idênticas. Isso é natural, uma vez que temos um efeito confundido por definição, \(D = ABC\), outros também devem estar. Para saber quem são os pares combinados desse experimento, basta multiplicar cada fator pela relação de definição, por exemplo,

\[ \begin{align} A \cdot I = A \cdot ABCD = A^2BCD = BCD \\ B \cdot I = B \cdot ABCD = AB^2CD = ACD \\ C \cdot I = C \cdot ABCD = ABC^2D = ABD \\ D \cdot I = D \cdot ABCD = ABCD^2 = ABC \\ \end{align} \]

Portanto, notamos que todos os efeitos principais estão confudidos com todas as interações de terceira ordem. Portanto, ao estimar um efeito principal, estamos também estimado a respectiva interação de terceira ordem. Como acreditamos que os efeitos das interações de ordem alta são desprezíveis, então podemos supor que o efeito estimado é de fato do efeto principal.

Vemos também que as interaçoes de segunda ordem estão associadas entre si, pois,

\[ \begin{align} AB \cdot I = AB \cdot ABCD = A^2B^2CD = CD \\ AC \cdot I = AC \cdot ABCD = A^2BC^2D = BD \\ AD \cdot I = AD \cdot ABCD = A^2BCD^2 = BC \\ \end{align} \]

Com isso, não faz sentido estimar o efeito para todas as interações de segunda ordem, assim como estimar os efeitos de terceira ordem. Estimar os efeitos principais (que é de interesse), é virtualmente igual a estimar os efeitos de terceira ordem. Estimar os efeitos de segunda ordem AB, AC, e AD, é o mesmo que estimar CD, BD, e BC, respectivamente.

Portanto, podemos identificar os fatores únicos na tabela de sinais

(tab <- unique(tab, MARGIN = 2))
     (Intercept)  A  B  C  D A:B A:C B:C
(1)            1 -1 -1 -1 -1   1   1   1
ab             1  1  1 -1 -1   1  -1  -1
ac             1  1 -1  1 -1  -1   1  -1
bc             1 -1  1  1 -1  -1  -1   1
ad             1  1 -1 -1  1  -1  -1   1
bd             1 -1  1 -1  1  -1   1  -1
cd             1 -1 -1  1  1   1  -1  -1
abcd           1  1  1  1  1   1   1   1

A partir desa tabela, podemos prosseguir a análise de forma usual, calculando os contrastes e efeitos, e verificando o tamanho dos efeitpos através do gráfico de probabilidades normais. Note apenas que o \(k\) para esse experimento fracionado é 3, e não 4 como o original, pois ele é um \(2^{4-1}\) = \(2^{3}\)

##----------------------------------------------------------------------
## Estimativa dos efeitos
y <- dados$y
## Contrastes
contr <- t(tab[, -1]) %*% y
## Efeitos = contraste/(n2^{k-1})
## NOTE que agora o k é 3, pois o experimento básico é 2^{4-1} = 2^3
k <- 3
(ef <- contr/(n * 2^(k-1)))
      [,1]
A   -127.0
B      4.0
C     11.5
D    290.5
A:B  -10.0
A:C  -25.5
B:C -197.5
##----------------------------------------------------------------------
## Gráfico de probabilidade normal dos efeitos estimados
qqaux <- qqnorm(ef, col = 2, pch = 19); qqline(ef)
text(qqaux$x, qqaux$y, rownames(qqaux$y), cex = 0.8, pos = 3)

Através destes resultados vemos que se destacam os fatores A, D,e a interação BC. Note que o par associado de BC é AD, portanto acreditamos que a interação que de fato é importante é AD. Esse resultado está de acordo com o resultado do experimento completo mostrado acima.

Continuamos a análise como antes, ajustando um modelo apenas com os efeitos que notamos serem importantes na análise visual acima.

##----------------------------------------------------------------------
## Ajuste do modelo
## Modelo com efeitos significativos apenas
## NOTE que os efeitos importantes foram A, D, e BC, mas como vimos,
## BC = AD, portanto podemos ajustar o modelo com A + D + A:D
m2 <- lm(y ~ A * D, data = dados)
anova(m2)
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
A          1  32258   32258  71.804 0.0010631 ** 
D          1 168781  168781 375.694 4.177e-05 ***
A:D        1  78013   78013 173.650 0.0001916 ***
Residuals  4   1797     449                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m2)
## Quantis normais
qqnorm(res); qqline(res)

Exerício

Considere os dados de um experimento completo \(2^5\)

url <- "http://www.leg.ufpr.br/~fernandomayer/dados/montgomery_14-41.txt"
dados <- read.table(url, header = TRUE)
dados
       A  B  C  D  E  Y
(1)   -1 -1 -1 -1 -1  7
a      1 -1 -1 -1 -1  9
b     -1  1 -1 -1 -1 34
ab     1  1 -1 -1 -1 55
c     -1 -1  1 -1 -1 16
ac     1 -1  1 -1 -1 20
bc    -1  1  1 -1 -1 40
abc    1  1  1 -1 -1 60
d     -1 -1 -1  1 -1  8
ad     1 -1 -1  1 -1 10
bd    -1  1 -1  1 -1 32
abd    1  1 -1  1 -1 50
cd    -1 -1  1  1 -1 18
acd    1 -1  1  1 -1 21
bcd   -1  1  1  1 -1 44
abcd   1  1  1  1 -1 61
e     -1 -1 -1 -1  1  8
ae     1 -1 -1 -1  1 12
be    -1  1 -1 -1  1 35
abe    1  1 -1 -1  1 52
ce    -1 -1  1 -1  1 15
ace    1 -1  1 -1  1 22
bce   -1  1  1 -1  1 45
abce   1  1  1 -1  1 65
de    -1 -1 -1  1  1  6
ade    1 -1 -1  1  1 10
bde   -1  1 -1  1  1 30
abde   1  1 -1  1  1 53
cde   -1 -1  1  1  1 15
acde   1 -1  1  1  1 20
bcde  -1  1  1  1  1 41
abcde  1  1  1  1  1 63

Suponha que somente metade das 32 corridas pode ser feita. Com isso:

  1. Escolha a metade que você pensa que deve ser corrida
  2. Determine as relações associadas para seu planejamento
  3. Estime os efeitos dos fatores
  4. Faça um gráfico de probabilidade normal para os efeitos estimados
  5. Faça a análise de variância para os fatores identificados como importantes no item anterior
  6. Analise os resíduos do modelo
  7. Forneça uma interpretação prática dos resultados
## a)
## Escolhendo a relação de definição
## I = ABCDE
## então usamos essa definição como o separador das meia frações. Cria a
## coluna ABCDE
dados$ABCDE <- with(dados, A*B*C*D*E)
## ordena pelos sinais dessa coluna
dados[order(dados$ABCDE, decreasing = TRUE), ]
       A  B  C  D  E  Y ABCDE
a      1 -1 -1 -1 -1  9     1
b     -1  1 -1 -1 -1 34     1
c     -1 -1  1 -1 -1 16     1
abc    1  1  1 -1 -1 60     1
d     -1 -1 -1  1 -1  8     1
abd    1  1 -1  1 -1 50     1
acd    1 -1  1  1 -1 21     1
bcd   -1  1  1  1 -1 44     1
e     -1 -1 -1 -1  1  8     1
abe    1  1 -1 -1  1 52     1
ace    1 -1  1 -1  1 22     1
bce   -1  1  1 -1  1 45     1
ade    1 -1 -1  1  1 10     1
bde   -1  1 -1  1  1 30     1
cde   -1 -1  1  1  1 15     1
abcde  1  1  1  1  1 63     1
(1)   -1 -1 -1 -1 -1  7    -1
ab     1  1 -1 -1 -1 55    -1
ac     1 -1  1 -1 -1 20    -1
bc    -1  1  1 -1 -1 40    -1
ad     1 -1 -1  1 -1 10    -1
bd    -1  1 -1  1 -1 32    -1
cd    -1 -1  1  1 -1 18    -1
abcd   1  1  1  1 -1 61    -1
ae     1 -1 -1 -1  1 12    -1
be    -1  1 -1 -1  1 35    -1
ce    -1 -1  1 -1  1 15    -1
abce   1  1  1 -1  1 65    -1
de    -1 -1 -1  1  1  6    -1
abde   1  1 -1  1  1 53    -1
acde   1 -1  1  1  1 20    -1
bcde  -1  1  1  1  1 41    -1
## As combações com sinal positivo de ABCDE podem ser a fração
## principal, portanto, fazemos a seleção desses elementos
dados <- subset(dados, ABCDE == 1)
## Portanto, as combinaçoes seriam
row.names(dados)
 [1] "a"     "b"     "c"     "abc"   "d"     "abd"   "acd"   "bcd"  
 [9] "e"     "abe"   "ace"   "bce"   "ade"   "bde"   "cde"   "abcde"
## b)
## As relaçãoes associadas devem ser encontradas manualmente, ou usando
## a função abaixo
aliases <- function(ef, cd){
    ali <- function(ef, cd){
        ef <- gsub("\\W", "", ef)
        ef <- unlist(strsplit(ef, ""))
        cd <- gsub("\\W", "", cd)
        cd <- unlist(strsplit(cd, ""))
        x <- sort(c(ef, cd))
        tb <- table(x)
        paste(names(tb[tb==1]), collapse="")
    }
    sapply(ef, simplify=FALSE,
           function(e){
               x <- sapply(cd, function(i){
                   x <- ali(e, i)
               })
               x[order(nchar(x))]
           })
}

## Cria a tabela de sinais
(tab <- model.matrix(~ A*B*C*D*E, data = dados))
      (Intercept)  A  B  C  D  E A:B A:C B:C A:D B:D C:D A:E B:E C:E D:E
a               1  1 -1 -1 -1 -1  -1  -1   1  -1   1   1  -1   1   1   1
b               1 -1  1 -1 -1 -1  -1   1  -1   1  -1   1   1  -1   1   1
c               1 -1 -1  1 -1 -1   1  -1  -1   1   1  -1   1   1  -1   1
abc             1  1  1  1 -1 -1   1   1   1  -1  -1  -1  -1  -1  -1   1
d               1 -1 -1 -1  1 -1   1   1   1  -1  -1  -1   1   1   1  -1
abd             1  1  1 -1  1 -1   1  -1  -1   1   1  -1  -1  -1   1  -1
acd             1  1 -1  1  1 -1  -1   1  -1   1  -1   1  -1   1  -1  -1
bcd             1 -1  1  1  1 -1  -1  -1   1  -1   1   1   1  -1  -1  -1
e               1 -1 -1 -1 -1  1   1   1   1   1   1   1  -1  -1  -1  -1
abe             1  1  1 -1 -1  1   1  -1  -1  -1  -1   1   1   1  -1  -1
ace             1  1 -1  1 -1  1  -1   1  -1  -1   1  -1   1  -1   1  -1
bce             1 -1  1  1 -1  1  -1  -1   1   1  -1  -1  -1   1   1  -1
ade             1  1 -1 -1  1  1  -1  -1   1   1  -1  -1   1  -1  -1   1
bde             1 -1  1 -1  1  1  -1   1  -1  -1   1  -1  -1   1  -1   1
cde             1 -1 -1  1  1  1   1  -1  -1  -1  -1   1  -1  -1   1   1
abcde           1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1
      A:B:C A:B:D A:C:D B:C:D A:B:E A:C:E B:C:E A:D:E B:D:E C:D:E A:B:C:D
a         1     1     1    -1     1     1    -1     1    -1    -1      -1
b         1     1    -1     1     1    -1     1    -1     1    -1      -1
c         1    -1     1     1    -1     1     1    -1    -1     1      -1
abc       1    -1    -1    -1    -1    -1    -1     1     1     1      -1
d        -1     1     1     1    -1    -1    -1     1     1     1      -1
abd      -1     1    -1    -1    -1     1     1    -1    -1     1      -1
acd      -1    -1     1    -1     1    -1     1    -1     1    -1      -1
bcd      -1    -1    -1     1     1     1    -1     1    -1    -1      -1
e        -1    -1    -1    -1     1     1     1     1     1     1       1
abe      -1    -1     1     1     1    -1    -1    -1    -1     1       1
ace      -1     1    -1     1    -1     1    -1    -1     1    -1       1
bce      -1     1     1    -1    -1    -1     1     1    -1    -1       1
ade       1    -1    -1     1    -1    -1     1     1    -1    -1       1
bde       1    -1     1    -1    -1     1    -1    -1     1    -1       1
cde       1     1    -1    -1     1    -1    -1    -1    -1     1       1
abcde     1     1     1     1     1     1     1     1     1     1       1
      A:B:C:E A:B:D:E A:C:D:E B:C:D:E A:B:C:D:E
a          -1      -1      -1       1         1
b          -1      -1       1      -1         1
c          -1       1      -1      -1         1
abc        -1       1       1       1         1
d           1      -1      -1      -1         1
abd         1      -1       1       1         1
acd         1       1      -1       1         1
bcd         1       1       1      -1         1
e          -1      -1      -1      -1         1
abe        -1      -1       1       1         1
ace        -1       1      -1       1         1
bce        -1       1       1      -1         1
ade         1      -1      -1       1         1
bde         1      -1       1      -1         1
cde         1       1      -1      -1         1
abcde       1       1       1       1         1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22
[24] 23 24 25 26 27 28 29 30 31
## Aplicando a função, encontrando os pares associados
aliases(colnames(tab)[-c(1, 32)], "ABCDE")
$A
 ABCDE 
"BCDE" 

$B
 ABCDE 
"ACDE" 

$C
 ABCDE 
"ABDE" 

$D
 ABCDE 
"ABCE" 

$E
 ABCDE 
"ABCD" 

$`A:B`
ABCDE 
"CDE" 

$`A:C`
ABCDE 
"BDE" 

$`B:C`
ABCDE 
"ADE" 

$`A:D`
ABCDE 
"BCE" 

$`B:D`
ABCDE 
"ACE" 

$`C:D`
ABCDE 
"ABE" 

$`A:E`
ABCDE 
"BCD" 

$`B:E`
ABCDE 
"ACD" 

$`C:E`
ABCDE 
"ABD" 

$`D:E`
ABCDE 
"ABC" 

$`A:B:C`
ABCDE 
 "DE" 

$`A:B:D`
ABCDE 
 "CE" 

$`A:C:D`
ABCDE 
 "BE" 

$`B:C:D`
ABCDE 
 "AE" 

$`A:B:E`
ABCDE 
 "CD" 

$`A:C:E`
ABCDE 
 "BD" 

$`B:C:E`
ABCDE 
 "AD" 

$`A:D:E`
ABCDE 
 "BC" 

$`B:D:E`
ABCDE 
 "AC" 

$`C:D:E`
ABCDE 
 "AB" 

$`A:B:C:D`
ABCDE 
  "E" 

$`A:B:C:E`
ABCDE 
  "D" 

$`A:B:D:E`
ABCDE 
  "C" 

$`A:C:D:E`
ABCDE 
  "B" 

$`B:C:D:E`
ABCDE 
  "A" 
## c)
## Como vimos que existem colunas confundidas, vamos estimar os efeitos
## apenas daquelas únicas
(tab <- unique(tab, MARGIN = 2))
      (Intercept)  A  B  C  D  E A:B A:C B:C A:D B:D C:D A:E B:E C:E D:E
a               1  1 -1 -1 -1 -1  -1  -1   1  -1   1   1  -1   1   1   1
b               1 -1  1 -1 -1 -1  -1   1  -1   1  -1   1   1  -1   1   1
c               1 -1 -1  1 -1 -1   1  -1  -1   1   1  -1   1   1  -1   1
abc             1  1  1  1 -1 -1   1   1   1  -1  -1  -1  -1  -1  -1   1
d               1 -1 -1 -1  1 -1   1   1   1  -1  -1  -1   1   1   1  -1
abd             1  1  1 -1  1 -1   1  -1  -1   1   1  -1  -1  -1   1  -1
acd             1  1 -1  1  1 -1  -1   1  -1   1  -1   1  -1   1  -1  -1
bcd             1 -1  1  1  1 -1  -1  -1   1  -1   1   1   1  -1  -1  -1
e               1 -1 -1 -1 -1  1   1   1   1   1   1   1  -1  -1  -1  -1
abe             1  1  1 -1 -1  1   1  -1  -1  -1  -1   1   1   1  -1  -1
ace             1  1 -1  1 -1  1  -1   1  -1  -1   1  -1   1  -1   1  -1
bce             1 -1  1  1 -1  1  -1  -1   1   1  -1  -1  -1   1   1  -1
ade             1  1 -1 -1  1  1  -1  -1   1   1  -1  -1   1  -1  -1   1
bde             1 -1  1 -1  1  1  -1   1  -1  -1   1  -1  -1   1  -1   1
cde             1 -1 -1  1  1  1   1  -1  -1  -1  -1   1  -1  -1   1   1
abcde           1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1
## Contrastes
contr <- t(tab[, -1]) %*% dados$Y
## Efeitos = contraste/(n2^{k-1})
## NOTE que agora o k é 4, pois o experimento básico é 2^{5-1} = 2^4
k <- 4
n <- 1
(ef <- contr/(n * 2^(k-1)))
      [,1]
A   10.875
B   33.625
C   10.625
D   -0.625
E    0.375
A:B  7.125
A:C  0.625
B:C  0.875
A:D  0.875
B:D -0.375
C:D  0.625
A:E  1.375
B:E  0.125
C:E  0.625
D:E -1.625
## d)
## Gráfico de probabilidade normal dos efeitos estimados
qqaux <- qqnorm(ef, col = 2, pch = 19); qqline(ef)
text(qqaux$x, qqaux$y, rownames(qqaux$y), cex = 0.8, pos = 3)

## e)
## Os efeitos iportantes foram A, B, C, AB, e DE. Note que
## DE . ADCDE = ABC
## entao podemos especificar o modelo com
m0 <- lm(Y ~ A + B + C + A:B + A:B:C, data = dados)
anova(m0)
Analysis of Variance Table

Response: Y
          Df Sum Sq Mean Sq  F value    Pr(>F)    
A          1  473.1   473.1  223.935 3.577e-08 ***
B          1 4522.6  4522.6 2140.858 5.357e-13 ***
C          1  451.6   451.6  213.757 4.472e-08 ***
A:B        1  203.1   203.1   96.124 1.905e-06 ***
A:B:C      1   10.6    10.6    5.000   0.04933 *  
Residuals 10   21.1     2.1                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## f)
res <- residuals(m0)
qqnorm(res); qqline(res)