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 \(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 ' ' 1
summary(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.02881

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.

## 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 ' ' 1
summary(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.02881

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,

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 ' ' 1

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.

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 ' ' 1
anova(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 ' ' 1

Porque 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.0009900542

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.

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.4609

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

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 ' ' 1

Assim 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.4609

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

## 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)
       })