Aula ENGD02 - Planejamento e Análise de Experimentos
Brenner Silva
2018
Planejamento e Análise de Experimentos \(2^{k}\)
Peças fabricadas em um processo de moldagem por injeção estão apresentando um encolhimento excessivo, que está causando problemas nas operações de montagem antes da área de moldagem por injeção.
Em um esforço de reduzir o encolhimento, uma equipe de melhoria da qualidade decidiu usar um experimento planejado para estudar o processo de moldagem por injeção.
A equipe investigou seis fatores - temperatura do molde (A), velocidade do parafuso (B), tempo de retenção (C), tempo do ciclo (D), tamanho do ponto de injeção (E) e pressão de retenção (F) - cada um com dois níveis, tendo como objetivo aprender como cada fator afeta o encolhimento e obter informações preliminares sobre como os fatores interagem.
MONTGOMERY, D. C.; RUNGER, G. C. Applied Statistics and Probability for Engineers. Third Edition. John Wiley & Son. 2003. EXAMPLE 14-8.
*Peço gentilmente que rodem previamente ao menos este bloco de códigos para que possamos obter bom desenvolvimento da aula.
if(!require("FrF2")) install.packages("FrF2") ; library(FrF2)
if(!require("lattice")) install.packages("lattice") ; library(lattice)
if(!require("scatterplot3d")) install.packages("scatterplot3d") ; library(scatterplot3d)
# Se houver problemas na instalação do FrF2, deve-se instalar primeiro o pacote DoE.base
# if(!require("DoE.base")) install.packages("DoE.base") ; library(DoE.base)
Agora vamos direcionar o script em desenvolvimento a um diretório.
# setwd("Mypath/Myfiles")GERANDO O DESIGN
Sabe-se que o experimento escolhido foi fatorial fracionário, sendo \(p=2\), planejamento básico nos fatores \(A\), \(B\), \(C\) e \(D\), e geradores \(I=ABCE=BCDF\). Assim temos que \(E=ABC\) e \(F=BCD\) e podemos montar nossa matriz de sinais.
plan.injecao <- FrF2(nfactors = 6,
nruns = 16,
replications = 1,
randomize = FALSE,
generators = c("ABC","BCD"),
factor.names = list(
'TM' = c(-1, +1),
'VP' = c(-1, +1),
'TR' = c(-1, +1),
'TC' = c(-1, +1),
'TPI' = c(-1, +1),
'PR' = c(-1, +1)))
# FATORES
# Temperatura do molde (A) => TM
# Velocidade do parafuso (B) => VP
# Tempo de retenção (C) => TR
# Tempo do ciclo (D) => TC
# Tamanho do ponto de injeção (E) => TPI
# Pressão de retenção (F) => PR
Feito e analisado a matriz de sinais, vamos adicionar nossas respostas aos respectivos tratamentos.
encolhimento <- c(6, 10, 32, 60, 4, 15, 26, 60, 8, 12, 34, 60, 16, 5, 37, 52)
new_plan <- add.response(design = plan.injecao, response = encolhimento)ANÁLISES GRÁFICAS INICIAIS
O próximo passo é começar a analisar as respostas obtidas através deste experimento.
# MEPlot(new_plan)
MEPlot(new_plan, abbrev = 5, cex.xax = 1.6, cex.main = 2, main = "Efeitos Principais")
# png('test.png', units="in", width=13, height=7, res=500)
# MEPlot(new_plan, abbrev = 5, cex.xax = 1.6, cex.main = 2)
# dev.off()
# É possível salvar a imagem em outros tipos de extensão, como .tiff. Basta
# substituir o png por tiff na função.
# IAPlot(new_plan)
IAPlot(new_plan, abbrev = 5, show.alias = TRUE, lwd = 2,
cex = 2, cex.xax = 1.2, cex.lab = 1.5, main = "Matriz dos Efeitos de Interação")
Utilizaremos a função xyplot do pacote laticce para gerar os gráficos de interação específico a fim de facilitar a análise gráfica. O uso da função split é para possibilitar múltiplos plots em uma janela.
plot1 <- xyplot(encolhimento ~ new_plan$TM,
groups = new_plan$VP,
type = "a",
ylab = "Encolhimento",
xlab = "Temperatura do molde ",
auto.key = list(space = "right",
points = FALSE,
title = "Velocidade do parafuso",
cex = 0.7,
lines = TRUE))
plot2 <- xyplot(encolhimento ~ new_plan$TR,
groups = new_plan$TPI,
type = "a",
ylab = "Encolhimento",
xlab = "Tempo de retenção",
auto.key = list(space = "right",
points = FALSE,
title = "Tamanho do ponto de injeção",
cex = 0.7,
lines = TRUE))
print(plot1, split = c(1, 1, 1, 2), more = TRUE)
print(plot2, split = c(1, 2, 1, 2))
#(new_plan,interaction.plot(new_plan$TM,
# new_plan$VP,
# new_plan$encolhimento,
# lwd=4,
# col=c("orange","purple"),
# lty=c(2,2)))
Perceba que o pacote laticce é mais arcaico e possui menos ferramentas que o o pacote ggplot2. Entretanto para análise de experimentos ele é extremamente usual em R.
ANÁLISE DA SIGNIFICÂNCIA DOS EFEITOS
Gerando gráficos de probabilidade normal e half-normal (semi-normal) com a função DanielPlot para identificar efeitos possivelmente ativos. O gráfico de probabilidade half normal dos efeitos mostra os valores absolutos (módulo) dos efeitos. Efeitos mais afastados de 0 que estão no eixo x têm maior magnitude. Efeitos mais afastados de 0 são estatisticamente mais significativos. A distância que os pontos devem estar de zero para serem estatisticamente significativos depende do nível de significância.
DanielPlot(new_plan,
code = TRUE,
half = T,
alpha = 0.05,
cex.main = 1.8, cex.sub = 1.2, cex.pch = 1.2, cex.lab = 1.4,
cex.fac = 1.4, cex.axis = 1.2, cex.legend = 1.2, main = "Gráfico half-normal")
# O argumento alpha é o nível de significância. ANOVA
anovaP <- aov(new_plan)
summary(anovaP)## Number of observations used: 16
## Formula:
## encolhimento ~ (TM + VP + TR + TC + TPI + PR)^2
## Df Sum Sq Mean Sq F value Pr(>F)
## TM 1 770 770 16.191 0.05657 .
## VP 1 5077 5077 106.735 0.00924 **
## TR 1 3 3 0.064 0.82339
## TC 1 8 8 0.159 0.72862
## TPI 1 1 1 0.012 0.92333
## PR 1 1 1 0.012 0.92333
## TM:VP 1 564 564 11.859 0.07496 .
## TM:TR 1 11 11 0.222 0.68387
## TM:TC 1 116 116 2.430 0.25939
## TM:TPI 1 14 14 0.296 0.64112
## TM:PR 1 2 2 0.033 0.87288
## VP:TC 1 0 0 0.001 0.97438
## VP:PR 1 0 0 0.001 0.97438
## Residuals 2 95 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Fazendo a estimativa dos coeficientes
# Função lm() é residente do R
lmP <- lm(new_plan)
summary(lmP)## Number of observations used: 16
## Formula:
## encolhimento ~ (TM + VP + TR + TC + TPI + PR)^2
##
## Call:
## lm.default(formula = fo, data = model.frame(fo, data = formula))
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10
## 2.375 -2.375 2.500 -2.500 -2.500 2.500 -2.375 2.375 -2.375 2.375
## 11 12 13 14 15 16
## -2.500 2.500 2.500 -2.500 2.375 -2.375
##
## Coefficients: (8 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.3125 1.7241 15.841 0.00396 **
## TM1 6.9375 1.7241 4.024 0.05657 .
## VP1 17.8125 1.7241 10.331 0.00924 **
## TR1 -0.4375 1.7241 -0.254 0.82339
## TC1 0.6875 1.7241 0.399 0.72862
## TPI1 0.1875 1.7241 0.109 0.92333
## PR1 0.1875 1.7241 0.109 0.92333
## TM1:VP1 5.9375 1.7241 3.444 0.07496 .
## TM1:TR1 -0.8125 1.7241 -0.471 0.68387
## TM1:TC1 -2.6875 1.7241 -1.559 0.25939
## TM1:TPI1 -0.9375 1.7241 -0.544 0.64112
## TM1:PR1 0.3125 1.7241 0.181 0.87288
## VP1:TR1 NA NA NA NA
## VP1:TC1 -0.0625 1.7241 -0.036 0.97438
## VP1:TPI1 NA NA NA NA
## VP1:PR1 -0.0625 1.7241 -0.036 0.97438
## TR1:TC1 NA NA NA NA
## TR1:TPI1 NA NA NA NA
## TR1:PR1 NA NA NA NA
## TC1:TPI1 NA NA NA NA
## TC1:PR1 NA NA NA NA
## TPI1:PR1 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.897 on 2 degrees of freedom
## Multiple R-squared: 0.9857, Adjusted R-squared: 0.8929
## F-statistic: 10.62 on 13 and 2 DF, p-value: 0.08928ANÁLISE DE RESÍDUOS
par(mfrow = c(2,2))
plot(lmP)
A função residente par e mfrow me permite gera uma espécie de matriz 2x2 na window do R, desta maneira, rodando primeiro o código par(mfrow=c(2,2)) e posteriormente o plot, temos 4 gráficos de 1x.
lmP2 <- lm(encolhimento ~ new_plan$VP, data = new_plan)
summary(lmP2)##
## Call:
## lm.default(formula = encolhimento ~ new_plan$VP, data = new_plan)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.125 -6.156 -0.500 6.594 14.875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.312 2.658 10.275 6.67e-08 ***
## new_plan$VP1 17.812 2.658 6.701 1.01e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.63 on 14 degrees of freedom
## Multiple R-squared: 0.7623, Adjusted R-squared: 0.7453
## F-statistic: 44.9 on 1 and 14 DF, p-value: 1.009e-05par(mfrow=c(1,3))
plot(lmP2)
A linearidade dos residuos pode ser verificada visualmente através do gráfico Residuals vs Fitted.
O QQ Plot (*Normal Q-Q) nos auxilia a verificar a normalidade dos resíduos, no nosso caso os resíduos seguem bem alinhados na linha tracejada reta, indicando a normalidade.
O gráfico de Scale-Location , mostra se os resíduos são distribuídos igualmente ao longo dos intervalos de preditores, possibilitando a análise de homo ou heterocedasticidade.
* Maiores discussões sobre estes gráficos você pode encontrar aqui ou aqui.
EFEITOS
efeitos <- na.omit(2*coef(lmP)[-1])
print(efeitos)## TM1 VP1 TR1 TC1 TPI1 PR1 TM1:VP1 TM1:TR1
## 13.875 35.625 -0.875 1.375 0.375 0.375 11.875 -1.625
## TM1:TC1 TM1:TPI1 TM1:PR1 VP1:TC1 VP1:PR1
## -5.375 -1.875 0.625 -0.125 -0.125REGRESSÃO
VP = seq(-1, 1, 0.1)
TM = seq(-1, 1, length.out = length(VP))
f = function(VP, TM) 27.312+6.9375*TM+17.812*VP+5.9375*TM*VP
grid <- expand.grid(VP = VP,
TM = TM)
grid[, "fit"] <- f(grid$VP, grid$TM)ANÁLISE GRÁFICA 2
A nível de ilustração e para melhor análise da interação dos fatores, realizaremos plots em 3D e superficies de resposta.
Inicialmente utilizaremos o scatterplot3d para plotar um gráfico de linha e pontos em 3D, e nosso superfície de contorno usaremos a função contourplotdo pacote lattice, além do wireframe para superfície em 3D.
scatterplot3d(x = TM, y = VP, z = f(VP, TM),
main = "Análise de Encolhimento",
zlab = "Encolhimento",
color = "blue",
type = "o")
scatterplot3d(x = TM, y = VP, z = f(VP, TM),
main = "Análise de Encolhimento",
zlab = "Encolhimento",
color = "#56B4E9",
pch = 16,
type = "o",
grid = F,
box = T)
contourplot(fit ~ TM * VP, data = grid,
cuts = 20,
region = TRUE,
xlab = "TM",
ylab = "VP",
main = "An. Encolhimento")
# levelplot é uma função similar a contourplot.
O R tem uma paleta de cores residente que nos permite alterar a estética do nosso gráfico. Para promover a modificação, basta acrescentar o argumento col.regionsa função do nosso gráfico.
# Paleta de cores residente do R - col.regions
# heat.colors(100, alpha = 1)
# col.redions =
# rainbow(xx)
# terrain.colors(100, alpha = 1)
# topo.colors(100, alpha = 1)
# cm.colors(100, alpha = 1)
contourplot(fit ~ TM * VP, data = grid,
cuts = 20,
region = TRUE,
xlab = "TM",
ylab = "VP",
main = "An. Encolhimento",
col.regions = heat.colors(100, alpha = 1))
contourplot(fit ~ TM * VP, data = grid,
cuts = 20,
region = TRUE,
xlab = "TM",
ylab = "VP",
main = "An. Encolhimento",
col.regions = terrain.colors(100, alpha = 1))
Com a função wireframe é possível gerar nosso modelo de superfície de resposta.
wireframe(fit ~ TM * VP, data = grid,
panel.aspect = 0.9,
zoom = 0.7,
scales = list(arrows = FALSE),
drape = TRUE,
colorkey = TRUE)
wireframe(fit ~ TM * VP, data = grid,
panel.aspect = 0.9,
zoom = 0.7,
screen = list(z = 330, x = -75),
scales = list(arrows = FALSE),
drape = TRUE,
colorkey = TRUE)
wireframe(fit ~ TM * VP, data = grid,
panel.aspect = 0.5,
zoom = 0.8,
scales = list(arrows = FALSE),
drape = TRUE,
col.regions = heat.colors(100, alpha = 1))
wireframe(fit ~ TM * VP, data = grid,
panel.aspect = 0.5,
zoom = 0.8,
screen = list(z = 330, x = -75),
scales = list(arrows = F, cex = .8, col = "black", font = 3),
drape = TRUE,
xlab = expression(paste('TM')),
ylab = expression(paste('VP')),
zlab = list(label = "An. Encolhimento", font = 12, cex = .9),
col.regions = topo.colors(100, alpha = 1))
wireframe(fit ~ TM * VP, data = grid,
panel.aspect = 0.5,
zoom = 0.8,
screen = list(z = 330, x = -75),
scales = list(arrows = F, cex = .9, col = "black", font = 3),
drape = TRUE,
xlab = expression(paste('TM')),
ylab = expression(paste('VP')),
zlab = list(label = "An. Encolhimento", font = 14, cex = 1),
col.regions = rainbow(100,
s = 1,
v = 1,
start = 0,
end = max(1,100 - 1)/100,
alpha = 1))
# Para a presente função o argumento screen nos permite rotacionar a superfície.
# As funções residentes persp, contour e varfcn também gera bons gráficos para análise.DICAS
Além dos pacotes aqui utilizados, existem diversos outros pacotes no R que possibilitam o planejamento e análise estatística experimental, alguns dos outros pacotes são:
* Pacotes
+ DoE.base;
+ BsMD;
+ pid;
+ AlgDesign;
+ skpr;
+ BHH2;
+ unrepx;
+ conf.design;
+ qcc;
+ agricolae;
+ DoE.wrapper;
+ DiceDesign;
+ DiceEval;
+ rsm;
+ daewr.
Além dos livros de Montgomery, vale a pena consultar também os livros Design and Analysis of Experiments with R de Lawson (2015), Design of Experiments for Engineers and Scientists de Antony (2014), Lattice: Multivariate Data Visualization with R de Sarkar (2008), A First Course in Design and Analysis of Experiments de Oehlert (2010) e Regression Modeling Strategies de Harrell (2013). Slides e notas de aula, como as do professor Zocchi (2015a, 2015b) também são de grande ajuda.
Brenner Biasi S. Silva
