```{r setup, include=FALSE, purl=FALSE, eval=TRUE} opts_chunk$set( ## knitr options cache = TRUE, tidy = FALSE, comment = NA, fig.width = 7, fig.height = 5, fig.align = "center", #dpi = 60, ## aumentar o dpi para exibir maior. dev = "png" # fig.path = "figures/", ) ``` # 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 [+] ``` ```{r} ## Importação url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-2.txt" dados <- read.table(url, header = TRUE) str(dados) ## 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) summary(m0) ``` Repare 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. ```{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) } ``` ```{r} ## 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) summary(m1) ``` ### Suposições do modelo Como todo modelo de ANOVA: 1. $\epsilon_{ijkl} \sim \text{N}(0, \sigma^2)$ 2. $\epsilon_{ijkl}$ são independentes 3. Homogeneidade de variâncias ### Avaliação do modelo Quando ajustamos um modelo a um conjunto de dados, podemos ter basicamente dois objetivos: 1. Verificar quais variáveis influenciam na média da variável resposta $Y$ 2. Fazer predições com base nos parâmetros estimados 3. Outros: estimar parâmetros e fazer inferência Para verificar as variáves que infuenciam na rugosidade, vemos a ANOVA do modelo, ```{r} anova(m1) ``` Se 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. ```{r} 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: ```{r} anova(m1) anova(m2) ``` Porque se alteram os p-valores? ```{r} ## 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) ## 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) ``` ```{r, eval=FALSE, include=FALSE} ## par(mfrow = c(1,2)) ## plot(seq(0, 10, length=100), type = "l", xlab = "X", ylab = "f(x)", ## y = df(x = seq(0, 10, length=100), df1 = 1, df2 = 8), ## main = expression(list(nu[1] == 1, nu[2] == 8))) ## abline(v = qf(0.05, df1 = anova(m1)$"Df"[1], ## df2 = anova(m1)$"Df"[8], lower.tail = FALSE)) ## plot(seq(0, 10, length=100), type = "l", xlab = "X", ylab = "f(x)", ## y = df(x = seq(0, 10, length=100), df1 = 1, df2 = 12), ## main = expression(list(nu[1] == 1, nu[2] == 12))) ## abline(v = qf(0.05, df1 = anova(m2)$"Df"[1], ## df2 = anova(m2)$"Df"[4], lower.tail = FALSE)) ``` Porque 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. ```{r} anova(m2, m1) ``` Como 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 ```{r} m3 <- update(m2, . ~ Alim.n + Prof.n) anova(m3) ``` Assim como antes, podemos usar o TRV para comparar os modelos encaixados ```{r} anova(m3, m2, m1) ``` Com 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 ```{r} ## 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 ```{r} summary(m3) ``` ```{r} ## 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 ## 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 ``` ```{r} ##---------------------------------------------------------------------- ## 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 ##---------------------------------------------------------------------- ## 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))) ``` Para 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: ```{r} anova(lm(y ~ (A * B * C * D), data = dados)) ``` 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. ```{r} ##---------------------------------------------------------------------- ## 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. ```{r} ##---------------------------------------------------------------------- ## 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) ``` Note 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 ```{r} m1 <- update(m0, . ~ A * D, data = dados) anova(m1) ## Fazendo o TRV para confirmar o modelo anova(m1, m0) ``` #### Análise dos resíduos ```{r} ##---------------------------------------------------------------------- ## 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 ```{r} ##---------------------------------------------------------------------- ## 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 ``` Note como fica a predição para o conjunto de dados ```{r} cbind(dados, predict(m1)) ``` As 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. ```{r} ## 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 ```{r} ##---------------------------------------------------------------------- ## 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)) ``` ```{r, eval=FALSE} ##---------------------------------------------------------------------- ## 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: 1. Fazer um novo experimento, agora com apenas estes dois fatores e com um número maior de repetições: **projeção** 2. 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$ ```{r, eval=FALSE, echo=FALSE, include=FALSE} ## k <- 6 ## input <- replicate(k, c(-1,1), simplify=FALSE) ## names(input) <- LETTERS[1:k] ## dados <- do.call(expand.grid, input) ## dados$y <- c(36,35,37,34,30,32,28,38,40,32,32,28,34,26,30,28, ## 40,30,35,35,32,35,36,32,20,35,35,28,27,36,36,35, ## 32,34,30,24,30,20,32,25,20,20,22,35,26,15,19,24, ## 22,23,22,18,20,20,20,36,20,11,20,35,16,32,20,20) ``` ```{r} ##---------------------------------------------------------------------- ## Dados url <- "http://www.leg.ufpr.br/~fernandomayer/dados/fatorial_2-6.txt" dados <- read.table(url, header = TRUE) str(dados) ## 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 ##---------------------------------------------------------------------- ## 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))) ##---------------------------------------------------------------------- ## 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) m1 <- update(m0, . ~ (A + B + C + D + E + F)^3) anova(m1) m2 <- update(m1, . ~ (A + B + C + D + E + F)^2) anova(m2) m3 <- update(m2, . ~ D + F + (C + E)^2) anova(m3) m4 <- update(m3, . ~ D + F) anova(m4) ## Fazendo o TRV para confirmar o modelo anova(m4, m3, m2, m1, m0) ##---------------------------------------------------------------------- ## 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 ## 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) ```