Aula Prática DoE - PEI
ENGD83 Matemática Instrumental para Engenharia
Aula prática sobre Planejamento e Análise de Experimentos da disciplina Matemática Instrumental para Engenharia (ENGD83) do programa de pós-graduação em Engenharia Industrial (PEI) da UFBA desenvolvida pelo discente Brenner Biasi S. Silva sob supervisão da docente Karla P. Esquerre.
Experimentos Fatoriais \(2^{k}\)
Um experimento fatorial \(2^{k}\) é uma metodologia experimental que permite observar os efeitos que vários fatores podem ter sobre uma resposta através de combinações de nível de fator de dois ou mais fatores, variando os níveis de todos os fatores ao mesmo tempo, permitindo assim estudar as interações (efeitos sinérgicos) entre os fatores.
A estrutura do tratamento fatorial proporciona os benefícios de ser eficiente e de uma análise simplificada, mas a estruturação dos k fatores pode se tornar o procedimento volumoso. Por exemplo, sete fatores, o experimento tem \(2^{7}=128\) tratamentos. Haverá 127 graus de liberdade para a modelagem, sendo 7 graus de liberdade para efeitos principais, 21 graus de liberdade para interações de dois fatores, 35 graus de liberdade para interações de três fatores. Percebendo-se assim que experimentos fatoriais “demandam” muitos graus de liberdade em interações de elevadas ordens.
Elevadas quantidades de interações de alta ordem implicam em muitos graus de liberdade, podendo aumentar a possibilidade de obtermos erros. Em alguns experimentos pode ser inviável a continuidade do modelo full dado a possibilidade de erro, o custo, desinteresse em interações de elevada ordem e/ou quando sabemos a possibilidade de negligencia-las.
Um experimento fatorial fracionário (\(2^{k-p}\)) é uma modificação de um fatorial padrão full, no qual é possível obter informações sobre os principais efeitos e interações de baixa ordem sem perda significativa de informação acerca dos efeitos principais e interações de baixa ordem do experimento full. Sua fundamentação é baseada no confounding, análogo ao procedimento com blocos porém extraindo informação de apenas um dos blocos, da ideia de que é impossível executar todas as observações em um experimento sob condições homogêneas.
Uma meia-fração de um experimento \(2^{k}\) contém \(2^{k-1}\) corridas e é chamado de um planejamento fatorial fracionário \(2^{k-1}\). Do modelo padrão, considerando um experimento \(2^{5}\), os tratamentos a serem analisados serão formados selecionando-se apenas as combinações de tratamento que geram um sinal positivo no efeito de maior ordem de um experimento \(2^{k}\), sendo assim a interação ABCDE é chamada de gerador desta fração particular. O elemento identidade \(I\) tem também o sinal mais (\(+\)), assim têm-se a relação de definição para o planejamento \(I=ABCDE\).
Aula Prática
Nesta aula desenvolveremos a análise e interpretação de um experimento no qual buscou-se analisar o desempenho de um novo material para adsorção de óleo em água. O presente experimento é referente ao artigo Oil removal from water by fungal biomass: A factorial design analysis e foi desenvolvido originalmente no software Minitab.
Contexto
A presença de óleo emulsionado nas águas residuais é uma preocupação real, pois muitas vezes resulta em incrustação do equipamento de processamento e se torna parte da carga orgânica que deve ser tratada no tratamento biológico de tais águas residuais. Convencionalmente, técnicas físico-químicas são usados para a remoção de óleo e graxa.
Objetivo
O objetivo do presente estudo foi realizar uma análise de planejamento fatorial para analisar os fatores significativos que influenciaram a remoção de óleo da água por adsorção para aplicação de um novo biomaterial e entender seu impacto sobre o processo.
Planejamento e Análise de Experimentos \(2^{k}\)
Para a presente aula prática serão utilizados os seguintes pacotes abaixo.
*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("qualityTools")) install.packages("qualityTools") ; library(qualityTools)
if(!require("lattice")) install.packages("lattice") ; library(lattice)
# 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")
O planejamento escolhido foi o experimento fatorial fracionado \(2^{k-p}\) com cinco fatores, sendo \(p =1\), e com uma repetição. Deste modo, o experimento pode ser classificado como Planejamento de Resolução V, que pode ser ratificado na tabela abaixo.
O presente planejamento busca analisar a influência dos fatores pH da solução, temperatura, dose de adsorvente, concentração do óleo, velocidade de rotação do agitador e seus efeitos sinérgicos de acordo com a tabela mostrada abaixo.
Gerando o design
Dado o conhecimento sobre a análise experimental a ser desenvolvida, utilizaremos a função plan.person do pacote FrF2 (link) para gerar o design.
A função comentada é a maneira geral para esta análise. Entretanto, é possível personaliza-lá de modo a incluir o nome e level dos fatores em análise.
# plan.person <- FrF2(nfactors = 5,
# resolution = 3,
# replications = 1,
# randomize = FALSE)
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
plan.person <- FrF2(nfactors = 5,
resolution = 5,
replications = 2,
randomize = FALSE,
factor.names = list(
pH = c(3, 9),
Temp = c(5, 30),
Dose = c(0.05, 0.5),
Conc = c(50, 350),
Velo = c(100, 200)))
print(plan.person)## run.no run.no.std.rp pH Temp Dose Conc Velo
## 1 1 1.1 3 5 0.05 50 200
## 2 2 2.1 9 5 0.05 50 100
## 3 3 3.1 3 30 0.05 50 100
## 4 4 4.1 9 30 0.05 50 200
## 5 5 5.1 3 5 0.5 50 100
## 6 6 6.1 9 5 0.5 50 200
## 7 7 7.1 3 30 0.5 50 200
## 8 8 8.1 9 30 0.5 50 100
## 9 9 9.1 3 5 0.05 350 100
## 10 10 10.1 9 5 0.05 350 200
## 11 11 11.1 3 30 0.05 350 200
## 12 12 12.1 9 30 0.05 350 100
## 13 13 13.1 3 5 0.5 350 200
## 14 14 14.1 9 5 0.5 350 100
## 15 15 15.1 3 30 0.5 350 100
## 16 16 16.1 9 30 0.5 350 200
## 17 17 1.2 3 5 0.05 50 200
## 18 18 2.2 9 5 0.05 50 100
## 19 19 3.2 3 30 0.05 50 100
## 20 20 4.2 9 30 0.05 50 200
## 21 21 5.2 3 5 0.5 50 100
## 22 22 6.2 9 5 0.5 50 200
## 23 23 7.2 3 30 0.5 50 200
## 24 24 8.2 9 30 0.5 50 100
## 25 25 9.2 3 5 0.05 350 100
## 26 26 10.2 9 5 0.05 350 200
## 27 27 11.2 3 30 0.05 350 200
## 28 28 12.2 9 30 0.05 350 100
## 29 29 13.2 3 5 0.5 350 200
## 30 30 14.2 9 5 0.5 350 100
## 31 31 15.2 3 30 0.5 350 100
## 32 32 16.2 9 30 0.5 350 200
## class=design, type= FrF2
## NOTE: columns run.no and run.no.std.rp are annotation,
## not part of the data frameInput de resposta
Criando os vetores de resposta.
A nível didático, realizaremos a análise apenas para a remoção de óleo do tipo mineral (SMO).
*Dica: Para facilitar a manipulação dos dados e criar os vetores ou planilhas de maneira correspondente ao experimento, basta visualizar os tratamentos e organizar os dados de modo correspondente!
Oleo1 <- c(94.80, 23.40, 99.60, 38.00, 80.20, 27.00, 96.40, 60.00, 99.90, 9.10,
98.40, 37.10, 99.77, 51.70, 97.40, 58.00, 95.10, 25.00, 99.50, 39.10,
80.60, 28.00, 96.70, 61.20, 99.70, 10.00, 98.80, 37.40, 99.40, 52.00,
97.70, 58.40)Unindo o design às respostas
plan.atualizado1 <- add.response(design = plan.person, response = Oleo1)
# Visualizando o DoE e respostas
# View(plan.atualizado1)
summary(plan.atualizado1)## Call:
## FrF2(nfactors = 5, resolution = 5, replications = 2, randomize = FALSE,
## factor.names = list(pH = c(3, 9), Temp = c(5, 30), Dose = c(0.05,
## 0.5), Conc = c(50, 350), Velo = c(100, 200)))
##
## Experimental design of type FrF2
## 16 runs
## each run independently conducted 2 times
##
## Factor settings (scale ends):
## pH Temp Dose Conc Velo
## 1 3 5 0.05 50 100
## 2 9 30 0.50 350 200
##
## Responses:
## [1] Oleo1
##
## Design generating information:
## $legend
## [1] A=pH B=Temp C=Dose D=Conc E=Velo
##
## $generators
## [1] E=ABCD
##
##
## Alias structure:
## [[1]]
## [1] no aliasing among main effects and 2fis
##
##
## The design itself:
## run.no run.no.std.rp pH Temp Dose Conc Velo Oleo1
## 1 1 1.1 3 5 0.05 50 200 94.80
## 2 2 2.1 9 5 0.05 50 100 23.40
## 3 3 3.1 3 30 0.05 50 100 99.60
## 4 4 4.1 9 30 0.05 50 200 38.00
## 5 5 5.1 3 5 0.5 50 100 80.20
## 6 6 6.1 9 5 0.5 50 200 27.00
## 7 7 7.1 3 30 0.5 50 200 96.40
## 8 8 8.1 9 30 0.5 50 100 60.00
## 9 9 9.1 3 5 0.05 350 100 99.90
## 10 10 10.1 9 5 0.05 350 200 9.10
## 11 11 11.1 3 30 0.05 350 200 98.40
## 12 12 12.1 9 30 0.05 350 100 37.10
## 13 13 13.1 3 5 0.5 350 200 99.77
## 14 14 14.1 9 5 0.5 350 100 51.70
## 15 15 15.1 3 30 0.5 350 100 97.40
## 16 16 16.1 9 30 0.5 350 200 58.00
## 17 17 1.2 3 5 0.05 50 200 95.10
## 18 18 2.2 9 5 0.05 50 100 25.00
## 19 19 3.2 3 30 0.05 50 100 99.50
## 20 20 4.2 9 30 0.05 50 200 39.10
## 21 21 5.2 3 5 0.5 50 100 80.60
## 22 22 6.2 9 5 0.5 50 200 28.00
## 23 23 7.2 3 30 0.5 50 200 96.70
## 24 24 8.2 9 30 0.5 50 100 61.20
## 25 25 9.2 3 5 0.05 350 100 99.70
## 26 26 10.2 9 5 0.05 350 200 10.00
## 27 27 11.2 3 30 0.05 350 200 98.80
## 28 28 12.2 9 30 0.05 350 100 37.40
## 29 29 13.2 3 5 0.5 350 200 99.40
## 30 30 14.2 9 5 0.5 350 100 52.00
## 31 31 15.2 3 30 0.5 350 100 97.70
## 32 32 16.2 9 30 0.5 350 200 58.40
## class=design, type= FrF2
## NOTE: columns run.no and run.no.std.rp are annotation,
## not part of the data frameGerando gráficos de efeitos principais
Utilizando a função padrão MEPlot para gerar gráficos dos efeitos principais para cada fator.
MEPlot(plan.atualizado1, abbrev = 5, cex.xax = 1.6, cex.main = 2, main = "Efeitos Principais")
Salvando a imagem automaticamente no diretório.
png('test.png', units = "in", width = 13, height = 7, res = 300)
MEPlot(plan.atualizado1, abbrev = 5, cex.xax = 1.6, cex.main = 2)
dev.off()png 2
# É possível salvar a imagem em outros tipos de extensão, como .tiff. Basta
# substituir o png por tiff na função.
Observe que o pH é o fator de maior impacto no experimento. Quando o pH está em nível baixo é possível obter altos valores de remoção de óleo.
Matriz de efeitos de interação
# IAPlot(plan.atualizado1)
IAPlot(plan.atualizado1, abbrev = 5, show.alias = TRUE, lwd = 2,
cex = 2, cex.xax = 1.2, cex.lab = 1.5, main = "Matriz dos Efeitos de Interação")
Observe que os efeitos entre os fatores não são tão impactantes, Conc:Dose e pH:Dose são as interações que aparentam ter maior efeito sinérgico.
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(Oleo1 ~ plan.atualizado1$Conc,
groups = plan.atualizado1$Dose,
type = "a",
ylab = "Remoção de óleo (%)",
xlab = "Concentração",
auto.key = list(space = "right",
points = FALSE,
title = "Dose [g]",
cex = 0.7,
lines = TRUE))
plot2 <- xyplot(Oleo1 ~ plan.atualizado1$pH,
groups = plan.atualizado1$Dose,
type = "a",
ylab = "Remoção de óleo (%)",
xlab = "pH",
auto.key = list(space = "right",
points = FALSE,
title = "Dose [g]",
cex = 0.7,
lines = TRUE))
print(plot1, split = c(1, 1, 1, 2), more = TRUE)
print(plot2, split = c(1, 2, 1, 2))
# As funções residentes plot(effect(x)) e interaction.plot também geram gráficos similares.
# with(plan.atualizado1,interaction.plot(plan.atualizado1$pH,
# plan.atualizado1$Dose,
# plan.atualizado1$Oleo1,
# lwd=4,
# col=c("orange","purple"),
# lty=c(2,2)))
A partir do gráfico da interação Concentração:Dose, observamos que o processo de remoção de óleo é muito pouco sensível à concentração do óleo (Conc) quando a dose está em seu nível mais baixo. E um pouco mais sensível a concentração quando a dose está em seu nível mais alto. Com a dose em seu nível baixo, o processo produz uma remoção de óleo médio de 60%, e com a dose alta e concentração alta, é possível remover aproximadamente 80% do óleo. Sendo assim, é interessante avaliar fatores como custo do adsorvente para verificar seu uso no nível alto ou não.
Do segundo gráfico, é possível extrair a seguinte tabela:
| FATOR | DOSE | |
|---|---|---|
| FATOR pH | -1 | +1 |
| -1 | 99 | 90 |
| +1 | 30 | 50 |
Onde:
\[
pH\cong \frac{30+50}{2}-\frac{99+90}{2}=\frac{80}{2}-\frac{189}{2}=-\frac{109}{2}=-54.5
\]
e a resposta da variação da dose é:
\[ Dose\cong \frac{90+50}{2}-\frac{99+30}{2}=5.5 \]
E a resposta da interação é:
\[ pH\times Dose\cong \frac{90+30}{2}-\frac{99+50}{2}=-14.5 \]
Note que no nível baixo de Dose o pH tem efeito de \(E_{pH}=30-99=-60\) e no nível alto de Dose é de \(E_{pH}=50-90=-40\), sendo assim, o efeito de pH depende claramente do Dose e a diferença entre os dois valores de \(E_{pH}\) ratifica a existência da interação.
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.
Daniel plot - Gráfico half normal de 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(plan.atualizado1,
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.
De acordo com o gráfico half-normal, apenas o Fator A (pH) apresenta significância.
Pareto
O gráfico de Pareto é um gráfico de efeitos e interações padronizados, permitindo que os fatores e interações mais importantes no processo em análise sejam identificados. Ele exibe os valores absolutos dos efeitos e desenha uma linha de referência (t crítico).
Para este desenvolvimento, é necessário o uso do pacote qualityTools, pois este modelo gráfico ainda não está disponível ao pacote FrF2. É interessante verificar o design dos planejamentos para ratificação da ordem das respostas no vetor resposta.
A função utilizada para utilizar o design é a fracDesign, observe que esta função é similar a função FrF2do pacote FrF2. Entretanto, agora devemos inserir nossa gerador.
fat5_1 <- fracDesign(k = 5,
gen = "E=ABCD",
centerCube = 0,
replicates = 2)
runOrd(fat5_1) <- standOrd(fat5_1)
# fat5_1
Aglutinando as informações de resposta do nosso vetor Oleo1 ao design através da função response. E gerando o gráfico através da função paretoPlot.
response(fat5_1) <- Oleo1
# fat5_1
paretoPlot(fat5_1,
ylim = c(0,350),
ylab = "Standadized Effects",
xlab = "Term",
main = "Pareto Chart of the Standadized Effects")
É possível verificar através deste modelo gráfico que todos os fatores principais e interações são significantes, pois excedem a linha de referência. Para maiores informações sobre este modelo gráfico acesse o descritivo do pacote qualityTools e ou o info do Minitab.
ANOVA
Para gerar a análise de variância (ANOVA) será necessário o uso da função residente aov. Entretanto, é conveniente facilitar a escrita da equação realizando algumas simples manipulações e modificações.
Para destrinchar o planejamento experimental (plan.atualizado1) utilizaremos o modo clássico utilizando o $ e uma maneira mais dinâmica através da função residente attach.
attach(plan.atualizado1)
# SEM UTILIZAR attach()
A <- pH
B <- Temp
# UTILIZANDO attach(plan.atualizado1)
C <- Dose
D <- Conc
E <- Velo
# Veja que criamos os vetores A, B, C, D e E, referentes as informações contidas em plan.atualizado1
Observe que cada coluna do data frame plan.atualizado1 representa também um vetor (a informação do data frame não foi perdida), agora podemos trabalhar com vetores!
anova1 <- aov(Oleo1 ~ A*B*C*D*E)
summary(anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 26368 26368 103696.84 < 2e-16 ***
## B 1 1225 1225 4819.51 < 2e-16 ***
## C 1 609 609 2394.01 < 2e-16 ***
## D 1 113 113 444.94 4.20e-13 ***
## E 1 96 96 377.60 1.49e-12 ***
## A:B 1 512 512 2012.60 < 2e-16 ***
## A:C 1 1442 1442 5671.93 < 2e-16 ***
## B:C 1 8 8 32.37 3.35e-05 ***
## A:D 1 41 41 160.78 9.22e-10 ***
## B:D 1 175 175 687.06 1.43e-14 ***
## C:D 1 367 367 1443.31 < 2e-16 ***
## A:E 1 344 344 1354.16 < 2e-16 ***
## B:E 1 58 58 229.67 6.55e-11 ***
## C:E 1 14 14 55.08 1.45e-06 ***
## D:E 1 22 22 87.15 7.09e-08 ***
## Residuals 16 4 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(2,2))
plot(anova1)
Observe que a partir da nossa ANOVA podemos afirmar com nível de significância de 5% que todos os fatores e as interações são significantes, diferentemente do que apresentado pelo gráfico half-normal.
Analisando os resíduos do modelo
A linearidade dos residuos pode ser verificada visualmente através do gráfico Residuals vs Fitted. Como os resíduos igualmente espalhados ao redor de uma linha horizontal sem padrões distintos, isso é uma boa indicação de que você não tem relacionamentos não-lineares.
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 (gráfico das raízes quadradas dos valores absolutos dos resíduos padronizados vs valores ajustados), mostra se os resíduos são distribuídos igualmente ao longo dos intervalos de preditores. Sendo assim, é possível verificar a suposição de variância igual (homocedasticidade), não sugerindo uma eventual transformação dos dados.
O gráfico Residuals vs Leverage (influência) nos ajuda a encontrar casos influentes (isto é, sujeitos), se houver. Nem todos os outliers são influentes na análise de regressão linear (o que quer dizer que outliers significam).
Vale a pena ressaltar que as análises gráficas devem ser acompanhadas de análises numéricas de testes, como o teste de normalidade de Shapiro-Wilk dos resíduos. O Shapiro test é baseado no seguinte teste de hipótese:
\[ \left\{\begin{matrix} H_{0}: A\ amostra\ provém\ de\ uma\ população\ normal\\ H_{1}: A\ amostra\ não\ provém\ de\ uma\ população\ normal \end{matrix}\right. \]
shapiro.test(anova1$residuals)##
## Shapiro-Wilk normality test
##
## data: anova1$residuals
## W = 0.9708, p-value = 0.5217
Como o p-value é maior que o alfa, podemos afirmar com nível de significância de 5% que a amostra provém de uma população normal.
Regressão
Gerando o modelo de regressão linear através da função residente lm para o teste de significância dos coeficientes.
rg <- lm(plan.atualizado1)
summary(lm(plan.atualizado1))## Number of observations used: 32
## Formula:
## Oleo1 ~ (pH + Temp + Dose + Conc + Velo)^2
##
## Call:
## lm.default(formula = fo, data = model.frame(fo, data = formula))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8000 -0.1888 0.0000 0.1888 0.8000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.16781 0.08914 753.497 < 2e-16 ***
## pH1 -28.70531 0.08914 -322.020 < 2e-16 ***
## Temp1 6.18844 0.08914 69.423 < 2e-16 ***
## Dose1 4.36156 0.08914 48.929 < 2e-16 ***
## Conc1 1.88031 0.08914 21.094 4.20e-13 ***
## Velo1 -1.73219 0.08914 -19.432 1.49e-12 ***
## pH1:Temp1 3.99906 0.08914 44.862 < 2e-16 ***
## pH1:Dose1 6.71344 0.08914 75.312 < 2e-16 ***
## pH1:Conc1 -1.13031 0.08914 -12.680 9.22e-10 ***
## pH1:Velo1 -3.28031 0.08914 -36.799 < 2e-16 ***
## Temp1:Dose1 0.50719 0.08914 5.690 3.35e-05 ***
## Temp1:Conc1 -2.33656 0.08914 -26.212 1.43e-14 ***
## Temp1:Velo1 1.35094 0.08914 15.155 6.55e-11 ***
## Dose1:Conc1 3.38656 0.08914 37.991 < 2e-16 ***
## Dose1:Velo1 0.66156 0.08914 7.421 1.45e-06 ***
## Conc1:Velo1 -0.83219 0.08914 -9.336 7.09e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5043 on 16 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9997
## F-statistic: 8231 on 15 and 16 DF, p-value: < 2.2e-16# par(mfrow=c(2,2))
# plot(rg)
# Os gráficos obtidos com o objeto gerado a partir da função lm
# são equivalentes aos gerados pela função rg.
# O comando par(mfrow=c(x, y)) é similar ao print com uso do split.
É possível verificar que todos os efeitos e interações são significantes, assim como visto na ANOVA com a estatística F. O efeito principal tem interpretação direta, por exemplo, o efeito do Fator pH é o dobro do coeficiente de regressão estimado que é mostrado, ou seja:
\[ EA=2\times \hat{\beta} _{A} =2\times \left (-28.70531 \right )=-57.41062 \]
Isso significa que, em média, quando o pH é elevado de 3 para 9 o percentual de remoção de óleo diminuirá em -57.41062 %, como verificado no plot de interações.
O modelo de regressão proposto é o do próprio artigo, sendo a equação:
\[ Eq.=113-9.57\times pH+0.495\times Temperature+19.4\times Dose+0.0125\times Concentration-0.0346\times Speed \]
Observe o quão importante é ter dominar o objeto de estudo. Por parcimônia e dado a complexidade de ajuste baseado na análise do teste T de lm(plan.atualizado1) não foi possível ajustar o modelo num primeiro instante. Pode ser comum ajustar vários modelos aos dados, resultando em múltiplas ANOVAs.
Logo, uma nova análise voltada apenas para os modelos de primeira ordem é essencial para verificação de ajuste do modelo.
Entretanto é necessário realizar um pequeno ajuste a nível de programação e manipulação de dados no nosso código. Transformaremos, primeiro, nossas informações referentes aos fatores em valores numéricos (do modo atual está como fatores).
d <- within(plan.atualizado1, {
pH <- as.numeric(as.character(pH))
Temp <- as.numeric(as.character(Temp))
Dose <- as.numeric(as.character(Dose))
Conc <- as.numeric(as.character(Conc))
Velo <- as.numeric(as.character(Velo)) })
Desenvolvendo nosso modelo de ajuste apenas para os efeitos principais (primeira ordem).
# model linear ~ (em função de) ..., considerando 'd'
nmodel <- (lm(Oleo1 ~ (pH + Temp + Dose + Conc + Velo), d))
# summary(nmodel)
Para estimar valores dado nosso modelo de ajuste, criaremos vetores sequenciais para os fatores em análise do experimento.
# pH_rg = seq(de, até, com passo de)
pH_rg = seq(3,9,1)
# Temperatura_rg = seq(de, até, com comprimento de vetor igual a 'pH_rg')
Temperatura_rg = seq(5,30,length.out = length(pH_rg))
Dose_rg = seq(0.05,0.5,length.out = length(pH_rg))
Conc_rg = seq(50,350,length.out = length(pH_rg))
Velo_rg = seq(100,200,length.out = length(pH_rg))
Agora vamos estruturar nosso data frame (dados) com os valores obtidos nos vetores sequenciais através da função expand.grid.
dados <- expand.grid(pH = pH_rg,
Temp = Temperatura_rg,
Dose = Dose_rg,
Conc = Conc_rg,
Velo = Velo_rg)
# View(dados)
Dado os valores dos fatores em análise e conhecendo nosso modelo de regressão, vamos gerar os valores através da função predict. Note que através de dados$p já criaremos a coluna p (referente a predição) no nosso data frame dados.
dados$p <- predict(nmodel, newdata=dados)
# View(dados)
A nível de ilustração e para melhor análise da interação dos fatores, realizaremos plots em 3D e superficies de resposta.
Para plotar nosso superfície de contorno usaremos a função contourplotdo pacote laticce.
contourplot(p ~ pH * Dose, data = dados,
cuts = 20,
region = TRUE,
xlab = "pH",
ylab = "Dose [g]",
main = "Percentual de Remoção de Óleo")
# 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(p ~ pH * Dose, data = dados,
cuts = 30,
region = TRUE,
xlab = "pH",
ylab = "Dose [g]",
main = "Percentual de Remoção de Óleo",
col.regions = heat.colors(100, alpha = 1))
contourplot(p ~ pH * Dose, data = dados,
cuts = 30,
region = TRUE,
xlab = "pH",
ylab = "Dose [g]",
main = "Percentual de Remoção de Óleo",
col.regions = terrain.colors(100, alpha = 1))
Com a função wireframe é possível gerar nosso modelo de superfície de resposta.
wireframe(p ~ pH * Dose, data = dados,
panel.aspect = 0.9,
zoom = 0.7,
screen = list(z = 230, x = -75),
scales = list(arrows = FALSE),
drape = TRUE,
colorkey = TRUE)
wireframe(p ~ pH * Dose, data = dados,
panel.aspect = 0.5,
zoom = 0.8,
screen = list(z = 230, x = -75),
scales = list(arrows = FALSE),
drape = TRUE,
col.regions = heat.colors(100, alpha = 1))
wireframe(p ~ pH * Dose, data = dados,
panel.aspect = 0.5,
zoom = 0.8,
screen = list(z = 230, x = -75),
scales = list(arrows = F, cex = .8, col = "black", font = 3),
drape = TRUE,
xlab = expression(paste(pH)),
ylab = expression(paste('Dose [g]')),
zlab = list(label = "% Remoção Óleo", font = 12, cex = .9),
col.regions = topo.colors(100, alpha = 1))
wireframe(p ~ pH * Dose, data = dados,
panel.aspect = 0.5,
zoom = 0.8,
screen = list(z = 230, x = -75),
scales= list(arrows = F, cex= .9, col = "black", font = 3),
drape = TRUE,
xlab = expression(paste(pH)),
ylab = expression(paste('Dose [g]')),
zlab = list(label = "% Remoção Óleo", 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.
Com esta superfície de resposta podemos visualizar que realmente os maiores percentuais de remoção, tendo em vista os fatores pH e Dose, são alcançados nos níveis baixos de pH.
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
Ajuste do modelo de equação
Outra maneira de realizar a predição é utilizando uma function. Entretanto este modo não é tão “didático” quando o uso do predict, pois é necessário inserir informações no script aumentando o trabalho necessário. Observe a equação abaixo, é um modelo geral, mas vamos considerar apenas dois fatores como significativos.
\[
Eq. regressão = \beta _{0}+\beta _{1}\cdot x_{1}+\beta _{2}\cdot x_{2}+...
\]
Por exemplo, considere que na análise da estatística t o coeficiente \(\beta _{0}\) é 281.6, -39.2 o valor de \(\beta _{1}\) e -2.901 o valor de \(\beta _{2}\). A função para realizar a predição ficará da seguinte maneira:
# f = function(FatorA, FatorB) 281.6-39.2*FatorA-2.901*FatorB
O entendimento de f é análogo a gerar um objeto, entretanto ele sempre será em função das variáveis em análise, que no exemplo são os Fatores A e B.
No ajuste de dados para substituir a função predict, cria-se esta função f e posteriormente atribuímos os valores calculados ao grid de valores a serem testados, de maneira análoga ao uso da predict.
# dados$p <- f(FatorA, FatorB)
Lembrando que o modelo com predict é:
# dados$p <- predict(nmodel, newdata=dados)
Outro ponto negativo do uso de uma function é a possível perda de informação dado o possível arredondamento dos valores exibidos na estatística t.
Obter apenas 01 imagem da análise de resíduos
Na análise de resíduos é possível obter um gráfico por vez utilizando a função plot e determinando qual o número (which) do gráfico a ser exibido.
# plot(rg, which = 1)
plot(rg, which = 3)
Lembrando que os gráficos principais são os de número 1, 2 e 3. que se solicitado apenas plot(rg) aparecerá no console a seguinte frase: Hit
E neste link tem um modelo de como organizar a matriz de sinais produzida através do pacote FrF2 e do pacote qualityTools.