Experimentos fatoriais \(2^k\)
Análise de um experimento \(2^3\)
Exemplo 14-2, Montgomery (EAPE), pg. 375.
Um engenheiro mecânico está estudando a rugosidade da superfície de uma peça usinada. Três fatores são de interesse: taxa de alimentação (A), profundidade do corte (B), e ângulo da ferramenta (C). Todos os três fatores tem 2 níveis, e um planejamento fatorial foi feito com 2 réplicas.
A representação geométrica é
nom A B C
1 1 - - -
2 a + - -
3 b - + -
4 ab + + -
5 c - - +
6 ac + - +
7 bc - + +
8 abc + + + bc ------- abc
.| .|
[+] c__|______ac |
| | | |
C | b -----|- ab [+]
| . | . B
[-] (1)________a [-]
[-] A [+]## Importação
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-2.txt"
dados <- read.table(url, header = TRUE)
str(dados)'data.frame': 16 obs. of 4 variables:
$ Alim : int 20 20 30 30 20 20 30 30 20 20 ...
$ Prof : num 0.025 0.025 0.025 0.025 0.04 0.04 0.04 0.04 0.025 0.025 ...
$ Angulo : int 15 15 15 15 15 15 15 15 25 25 ...
$ Rugosidade: int 9 7 10 12 9 11 12 15 11 10 ...## Visualizacao
par(mfrow=c(2,2))
with(dados, interaction.plot(Alim, Prof, Rugosidade))
with(dados, interaction.plot(Alim, Angulo, Rugosidade))
with(dados, interaction.plot(Prof, Angulo, Rugosidade))
par(mfrow=c(1,1))
wireframe(Rugosidade ~ Alim * Prof, data = dados)
wireframe(Rugosidade ~ Alim * Angulo, data = dados)
wireframe(Rugosidade ~ Prof * Angulo, data = dados)
xyplot(Rugosidade ~ Alim | Prof, groups = Angulo, data = dados,
type = c("p", "a"), auto.key=TRUE)
xyplot(Rugosidade ~ Alim, groups = interaction(Prof, Angulo, sep = ":"),
data = dados, type = c("p", "a"), auto.key = TRUE)
##----------------------------------------------------------------------
## Analise usando as funcoes do R
## Transforma as variáveis em fator
dados <- transform(dados, Alim = factor(Alim),
Prof = factor(Prof),
Angulo = factor(Angulo))
## Ajusta o modelo, especificando o conraste soma zero
m0 <- lm(Rugosidade ~ Alim * Prof * Angulo, data = dados,
contrasts = list(Alim = "contr.sum",
Prof = "contr.sum",
Angulo = "contr.sum"))
anova(m0)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim 1 45.562 45.562 18.6923 0.002534 **
Prof 1 10.562 10.562 4.3333 0.070931 .
Angulo 1 3.062 3.062 1.2564 0.294849
Alim:Prof 1 7.563 7.563 3.1026 0.116197
Alim:Angulo 1 0.063 0.063 0.0256 0.876749
Prof:Angulo 1 1.563 1.563 0.6410 0.446463
Alim:Prof:Angulo 1 5.062 5.062 2.0769 0.187512
Residuals 8 19.500 2.438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)
Call:
lm(formula = Rugosidade ~ Alim * Prof * Angulo, data = dados,
contrasts = list(Alim = "contr.sum", Prof = "contr.sum",
Angulo = "contr.sum"))
Residuals:
Min 1Q Median 3Q Max
-1.5 -1.0 0.0 1.0 1.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.0625 0.3903 28.343 2.6e-09 ***
Alim1 -1.6875 0.3903 -4.323 0.00253 **
Prof1 -0.8125 0.3903 -2.082 0.07093 .
Angulo1 -0.4375 0.3903 -1.121 0.29485
Alim1:Prof1 0.6875 0.3903 1.761 0.11620
Alim1:Angulo1 0.0625 0.3903 0.160 0.87675
Prof1:Angulo1 -0.3125 0.3903 -0.801 0.44646
Alim1:Prof1:Angulo1 -0.5625 0.3903 -1.441 0.18751
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.561 on 8 degrees of freedom
Multiple R-squared: 0.7902, Adjusted R-squared: 0.6066
F-statistic: 4.304 on 7 and 8 DF, p-value: 0.02881Repare que os sinais dos coeficientes estimados pela função lm() são contrários aos sinais que calculamos a mão. Isso ocorre porque no contraste soma zero do R, é assumido que o último nível é igual ao negativo do primeiro, ou seja, \[
\alpha_2 = -\alpha_1
\]
Como, por padrão, assumimos calcular os efeitos usando a diferença de médias sempre no formato efeito = nivel_alto - nivel_baixo então precisamos assumir que \[
\alpha_1 = -\alpha_2
\]
Do ponto de vista prático, o inversão dos sinais não altera em nada a interpretação e os resultados do modelo, mas podemos fazer a saída do R ficar compatível com a nossa notação. Para isso, basta utilizarmos a codificação numérica -1 e 1 para os níveis baixo e alto dos fatores, ao invés de utilizarmos a codificação padrão do R.
## Funcao para converter os niveis de um fator na codificacao -1 (baixo)
## e 1 (alto). Valido apenas para fatores com 2 niveis.
convert <- function(x){
if(!is.factor(x)){
x <- as.factor(x)
}
if(!identical(nlevels(x), 2L)){
stop("O fator deve ter exatamente 2 niveis")
}
x <- as.numeric(x)
cod <- (x - mean(x))/0.5
return(cod)
}## Convertendo os níveis dos fatores para -1 e 1
dados <- transform(dados, Alim.n = convert(Alim),
Prof.n = convert(Prof),
Angulo.n = convert(Angulo))
## Ajusta o modelo com a nova codificação
m1 <- lm(Rugosidade ~ Alim.n * Prof.n * Angulo.n, data = dados)
anova(m1)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim.n 1 45.562 45.562 18.6923 0.002534 **
Prof.n 1 10.562 10.562 4.3333 0.070931 .
Angulo.n 1 3.062 3.062 1.2564 0.294849
Alim.n:Prof.n 1 7.563 7.563 3.1026 0.116197
Alim.n:Angulo.n 1 0.063 0.063 0.0256 0.876749
Prof.n:Angulo.n 1 1.563 1.563 0.6410 0.446463
Alim.n:Prof.n:Angulo.n 1 5.062 5.062 2.0769 0.187512
Residuals 8 19.500 2.438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)
Call:
lm(formula = Rugosidade ~ Alim.n * Prof.n * Angulo.n, data = dados)
Residuals:
Min 1Q Median 3Q Max
-1.5 -1.0 0.0 1.0 1.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.0625 0.3903 28.343 2.6e-09 ***
Alim.n 1.6875 0.3903 4.323 0.00253 **
Prof.n 0.8125 0.3903 2.082 0.07093 .
Angulo.n 0.4375 0.3903 1.121 0.29485
Alim.n:Prof.n 0.6875 0.3903 1.761 0.11620
Alim.n:Angulo.n 0.0625 0.3903 0.160 0.87675
Prof.n:Angulo.n -0.3125 0.3903 -0.801 0.44646
Alim.n:Prof.n:Angulo.n 0.5625 0.3903 1.441 0.18751
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.561 on 8 degrees of freedom
Multiple R-squared: 0.7902, Adjusted R-squared: 0.6066
F-statistic: 4.304 on 7 and 8 DF, p-value: 0.02881Suposições do modelo
Como todo modelo de ANOVA:
- \(\epsilon_{ijkl} \sim \text{N}(0, \sigma^2)\)
- \(\epsilon_{ijkl}\) são independentes
- Homogeneidade de variâncias
Avaliação do modelo
Quando ajustamos um modelo a um conjunto de dados, podemos ter basicamente dois objetivos:
- Verificar quais variáveis influenciam na média da variável resposta \(Y\)
- Fazer predições com base nos parâmetros estimados
- Outros: estimar parâmetros e fazer inferência
Para verificar as variáves que infuenciam na rugosidade, vemos a ANOVA do modelo,
anova(m1)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim.n 1 45.562 45.562 18.6923 0.002534 **
Prof.n 1 10.562 10.562 4.3333 0.070931 .
Angulo.n 1 3.062 3.062 1.2564 0.294849
Alim.n:Prof.n 1 7.563 7.563 3.1026 0.116197
Alim.n:Angulo.n 1 0.063 0.063 0.0256 0.876749
Prof.n:Angulo.n 1 1.563 1.563 0.6410 0.446463
Alim.n:Prof.n:Angulo.n 1 5.062 5.062 2.0769 0.187512
Residuals 8 19.500 2.438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se o nosso objetivo fosse simplesmente construir um modelo para avaliar quais efeitos e interações influemciam no processo, então a análise poderia parar aqui, pois já temos esse resultado. Então deveria se proceder para a análise dos resíduos para verificar as suposições do modelo.
Caso o nosso objetivo seja fazer uma predição para a variável resposta, então devemos ajustar o modelo mais parcimonioso possível. Um modelo parcimonioso é aquele que consegue explicar a maior parte da variabilidade dos dados, com o menor número de parâmetros possível.
Asumindo uma avaliação não muito “fechada” (ou restrita) em um ponto de corte de p-valor especifico, notamos que Alim e Prof são fatores influentes. Como os dois são influentes, e a interação entre eles é “marginalmente” significativa, vamos manter estes dois fatores principais e a interação para construir um novo modelo.
m2 <- update(m1, . ~ Alim.n * Prof.n)Como, por construção, o modelo é ortogonal, então as somas de quadrados e os coeficientes estimados também serão ortogonais. Dessa forma, a retirada de termos do modelo não interfere nos coeficientes estimados. No caso da ANOVA, as somas de quadrados permanecem as mesmas, mas a \(SQRes\) e os \(GL\) associados aumentam, pois as \(SQ\) dos fatores retirados são incorporados aos resíduos. Veja:
anova(m1)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim.n 1 45.562 45.562 18.6923 0.002534 **
Prof.n 1 10.562 10.562 4.3333 0.070931 .
Angulo.n 1 3.062 3.062 1.2564 0.294849
Alim.n:Prof.n 1 7.563 7.563 3.1026 0.116197
Alim.n:Angulo.n 1 0.063 0.063 0.0256 0.876749
Prof.n:Angulo.n 1 1.563 1.563 0.6410 0.446463
Alim.n:Prof.n:Angulo.n 1 5.062 5.062 2.0769 0.187512
Residuals 8 19.500 2.438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim.n 1 45.562 45.562 18.6923 0.0009901 ***
Prof.n 1 10.562 10.562 4.3333 0.0594472 .
Alim.n:Prof.n 1 7.563 7.563 3.1026 0.1036022
Residuals 12 29.250 2.437
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Porque se alteram os p-valores?
## Exemplo para Alim.n
## Usando o modelo completo
pf(anova(m1)$"F value"[1], df1 = anova(m1)$"Df"[1],
df2 = anova(m1)$"Df"[8], lower.tail = FALSE)[1] 0.002534218## Exemplo para o modelo sob avaliação
pf(anova(m2)$"F value"[1], df1 = anova(m2)$"Df"[1],
df2 = anova(m2)$"Df"[4], lower.tail = FALSE)[1] 0.0009900542Porque mudam os graus de liberdade do numerador e do denominador, altera a distribuição.
Para comparar se o modelo reduzido pode de fato ser utilizado no lugar do modelo completo, podemos fazer o teste da razão de verossimilhança (TRV) entre os dois modelos encaixados.
anova(m2, m1)Analysis of Variance Table
Model 1: Rugosidade ~ Alim.n + Prof.n + Alim.n:Prof.n
Model 2: Rugosidade ~ Alim.n * Prof.n * Angulo.n
Res.Df RSS Df Sum of Sq F Pr(>F)
1 12 29.25
2 8 19.50 4 9.75 1 0.4609Como a diferença das somas de quadrados residuais entre os dois modelos não é significativa, então concluímos que o modelo com menos parâmetros pode ser utilizado no lugar do modelo com mais parâmetros - este é o princípio da parcimônia.
Podemos ainda estar em dúvida sobre a inclusão do termo de interação, que é “marginalmente” significativo. Para avaliar se de fato esse termo é importante, podemos ajustar um novo modelo sem esse termo
m3 <- update(m2, . ~ Alim.n + Prof.n)
anova(m3)Analysis of Variance Table
Response: Rugosidade
Df Sum Sq Mean Sq F value Pr(>F)
Alim.n 1 45.562 45.562 16.0900 0.00148 **
Prof.n 1 10.562 10.562 3.7301 0.07554 .
Residuals 13 36.812 2.832
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Assim como antes, podemos usar o TRV para comparar os modelos encaixados
anova(m3, m2, m1)Analysis of Variance Table
Model 1: Rugosidade ~ Alim.n + Prof.n
Model 2: Rugosidade ~ Alim.n + Prof.n + Alim.n:Prof.n
Model 3: Rugosidade ~ Alim.n * Prof.n * Angulo.n
Res.Df RSS Df Sum of Sq F Pr(>F)
1 13 36.812
2 12 29.250 1 7.5625 3.1026 0.1162
3 8 19.500 4 9.7500 1.0000 0.4609Com esse resultado, podemos então concluir que a interação não é importante, e vamos dizer que o último modelo é o modelo escolhido.
Análise dos resíduos
Para verificar as suposições feitas para o modelo, procedemos com a análise dos resíduos do modelo final
## Residuos de modelos
res <- residuals(m3)
## Quantis normais
qqnorm(res); qqline(res)
## Residuos vs fatores
library(lattice)
xyplot(res ~ Alim, data = dados,
panel = function(...){
panel.xyplot(...)
panel.abline(h = 0, col = 2, lty = 2)
})
xyplot(res ~ Prof, data = dados,
panel = function(...){
panel.xyplot(...)
panel.abline(h = 0, col = 2, lty = 2)
})
Predição
A partir dessa avaliação, podemos fazer a predição usando os coeficientes estimados do modelo
summary(m3)
Call:
lm(formula = Rugosidade ~ Alim.n + Prof.n, data = dados)
Residuals:
Min 1Q Median 3Q Max
-2.188 -1.562 0.250 1.156 2.438
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.0625 0.4207 26.296 1.17e-12 ***
Alim.n 1.6875 0.4207 4.011 0.00148 **
Prof.n 0.8125 0.4207 1.931 0.07554 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.683 on 13 degrees of freedom
Multiple R-squared: 0.6039, Adjusted R-squared: 0.543
F-statistic: 9.91 on 2 and 13 DF, p-value: 0.002431## Predição para as combinações únicas dos fatores
pred <- expand.grid(Alim.n = unique(dados$Alim.n),
Prof.n = unique(dados$Prof.n))
pred$Rugosidade <- predict(m3, newdata = pred)
pred Alim.n Prof.n Rugosidade
1 -1 -1 8.5625
2 1 -1 11.9375
3 -1 1 10.1875
4 1 1 13.5625## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(Alim.n = seq(-1, 1, length.out = 20),
Prof.n = seq(-1, 1 ,length.out = 20))
pred$Rugosidade <- predict(m3, newdata = pred)
## Vários formas de visualizar
wireframe(Rugosidade ~ Alim.n + Prof.n, data = pred)
wireframe(Rugosidade ~ Alim.n + Prof.n, data = pred, drape = TRUE)
levelplot(Rugosidade ~ Alim.n + Prof.n, data = pred)
levelplot(Rugosidade ~ Alim.n + Prof.n, data = pred, cuts = 90,
col.regions = heat.colors)
Fatoriais \(2^k\) para \(k > 3\)
Á medida que o número de fatores cresce em um experimento fatorial, o número de efeitos que podem ser estimados também cresce. Por exemplo,
- \(2^4\) possui 4 efeitos principais, 6 interações de segunda ordem, 4 interações de terceira ordem, e 1 interação de quarta ordem
- \(2^6\) possui 6 efeitos principais, 15 interações de segunda ordem, 20 interações de terceira ordem, 15 interações de quarta ordem, 6 interações de quinta ordem, e 1 interação de sexta ordem
Para um experimento \(2^4\), a representação geométrica fica
nom A B C D
1 (1) - - - -
2 a + - - -
3 b - + - -
4 ab + + - -
5 c - - + -
6 ac + - + -
7 bc - + + -
8 abc + + + -
9 d - - - +
10 ad + - - +
11 bd - + - +
12 abd + + - +
13 cd - - + +
14 acd + - + +
15 bcd - + + +
16 abcd + + + + D [-] | D [+]
|
bc ------- abc | bcd ------- abcd
.| .| | .| .|
[+] c__|______ac | | [+] cd__|_____acd |
| | | | | | | | |
C | b -----|- ab [+] | C | bd -----|- abd [+]
| . | . B | | . | . B
[-] (1)________a [-] | [-] d _______ ad [-]
[-] A [+] | [-] A [+]Em muitas situações, o princípio da esparsidade dos efeitos se aplica, ou seja, o sistema geralmente é dominado pelos efeitos principais e pelas interações de ordens baixas. As interações de terceira ordem ou superiores geralmente são negligenciadas.
Como consequência, quando o número de fatores for moderadamente grande, como \(k \geq 4\) ou \(5\), uma prática comum é correr somente uma réplica do planejamento, e então combinar as interações de ordem alta como uma estimativa do erro. No entanto, se eventualmente alguma interação de ordem alta for significativa, então esse procedimento não é adequado.
Um método simples para verificar o tamanho dos efeitos de um fator é construindo um gráfico das estimativas dos efeitos em uma escala de probabilidade normal. Por construção, a distribuição dos efeitos estimados possui média zero e variância \(\sigma^2\). Portanto, os efeitos que forem desprezíveis estarão em cima de uma linha reta nesse gráfico, enquanto que efeitos significativos não terão média zero e estarão mais afastados dessa linha.
Exemplo de um experimento \(2^4\)
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##----------------------------------------------------------------------
## Importação
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-5.txt"
dados <- read.table(url, header = TRUE)
## Definições do experimento
k <- 4
n <- 1
##----------------------------------------------------------------------
## Visualização
xyplot(y ~ A, groups = interaction(B, C, D, sep = ":"), data = dados,
type=c("p", "a"), auto.key = TRUE)
xyplot(y ~ A | D, groups = interaction(B, C, sep = ":"), data = dados,
type = c("p", "a"), auto.key = TRUE)
xyplot(y ~ A | B + C, groups = D, data = dados,
type = c("p", "a"), auto.key = TRUE)
##----------------------------------------------------------------------
## Montando a tabela de sinais
tab <- model.matrix(~A * B * C * D, data = dados)
tab (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
9 1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1
10 1 1 -1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1 1
11 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1
12 1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1 -1 -1
13 1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 -1 -1
14 1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 -1
15 1 -1 1 1 1 -1 -1 1 -1 1 1 -1 -1 -1 1
16 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
9 -1
10 1
11 1
12 -1
13 1
14 -1
15 -1
16 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) [,1]
A -813
B -13
C 59
D 2449
A:B -63
A:C -199
B:C -351
A:D -1229
B:D -5
C:D -17
A:B:C -125
A:B:D 33
A:C:D 45
B:C:D -203
A:B:C:D -321## 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.125Para testar quais efeitos são significativos, poderíamos pensar inicialmente em fazer uma ANOVA. No entanto, este experimento não possui repetição, e temos um parâmetro para cada observação, fazendo ocm que o modelo seja saturado. Dessa forma, não temos como avaliar a significância dos efeitos. Veja o resultado da ANOVA para o modelo completo:
anova(lm(y ~ (A * B * C * D), data = dados))Warning in anova.lm(lm(y ~ (A * B * C * D), data = dados)): ANOVA F-tests
on an essentially perfect fit are unreliableAnalysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
A 1 41311 41311
B 1 11 11
C 1 218 218
D 1 374850 374850
A:B 1 248 248
A:C 1 2475 2475
B:C 1 7700 7700
A:D 1 94403 94403
B:D 1 2 2
C:D 1 18 18
A:B:C 1 977 977
A:B:D 1 68 68
A:C:D 1 127 127
B:C:D 1 2576 2576
A:B:C:D 1 6440 6440
Residuals 0 0 Uma forma de verificar quais efeitos são importantes é através de um gráfico de probabilidade normal para a estimativa dos esfeitos. Lembre-se que, sob a hipótses nula, todos os parâmetros do modelo são iguais a zero, e pela definição do termo de erro do modelo, se os parâmetros são iguais a zero, então sobra apenas o erro que possui média 0 e variância contante \(\sigma^2\). Dessa forma, se a hipótese nula for verdadeira, esperamos que os efeitos tenham também média 0, e fiquem em cima da linha em um gráfico de probabilidade normal. Os efeitos que se afastarem muito da linha são aqueles que possuem média diferente de zero, e portanto, são aqueles que temos interesse.
##----------------------------------------------------------------------
## 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 do gráfico acima, vemos que termos de ordem maior que 2 ficaram no centro, o que mostra que interações de ordem grande podem ser consideradas como nulas. Portanto podemos ajustar agora um modelo considerando apenas as interações de segunda ordem.
##----------------------------------------------------------------------
## 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 ' ' 1Note que a variabilidade das interações não consideradas foi integrada na soma de quadrado dos resíduos, e assim podemos ter uma estimativa para \(\sigma^2\).
Com esse resultado, confirmamos a hipótese leventada no gráfico de probabilidades normais de que apenas 2 fatores e uma interação são importantes. Por isso, podemos atualizar o modelo apenas com estes termos
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## Fazendo o TRV para confirmar o modelo
anova(m1, m0)Analysis of Variance Table
Model 1: y ~ A + D + A:D
Model 2: y ~ (A + B + C + D)^2
Res.Df RSS Df Sum of Sq F Pr(>F)
1 12 20858
2 5 10187 7 10671 0.7482 0.6498Análise dos resíduos
##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m1)
## Quantis normais
qqnorm(res); qqline(res)
## Residuos vs fatores
library(lattice)
xyplot(res ~ A, data = dados,
panel = function(...){
panel.xyplot(...)
panel.abline(h = 0, col = 2, lty = 2)
})
xyplot(res ~ D, data = dados,
panel = function(...){
panel.xyplot(...)
panel.abline(h = 0, col = 2, lty = 2)
})
xyplot(res ~ interaction(A, D, sep = ":"), data = dados,
panel = function(...){
panel.xyplot(...)
panel.abline(h = 0, col = 2, lty = 2)
})
Predições
##----------------------------------------------------------------------
## Predições
## Predição para as combinações únicas dos fatores
pred <- expand.grid(A = unique(dados$A),
D = unique(dados$B))
pred$y <- predict(m1, newdata = pred)
pred A D y
1 -1 -1 597.00
2 1 -1 649.00
3 -1 1 1056.75
4 1 1 801.50Note como fica a predição para o conjunto de dados
cbind(dados, predict(m1)) A B C D y predict(m1)
1 -1 -1 -1 -1 550 597.00
2 1 -1 -1 -1 669 649.00
3 -1 1 -1 -1 604 597.00
4 1 1 -1 -1 650 649.00
5 -1 -1 1 -1 633 597.00
6 1 -1 1 -1 642 649.00
7 -1 1 1 -1 601 597.00
8 1 1 1 -1 635 649.00
9 -1 -1 -1 1 1037 1056.75
10 1 -1 -1 1 749 801.50
11 -1 1 -1 1 1052 1056.75
12 1 1 -1 1 868 801.50
13 -1 -1 1 1 1075 1056.75
14 1 -1 1 1 860 801.50
15 -1 1 1 1 1063 1056.75
16 1 1 1 1 729 801.50As médias estimadas só levam em consideração a mudança nos níveis dos fatores importantes para o processo (nesse caso A, D e A:D), independente dos níveis dos outros fatores.
## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(A = seq(-1, 1, length.out = 30),
D = seq(-1, 1 ,length.out = 30))
pred$y <- predict(m1, newdata = pred)
## Vários formas de visualizar
wireframe(y ~ A + D, data = pred, drape = TRUE)
levelplot(y ~ A + D, data = pred)
levelplot(y ~ A + D, data = pred, cuts = 90,
col.regions = heat.colors)
EXTRA: visualização dos resultados
##----------------------------------------------------------------------
## Carregando função para fazer contornos de nível
source("panel.3d.contour.R")
##----------------------------------------------------------------------
## Gráfico com contornos de nível
display.brewer.all() # palhetas de cores predefinidas
colr <- brewer.pal(11, "Spectral") # escolhe uma palheta
colr <- colorRampPalette(colr, space="rgb") # interpola para ter mais níveis
wireframe(y ~ A + D, data = pred,
scales = list(arrows = FALSE),
zlim = c(300, 1100),
col = "gray50", col.contour = 1,
panel.3d.wireframe = "panel.3d.contour",
type = c("on", "bottom"),
col.regions = colr(100), drape = TRUE, alpha.regions = 0.5,
screen = list(z = 40, x = -70))
##----------------------------------------------------------------------
## Escolhendo a posição de observação do gráfico iterativamente
require(rpanel)
wire <- function(panel){
print(
wireframe(y ~ A + D, data = pred,
scales = list(arrows = FALSE),
zlim = c(300, 1100),
col = "gray50", col.contour = 1,
panel.3d.wireframe = "panel.3d.contour",
type = c("on","bottom"),
col.regions = colr(100), drape = TRUE, alpha.regions = 0.5,
screen = list(x = panel$x, z = panel$z, y = panel$y))
)
panel
}
panel <- rp.control()
rp.slider(panel, x, -90, 90, initval=-70, showvalue=TRUE, action=wire)
rp.slider(panel, z, -90, 90, initval=40, showvalue=TRUE, action=wire)
rp.slider(panel, y, -90, 90, initval=0, showvalue=TRUE, action=wire)A partir desse ponto, o experimentador possui duas opções:
- Fazer um novo experimento, agora com apenas estes dois fatores e com um número maior de repetições: projeção
- Fazer um novo experimento com estes dois fatores e, se desejado, outros fatores que não foram considerados no primeiro experimento. Nesse caso é preciso lembrar que as podem haver interações entre estes novos fatores e aqueles que foram descartados no primeiro experimento, e isto não está sendo considerado neste caso.
Análise de um experimento fatorial \(2^6\)
##----------------------------------------------------------------------
## Dados
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/fatorial_2-6.txt"
dados <- read.table(url, header = TRUE)
str(dados)'data.frame': 64 obs. of 7 variables:
$ A: int -1 1 -1 1 -1 1 -1 1 -1 1 ...
$ B: int -1 -1 1 1 -1 -1 1 1 -1 -1 ...
$ C: int -1 -1 -1 -1 1 1 1 1 -1 -1 ...
$ D: int -1 -1 -1 -1 -1 -1 -1 -1 1 1 ...
$ E: int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
$ F: int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
$ y: int 36 35 37 34 30 32 28 38 40 32 ...## Definições
k <- 6
n <- 1
##----------------------------------------------------------------------
## Visualização
xyplot(y ~ A | B * C, groups = interaction(D, E, F), data = dados,
type = c("p", "a"))
xyplot( y ~A | B * C * D, groups = interaction(E, F), data = dados,
type = c("p", "a"))
##----------------------------------------------------------------------
## Montando a tabela de sinais
tab <- model.matrix(~A * B * C * D * E * F, data = dados)
tab (Intercept) A B C D E F A:B A:C B:C A:D B:D C:D A:E B:E C:E D:E
1 1 -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 1 1
3 1 -1 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 1 1
5 1 -1 -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 -1 1
7 1 -1 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 -1 1
9 1 -1 -1 -1 1 -1 -1 1 1 1 -1 -1 -1 1 1 1 -1
10 1 1 -1 -1 1 -1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1
11 1 -1 1 -1 1 -1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1
12 1 1 1 -1 1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1
13 1 -1 -1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1
14 1 1 -1 1 1 -1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1
15 1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1
16 1 1 1 1 1 -1 -1 1 1 1 1 1 1 -1 -1 -1 -1
17 1 -1 -1 -1 -1 1 -1 1 1 1 1 1 1 -1 -1 -1 -1
18 1 1 -1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1
19 1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1
20 1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1
21 1 -1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1
22 1 1 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1
23 1 -1 1 1 -1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1
24 1 1 1 1 -1 1 -1 1 1 1 -1 -1 -1 1 1 1 -1
25 1 -1 -1 -1 1 1 -1 1 1 1 -1 -1 -1 -1 -1 -1 1
26 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 -1 -1 1
27 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 -1 -1 1 -1 1
28 1 1 1 -1 1 1 -1 1 -1 -1 1 1 -1 1 1 -1 1
29 1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1
30 1 1 -1 1 1 1 -1 -1 1 -1 1 -1 1 1 -1 1 1
31 1 -1 1 1 1 1 -1 -1 -1 1 -1 1 1 -1 1 1 1
32 1 1 1 1 1 1 -1 1 1 1 1 1 1 1 1 1 1
33 1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 1 1 1
34 1 1 -1 -1 -1 -1 1 -1 -1 1 -1 1 1 -1 1 1 1
35 1 -1 1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 1
36 1 1 1 -1 -1 -1 1 1 -1 -1 -1 -1 1 -1 -1 1 1
37 1 -1 -1 1 -1 -1 1 1 -1 -1 1 1 -1 1 1 -1 1
38 1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 1
39 1 -1 1 1 -1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1 1
40 1 1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 1
41 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 1 -1
42 1 1 -1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 1 -1
43 1 -1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1
44 1 1 1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1
45 1 -1 -1 1 1 -1 1 1 -1 -1 -1 -1 1 1 1 -1 -1
46 1 1 -1 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 -1
47 1 -1 1 1 1 -1 1 -1 -1 1 -1 1 1 1 -1 -1 -1
48 1 1 1 1 1 -1 1 1 1 1 1 1 1 -1 -1 -1 -1
49 1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 -1 -1 -1 -1
50 1 1 -1 -1 -1 1 1 -1 -1 1 -1 1 1 1 -1 -1 -1
51 1 -1 1 -1 -1 1 1 -1 1 -1 1 -1 1 -1 1 -1 -1
52 1 1 1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 -1 -1
53 1 -1 -1 1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 1 -1
54 1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1
55 1 -1 1 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1
56 1 1 1 1 -1 1 1 1 1 1 -1 -1 -1 1 1 1 -1
57 1 -1 -1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1
58 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 1 -1 -1 1
59 1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 -1 1 -1 1
60 1 1 1 -1 1 1 1 1 -1 -1 1 1 -1 1 1 -1 1
61 1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 1 -1 -1 1 1
62 1 1 -1 1 1 1 1 -1 1 -1 1 -1 1 1 -1 1 1
63 1 -1 1 1 1 1 1 -1 -1 1 -1 1 1 -1 1 1 1
64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
A:F B:F C:F D:F E:F A:B:C A:B:D A:C:D B:C:D A:B:E A:C:E B:C:E A:D:E
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
3 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
5 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
7 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
9 1 1 1 -1 1 -1 1 1 1 -1 -1 -1 1
10 -1 1 1 -1 1 1 -1 -1 1 1 1 -1 -1
11 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 1
12 -1 -1 1 -1 1 -1 1 -1 -1 -1 1 1 -1
13 1 1 -1 -1 1 1 1 -1 -1 -1 1 1 1
14 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1
15 1 -1 -1 -1 1 -1 -1 -1 1 1 1 -1 1
16 -1 -1 -1 -1 1 1 1 1 1 -1 -1 -1 -1
17 1 1 1 1 -1 -1 -1 -1 -1 1 1 1 1
18 -1 1 1 1 -1 1 1 1 -1 -1 -1 1 -1
19 1 -1 1 1 -1 1 1 -1 1 -1 1 -1 1
20 -1 -1 1 1 -1 -1 -1 1 1 1 -1 -1 -1
21 1 1 -1 1 -1 1 -1 1 1 1 -1 -1 1
22 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1
23 1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 1
24 -1 -1 -1 1 -1 1 -1 -1 -1 1 1 1 -1
25 1 1 1 -1 -1 -1 1 1 1 1 1 1 -1
26 -1 1 1 -1 -1 1 -1 -1 1 -1 -1 1 1
27 1 -1 1 -1 -1 1 -1 1 -1 -1 1 -1 -1
28 -1 -1 1 -1 -1 -1 1 -1 -1 1 -1 -1 1
29 1 1 -1 -1 -1 1 1 -1 -1 1 -1 -1 -1
30 -1 1 -1 -1 -1 -1 -1 1 -1 -1 1 -1 1
31 1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 1 -1
32 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1
33 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
34 1 -1 -1 -1 -1 1 1 1 -1 1 1 -1 1
35 -1 1 -1 -1 -1 1 1 -1 1 1 -1 1 -1
36 1 1 -1 -1 -1 -1 -1 1 1 -1 1 1 1
37 -1 -1 1 -1 -1 1 -1 1 1 -1 1 1 -1
38 1 -1 1 -1 -1 -1 1 -1 1 1 -1 1 1
39 -1 1 1 -1 -1 -1 1 1 -1 1 1 -1 -1
40 1 1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1
41 -1 -1 -1 1 -1 -1 1 1 1 -1 -1 -1 1
42 1 -1 -1 1 -1 1 -1 -1 1 1 1 -1 -1
43 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1
44 1 1 -1 1 -1 -1 1 -1 -1 -1 1 1 -1
45 -1 -1 1 1 -1 1 1 -1 -1 -1 1 1 1
46 1 -1 1 1 -1 -1 -1 1 -1 1 -1 1 -1
47 -1 1 1 1 -1 -1 -1 -1 1 1 1 -1 1
48 1 1 1 1 -1 1 1 1 1 -1 -1 -1 -1
49 -1 -1 -1 -1 1 -1 -1 -1 -1 1 1 1 1
50 1 -1 -1 -1 1 1 1 1 -1 -1 -1 1 -1
51 -1 1 -1 -1 1 1 1 -1 1 -1 1 -1 1
52 1 1 -1 -1 1 -1 -1 1 1 1 -1 -1 -1
53 -1 -1 1 -1 1 1 -1 1 1 1 -1 -1 1
54 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 -1
55 -1 1 1 -1 1 -1 1 1 -1 -1 -1 1 1
56 1 1 1 -1 1 1 -1 -1 -1 1 1 1 -1
57 -1 -1 -1 1 1 -1 1 1 1 1 1 1 -1
58 1 -1 -1 1 1 1 -1 -1 1 -1 -1 1 1
59 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 -1
60 1 1 -1 1 1 -1 1 -1 -1 1 -1 -1 1
61 -1 -1 1 1 1 1 1 -1 -1 1 -1 -1 -1
62 1 -1 1 1 1 -1 -1 1 -1 -1 1 -1 1
63 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 -1
64 1 1 1 1 1 1 1 1 1 1 1 1 1
B:D:E C:D:E A:B:F A:C:F B:C:F A:D:F B:D:F C:D:F A:E:F B:E:F C:E:F D:E:F
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
3 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
5 -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
7 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
9 1 1 -1 -1 -1 1 1 1 -1 -1 -1 1
10 1 1 1 1 -1 -1 1 1 1 -1 -1 1
11 -1 1 1 -1 1 1 -1 1 -1 1 -1 1
12 -1 1 -1 1 1 -1 -1 1 1 1 -1 1
13 1 -1 -1 1 1 1 1 -1 -1 -1 1 1
14 1 -1 1 -1 1 -1 1 -1 1 -1 1 1
15 -1 -1 1 1 -1 1 -1 -1 -1 1 1 1
16 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1
17 1 1 -1 -1 -1 -1 -1 -1 1 1 1 1
18 1 1 1 1 -1 1 -1 -1 -1 1 1 1
19 -1 1 1 -1 1 -1 1 -1 1 -1 1 1
20 -1 1 -1 1 1 1 1 -1 -1 -1 1 1
21 1 -1 -1 1 1 -1 -1 1 1 1 -1 1
22 1 -1 1 -1 1 1 -1 1 -1 1 -1 1
23 -1 -1 1 1 -1 -1 1 1 1 -1 -1 1
24 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1 1
25 -1 -1 -1 -1 -1 1 1 1 1 1 1 -1
26 -1 -1 1 1 -1 -1 1 1 -1 1 1 -1
27 1 -1 1 -1 1 1 -1 1 1 -1 1 -1
28 1 -1 -1 1 1 -1 -1 1 -1 -1 1 -1
29 -1 1 -1 1 1 1 1 -1 1 1 -1 -1
30 -1 1 1 -1 1 -1 1 -1 -1 1 -1 -1
31 1 1 1 1 -1 1 -1 -1 1 -1 -1 -1
32 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
33 -1 -1 1 1 1 1 1 1 1 1 1 1
34 -1 -1 -1 -1 1 -1 1 1 -1 1 1 1
35 1 -1 -1 1 -1 1 -1 1 1 -1 1 1
36 1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1
37 -1 1 1 -1 -1 1 1 -1 1 1 -1 1
38 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 1
39 1 1 -1 -1 1 1 -1 -1 1 -1 -1 1
40 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1
41 1 1 1 1 1 -1 -1 -1 1 1 1 -1
42 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1
43 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1
44 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1
45 1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1
46 1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1
47 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1
48 -1 -1 1 1 1 1 1 1 -1 -1 -1 -1
49 1 1 1 1 1 1 1 1 -1 -1 -1 -1
50 1 1 -1 -1 1 -1 1 1 1 -1 -1 -1
51 -1 1 -1 1 -1 1 -1 1 -1 1 -1 -1
52 -1 1 1 -1 -1 -1 -1 1 1 1 -1 -1
53 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1
54 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1
55 -1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1
56 -1 -1 1 1 1 -1 -1 -1 1 1 1 -1
57 -1 -1 1 1 1 -1 -1 -1 -1 -1 -1 1
58 -1 -1 -1 -1 1 1 -1 -1 1 -1 -1 1
59 1 -1 -1 1 -1 -1 1 -1 -1 1 -1 1
60 1 -1 1 -1 -1 1 1 -1 1 1 -1 1
61 -1 1 1 -1 -1 -1 -1 1 -1 -1 1 1
62 -1 1 -1 1 -1 1 -1 1 1 -1 1 1
63 1 1 -1 -1 1 -1 1 1 -1 1 1 1
64 1 1 1 1 1 1 1 1 1 1 1 1
A:B:C:D A:B:C:E A:B:D:E A:C:D:E B:C:D:E A:B:C:F A:B:D:F A:C:D:F B:C:D:F
1 1 1 1 1 1 1 1 1 1
2 -1 -1 -1 -1 1 -1 -1 -1 1
3 -1 -1 -1 1 -1 -1 -1 1 -1
4 1 1 1 -1 -1 1 1 -1 -1
5 -1 -1 1 -1 -1 -1 1 -1 -1
6 1 1 -1 1 -1 1 -1 1 -1
7 1 1 -1 -1 1 1 -1 -1 1
8 -1 -1 1 1 1 -1 1 1 1
9 -1 1 -1 -1 -1 1 -1 -1 -1
10 1 -1 1 1 -1 -1 1 1 -1
11 1 -1 1 -1 1 -1 1 -1 1
12 -1 1 -1 1 1 1 -1 1 1
13 1 -1 -1 1 1 -1 -1 1 1
14 -1 1 1 -1 1 1 1 -1 1
15 -1 1 1 1 -1 1 1 1 -1
16 1 -1 -1 -1 -1 -1 -1 -1 -1
17 1 -1 -1 -1 -1 1 1 1 1
18 -1 1 1 1 -1 -1 -1 -1 1
19 -1 1 1 -1 1 -1 -1 1 -1
20 1 -1 -1 1 1 1 1 -1 -1
21 -1 1 -1 1 1 -1 1 -1 -1
22 1 -1 1 -1 1 1 -1 1 -1
23 1 -1 1 1 -1 1 -1 -1 1
24 -1 1 -1 -1 -1 -1 1 1 1
25 -1 -1 1 1 1 1 -1 -1 -1
26 1 1 -1 -1 1 -1 1 1 -1
27 1 1 -1 1 -1 -1 1 -1 1
28 -1 -1 1 -1 -1 1 -1 1 1
29 1 1 1 -1 -1 -1 -1 1 1
30 -1 -1 -1 1 -1 1 1 -1 1
31 -1 -1 -1 -1 1 1 1 1 -1
32 1 1 1 1 1 -1 -1 -1 -1
33 1 1 1 1 1 -1 -1 -1 -1
34 -1 -1 -1 -1 1 1 1 1 -1
35 -1 -1 -1 1 -1 1 1 -1 1
36 1 1 1 -1 -1 -1 -1 1 1
37 -1 -1 1 -1 -1 1 -1 1 1
38 1 1 -1 1 -1 -1 1 -1 1
39 1 1 -1 -1 1 -1 1 1 -1
40 -1 -1 1 1 1 1 -1 -1 -1
41 -1 1 -1 -1 -1 -1 1 1 1
42 1 -1 1 1 -1 1 -1 -1 1
43 1 -1 1 -1 1 1 -1 1 -1
44 -1 1 -1 1 1 -1 1 -1 -1
45 1 -1 -1 1 1 1 1 -1 -1
46 -1 1 1 -1 1 -1 -1 1 -1
47 -1 1 1 1 -1 -1 -1 -1 1
48 1 -1 -1 -1 -1 1 1 1 1
49 1 -1 -1 -1 -1 -1 -1 -1 -1
50 -1 1 1 1 -1 1 1 1 -1
51 -1 1 1 -1 1 1 1 -1 1
52 1 -1 -1 1 1 -1 -1 1 1
53 -1 1 -1 1 1 1 -1 1 1
54 1 -1 1 -1 1 -1 1 -1 1
55 1 -1 1 1 -1 -1 1 1 -1
56 -1 1 -1 -1 -1 1 -1 -1 -1
57 -1 -1 1 1 1 -1 1 1 1
58 1 1 -1 -1 1 1 -1 -1 1
59 1 1 -1 1 -1 1 -1 1 -1
60 -1 -1 1 -1 -1 -1 1 -1 -1
61 1 1 1 -1 -1 1 1 -1 -1
62 -1 -1 -1 1 -1 -1 -1 1 -1
63 -1 -1 -1 -1 1 -1 -1 -1 1
64 1 1 1 1 1 1 1 1 1
A:B:E:F A:C:E:F B:C:E:F A:D:E:F B:D:E:F C:D:E:F A:B:C:D:E A:B:C:D:F
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
9 1 1 1 -1 -1 -1 1 1
10 -1 -1 1 1 -1 -1 -1 -1
11 -1 1 -1 -1 1 -1 -1 -1
12 1 -1 -1 1 1 -1 1 1
13 1 -1 -1 -1 -1 1 -1 -1
14 -1 1 -1 1 -1 1 1 1
15 -1 -1 1 -1 1 1 1 1
16 1 1 1 1 1 1 -1 -1
17 -1 -1 -1 -1 -1 -1 1 -1
18 1 1 -1 1 -1 -1 -1 1
19 1 -1 1 -1 1 -1 -1 1
20 -1 1 1 1 1 -1 1 -1
21 -1 1 1 -1 -1 1 -1 1
22 1 -1 1 1 -1 1 1 -1
23 1 1 -1 -1 1 1 1 -1
24 -1 -1 -1 1 1 1 -1 1
25 -1 -1 -1 1 1 1 -1 1
26 1 1 -1 -1 1 1 1 -1
27 1 -1 1 1 -1 1 1 -1
28 -1 1 1 -1 -1 1 -1 1
29 -1 1 1 1 1 -1 1 -1
30 1 -1 1 -1 1 -1 -1 1
31 1 1 -1 1 -1 -1 -1 1
32 -1 -1 -1 -1 -1 -1 1 -1
33 -1 -1 -1 -1 -1 -1 -1 1
34 1 1 -1 1 -1 -1 1 -1
35 1 -1 1 -1 1 -1 1 -1
36 -1 1 1 1 1 -1 -1 1
37 -1 1 1 -1 -1 1 1 -1
38 1 -1 1 1 -1 1 -1 1
39 1 1 -1 -1 1 1 -1 1
40 -1 -1 -1 1 1 1 1 -1
41 -1 -1 -1 1 1 1 1 -1
42 1 1 -1 -1 1 1 -1 1
43 1 -1 1 1 -1 1 -1 1
44 -1 1 1 -1 -1 1 1 -1
45 -1 1 1 1 1 -1 -1 1
46 1 -1 1 -1 1 -1 1 -1
47 1 1 -1 1 -1 -1 1 -1
48 -1 -1 -1 -1 -1 -1 -1 1
49 1 1 1 1 1 1 1 1
50 -1 -1 1 -1 1 1 -1 -1
51 -1 1 -1 1 -1 1 -1 -1
52 1 -1 -1 -1 -1 1 1 1
53 1 -1 -1 1 1 -1 -1 -1
54 -1 1 -1 -1 1 -1 1 1
55 -1 -1 1 1 -1 -1 1 1
56 1 1 1 -1 -1 -1 -1 -1
57 1 1 1 -1 -1 -1 -1 -1
58 -1 -1 1 1 -1 -1 1 1
59 -1 1 -1 -1 1 -1 1 1
60 1 -1 -1 1 1 -1 -1 -1
61 1 -1 -1 -1 -1 1 1 1
62 -1 1 -1 1 -1 1 -1 -1
63 -1 -1 1 -1 1 1 -1 -1
64 1 1 1 1 1 1 1 1
A:B:C:E:F A:B:D:E:F A:C:D:E:F B:C:D:E:F A:B:C:D:E:F
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
17 1 1 1 1 -1
18 -1 -1 -1 1 1
19 -1 -1 1 -1 1
20 1 1 -1 -1 -1
21 -1 1 -1 -1 1
22 1 -1 1 -1 -1
23 1 -1 -1 1 -1
24 -1 1 1 1 1
25 1 -1 -1 -1 1
26 -1 1 1 -1 -1
27 -1 1 -1 1 -1
28 1 -1 1 1 1
29 -1 -1 1 1 -1
30 1 1 -1 1 1
31 1 1 1 -1 1
32 -1 -1 -1 -1 -1
33 1 1 1 1 -1
34 -1 -1 -1 1 1
35 -1 -1 1 -1 1
36 1 1 -1 -1 -1
37 -1 1 -1 -1 1
38 1 -1 1 -1 -1
39 1 -1 -1 1 -1
40 -1 1 1 1 1
41 1 -1 -1 -1 1
42 -1 1 1 -1 -1
43 -1 1 -1 1 -1
44 1 -1 1 1 1
45 -1 -1 1 1 -1
46 1 1 -1 1 1
47 1 1 1 -1 1
48 -1 -1 -1 -1 -1
49 -1 -1 -1 -1 1
50 1 1 1 -1 -1
51 1 1 -1 1 -1
52 -1 -1 1 1 1
53 1 -1 1 1 -1
54 -1 1 -1 1 1
55 -1 1 1 -1 1
56 1 -1 -1 -1 -1
57 -1 1 1 1 -1
58 1 -1 -1 1 1
59 1 -1 1 -1 1
60 -1 1 -1 -1 -1
61 1 1 -1 -1 1
62 -1 -1 1 -1 -1
63 -1 -1 -1 1 -1
64 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 32 33 34 35 36 37 38 39 40 41 42 43 44 45
[47] 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63##----------------------------------------------------------------------
## Estimativa dos efeitos
## Podemos usar as colunas da tabela de sinais para calcular os
## contrastes
y <- dados$y
## Contrastes
(contr <- t(tab[, -1]) %*% y) [,1]
A 12
B 48
C -30
D -96
E -46
F -284
A:B 30
A:C 24
B:C 8
A:D 34
B:D 26
C:D 12
A:E 68
B:E 40
C:E 78
D:E 36
A:F 30
B:F 34
C:F 4
D:F -10
E:F -60
A:B:C 2
A:B:D 0
A:C:D -38
B:C:D -82
A:B:E -50
A:C:E 52
B:C:E -28
A:D:E 38
B:D:E 14
C:D:E 12
A:B:F 56
A:C:F -30
B:C:F -22
A:D:F 40
B:D:F 32
C:D:F -30
A:E:F 30
B:E:F -14
C:E:F 0
D:E:F 42
A:B:C:D -48
A:B:C:E -50
A:B:D:E -76
A:C:D:E 2
B:C:D:E -26
A:B:C:F -12
A:B:D:F 62
A:C:D:F 8
B:C:D:F -92
A:B:E:F 40
A:C:E:F 70
B:C:E:F -2
A:D:E:F -76
B:D:E:F -52
C:D:E:F -30
A:B:C:D:E -16
A:B:C:D:F -90
A:B:C:E:F -16
A:B:D:E:F 2
A:C:D:E:F -8
B:C:D:E:F -12
A:B:C:D:E:F -90## Efeitos = contraste/(n2^{k-1})
(ef <- contr/(n * 2^(k - 1))) [,1]
A 0.3750
B 1.5000
C -0.9375
D -3.0000
E -1.4375
F -8.8750
A:B 0.9375
A:C 0.7500
B:C 0.2500
A:D 1.0625
B:D 0.8125
C:D 0.3750
A:E 2.1250
B:E 1.2500
C:E 2.4375
D:E 1.1250
A:F 0.9375
B:F 1.0625
C:F 0.1250
D:F -0.3125
E:F -1.8750
A:B:C 0.0625
A:B:D 0.0000
A:C:D -1.1875
B:C:D -2.5625
A:B:E -1.5625
A:C:E 1.6250
B:C:E -0.8750
A:D:E 1.1875
B:D:E 0.4375
C:D:E 0.3750
A:B:F 1.7500
A:C:F -0.9375
B:C:F -0.6875
A:D:F 1.2500
B:D:F 1.0000
C:D:F -0.9375
A:E:F 0.9375
B:E:F -0.4375
C:E:F 0.0000
D:E:F 1.3125
A:B:C:D -1.5000
A:B:C:E -1.5625
A:B:D:E -2.3750
A:C:D:E 0.0625
B:C:D:E -0.8125
A:B:C:F -0.3750
A:B:D:F 1.9375
A:C:D:F 0.2500
B:C:D:F -2.8750
A:B:E:F 1.2500
A:C:E:F 2.1875
B:C:E:F -0.0625
A:D:E:F -2.3750
B:D:E:F -1.6250
C:D:E:F -0.9375
A:B:C:D:E -0.5000
A:B:C:D:F -2.8125
A:B:C:E:F -0.5000
A:B:D:E:F 0.0625
A:C:D:E:F -0.2500
B:C:D:E:F -0.3750
A:B:C:D:E:F -2.8125##----------------------------------------------------------------------
## 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
m0 <- lm(y ~ (A + B + C + D + E + F)^4, data = dados)
anova(m0)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
A 1 2.25 2.25 0.0596 0.8141907
B 1 36.00 36.00 0.9530 0.3614869
C 1 14.06 14.06 0.3723 0.5610556
D 1 144.00 144.00 3.8119 0.0918423 .
E 1 33.06 33.06 0.8752 0.3806600
F 1 1260.25 1260.25 33.3604 0.0006801 ***
A:B 1 14.06 14.06 0.3723 0.5610556
A:C 1 9.00 9.00 0.2382 0.6403960
A:D 1 18.06 18.06 0.4781 0.5115413
A:E 1 72.25 72.25 1.9126 0.2091841
A:F 1 14.06 14.06 0.3723 0.5610556
B:C 1 1.00 1.00 0.0265 0.8753519
B:D 1 10.56 10.56 0.2796 0.6133072
B:E 1 25.00 25.00 0.6618 0.4427188
B:F 1 18.06 18.06 0.4781 0.5115413
C:D 1 2.25 2.25 0.0596 0.8141907
C:E 1 95.06 95.06 2.5164 0.1566863
C:F 1 0.25 0.25 0.0066 0.9374407
D:E 1 20.25 20.25 0.5360 0.4878603
D:F 1 1.56 1.56 0.0414 0.8446276
E:F 1 56.25 56.25 1.4890 0.2618788
A:B:C 1 0.06 0.06 0.0017 0.9686908
A:B:D 1 0.00 0.00 0.0000 1.0000000
A:B:E 1 39.06 39.06 1.0340 0.3430644
A:B:F 1 49.00 49.00 1.2971 0.2922155
A:C:D 1 22.56 22.56 0.5973 0.4649161
A:C:E 1 42.25 42.25 1.1184 0.3253850
A:C:F 1 14.06 14.06 0.3723 0.5610556
A:D:E 1 22.56 22.56 0.5973 0.4649161
A:D:F 1 25.00 25.00 0.6618 0.4427188
A:E:F 1 14.06 14.06 0.3723 0.5610556
B:C:D 1 105.06 105.06 2.7811 0.1393158
B:C:E 1 12.25 12.25 0.3243 0.5868510
B:C:F 1 7.56 7.56 0.2002 0.6680857
B:D:E 1 3.06 3.06 0.0811 0.7840941
B:D:F 1 16.00 16.00 0.4235 0.5359454
B:E:F 1 3.06 3.06 0.0811 0.7840941
C:D:E 1 2.25 2.25 0.0596 0.8141907
C:D:F 1 14.06 14.06 0.3723 0.5610556
C:E:F 1 0.00 0.00 0.0000 1.0000000
D:E:F 1 27.56 27.56 0.7296 0.4212751
A:B:C:D 1 36.00 36.00 0.9530 0.3614869
A:B:C:E 1 39.06 39.06 1.0340 0.3430644
A:B:C:F 1 2.25 2.25 0.0596 0.8141907
A:B:D:E 1 90.25 90.25 2.3890 0.1661071
A:B:D:F 1 60.06 60.06 1.5899 0.2477344
A:B:E:F 1 25.00 25.00 0.6618 0.4427188
A:C:D:E 1 0.06 0.06 0.0017 0.9686908
A:C:D:F 1 1.00 1.00 0.0265 0.8753519
A:C:E:F 1 76.56 76.56 2.0267 0.1975639
A:D:E:F 1 90.25 90.25 2.3890 0.1661071
B:C:D:E 1 10.56 10.56 0.2796 0.6133072
B:C:D:F 1 132.25 132.25 3.5008 0.1035181
B:C:E:F 1 0.06 0.06 0.0017 0.9686908
B:D:E:F 1 42.25 42.25 1.1184 0.3253850
C:D:E:F 1 14.06 14.06 0.3723 0.5610556
Residuals 7 264.44 37.78
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0, . ~ (A + B + C + D + E + F)^3)
anova(m1)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
A 1 2.25 2.25 0.0560 0.8151
B 1 36.00 36.00 0.8958 0.3542
C 1 14.06 14.06 0.3499 0.5602
D 1 144.00 144.00 3.5832 0.0716 .
E 1 33.06 33.06 0.8227 0.3742
F 1 1260.25 1260.25 31.3593 1.249e-05 ***
A:B 1 14.06 14.06 0.3499 0.5602
A:C 1 9.00 9.00 0.2240 0.6407
A:D 1 18.06 18.06 0.4495 0.5096
A:E 1 72.25 72.25 1.7978 0.1937
A:F 1 14.06 14.06 0.3499 0.5602
B:C 1 1.00 1.00 0.0249 0.8761
B:D 1 10.56 10.56 0.2628 0.6133
B:E 1 25.00 25.00 0.6221 0.4387
B:F 1 18.06 18.06 0.4495 0.5096
C:D 1 2.25 2.25 0.0560 0.8151
C:E 1 95.06 95.06 2.3655 0.1383
C:F 1 0.25 0.25 0.0062 0.9378
D:E 1 20.25 20.25 0.5039 0.4853
D:F 1 1.56 1.56 0.0389 0.8455
E:F 1 56.25 56.25 1.3997 0.2494
A:B:C 1 0.06 0.06 0.0016 0.9689
A:B:D 1 0.00 0.00 0.0000 1.0000
A:B:E 1 39.06 39.06 0.9720 0.3349
A:B:F 1 49.00 49.00 1.2193 0.2814
A:C:D 1 22.56 22.56 0.5614 0.4616
A:C:E 1 42.25 42.25 1.0513 0.3163
A:C:F 1 14.06 14.06 0.3499 0.5602
A:D:E 1 22.56 22.56 0.5614 0.4616
A:D:F 1 25.00 25.00 0.6221 0.4387
A:E:F 1 14.06 14.06 0.3499 0.5602
B:C:D 1 105.06 105.06 2.6143 0.1202
B:C:E 1 12.25 12.25 0.3048 0.5864
B:C:F 1 7.56 7.56 0.1882 0.6687
B:D:E 1 3.06 3.06 0.0762 0.7851
B:D:F 1 16.00 16.00 0.3981 0.5346
B:E:F 1 3.06 3.06 0.0762 0.7851
C:D:E 1 2.25 2.25 0.0560 0.8151
C:D:F 1 14.06 14.06 0.3499 0.5602
C:E:F 1 0.00 0.00 0.0000 1.0000
D:E:F 1 27.56 27.56 0.6858 0.4165
Residuals 22 884.12 40.19
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2 <- update(m1, . ~ (A + B + C + D + E + F)^2)
anova(m2)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
A 1 2.25 2.25 0.0725 0.78906
B 1 36.00 36.00 1.1598 0.28765
C 1 14.06 14.06 0.4531 0.50457
D 1 144.00 144.00 4.6394 0.03703 *
E 1 33.06 33.06 1.0652 0.30794
F 1 1260.25 1260.25 40.6026 1.158e-07 ***
A:B 1 14.06 14.06 0.4531 0.50457
A:C 1 9.00 9.00 0.2900 0.59309
A:D 1 18.06 18.06 0.5819 0.44982
A:E 1 72.25 72.25 2.3277 0.13458
A:F 1 14.06 14.06 0.4531 0.50457
B:C 1 1.00 1.00 0.0322 0.85841
B:D 1 10.56 10.56 0.3403 0.56278
B:E 1 25.00 25.00 0.8054 0.37459
B:F 1 18.06 18.06 0.5819 0.44982
C:D 1 2.25 2.25 0.0725 0.78906
C:E 1 95.06 95.06 3.0627 0.08741 .
C:F 1 0.25 0.25 0.0081 0.92892
D:E 1 20.25 20.25 0.6524 0.42380
D:F 1 1.56 1.56 0.0503 0.82356
E:F 1 56.25 56.25 1.8123 0.18546
Residuals 42 1303.63 31.04
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m3 <- update(m2, . ~ D + F + (C + E)^2)
anova(m3)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
D 1 144.00 144.00 5.2054 0.02621 *
F 1 1260.25 1260.25 45.5559 7.713e-09 ***
C 1 14.06 14.06 0.5083 0.47872
E 1 33.06 33.06 1.1952 0.27881
C:E 1 95.06 95.06 3.4364 0.06886 .
Residuals 58 1604.50 27.66
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m4 <- update(m3, . ~ D + F)
anova(m4)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
D 1 144.0 144.00 5.0289 0.02857 *
F 1 1260.2 1260.25 44.0120 9.818e-09 ***
Residuals 61 1746.7 28.63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Fazendo o TRV para confirmar o modelo
anova(m4, m3, m2, m1, m0)Analysis of Variance Table
Model 1: y ~ D + F
Model 2: y ~ D + F + C + E + C:E
Model 3: y ~ A + B + C + D + E + F + A:B + A:C + A:D + A:E + A:F + B:C +
B:D + B:E + B:F + C:D + C:E + C:F + D:E + D:F + E:F
Model 4: y ~ A + B + C + D + E + F + A:B + A:C + A:D + A:E + A:F + B:C +
B:D + B:E + B:F + C:D + C:E + C:F + D:E + D:F + E:F + A:B:C +
A:B:D + A:B:E + A:B:F + A:C:D + A:C:E + A:C:F + A:D:E + A:D:F +
A:E:F + B:C:D + B:C:E + B:C:F + B:D:E + B:D:F + B:E:F + C:D:E +
C:D:F + C:E:F + D:E:F
Model 5: y ~ (A + B + C + D + E + F)^4
Res.Df RSS Df Sum of Sq F Pr(>F)
1 61 1746.69
2 58 1604.50 3 142.19 1.2546 0.3607
3 42 1303.63 16 300.88 0.4978 0.8825
4 22 884.12 20 419.50 0.5552 0.8573
5 7 264.44 15 619.69 1.0936 0.4782##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m4)
## Quantis normais
qqnorm(res); qqline(res)
##----------------------------------------------------------------------
## Predições
## Predição para as combinações únicas dos fatores
pred <- expand.grid(D = unique(dados$D),
F = unique(dados$F))
pred$y <- predict(m4, newdata = pred)
pred D F y
1 -1 -1 34.21875
2 1 -1 31.21875
3 -1 1 25.34375
4 1 1 22.34375## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(D = seq(-1, 1, length.out = 30),
F = seq(-1, 1 ,length.out = 30))
pred$y <- predict(m4, newdata = pred)
## Vários formas de visualizar
wireframe(y ~ D + F, data = pred, drape = TRUE)
levelplot(y ~ D + F, data = pred)
levelplot(y ~ D + F, data = pred, cuts = 90,
col.regions = heat.colors)