Vamos utilizar um exemplo que se encontra no livro “Design and Analysis of Experiments” (MONTGOMERY, 2001, pg 482), em que uma engenheira química está interessada em determinar em que condições de operações a produção de um processo é maximizado. Duas variáveis controláveis influencia o processo de produção: reação do tempo (x1) e reação de temperatura (x2). A engenheira está atualmente operando o processo com uma reação do tempo de 35 minutos e uma temperatura de 155ºF, com um resultado na produção de 40%. É improvável que esta contenha a região ótima, ela ajusta um modelo de primeira ordem.

#install.packages("rsm")
require(rsm)
## Loading required package: rsm
## Warning: package 'rsm' was built under R version 3.5.3
attach(ChemReact)

Primeiramente, ajustou-se um modelo de primeira de primeira ordem, com intuito de verificar se a verdadeira função que relaciona a variável resposta com os fatores, Temperatura e Tempo aproxima-se de um função linear.

Ajustando um modelo de superfície de resposta

Neste experimento, os dados no bloco B1 foram coletados primeiro e analisados, após os quais o bloco B2 foi adicionado e uma nova análise foi feita.

Na maioria dos casos, os cálculos para estimação dos parâmetros podem ser simplificados codificando os níveis das k variáveis independentes \(X_{i}\) Para criar um conjunto de dados codificado com o códigos, utilizaremos a fórmula:

CR <- coded.data(ChemReact, x1 ~(Time - 85)/5, x2 ~(Temp - 175)/5)
CR[1:7,]
##   Time Temp Block Yield
## 1   80  170    B1  80.5
## 2   80  180    B1  81.5
## 3   90  170    B1  82.0
## 4   90  180    B1  83.5
## 5   85  175    B1  83.9
## 6   85  175    B1  84.3
## 7   85  175    B1  84.0
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Time - 85)/5
## x2 ~ (Temp - 175)/5

coded.data (transforma os valores preditos e substitui essas variáveis para a versão codificada)

Na primeira etapa, ajusta-se um modelo de primeira ordem, utilizando os seguintes comandos:

CR.rsm1 <- rsm(Yield~FO(x1,x2),data = CR,subset = (Block =="B1"))
summary(CR.rsm1)
## 
## Call:
## rsm(formula = Yield ~ FO(x1, x2), data = CR, subset = (Block == 
##     "B1"))
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 82.81429    0.54719 151.3456 1.143e-08 ***
## x1           0.87500    0.72386   1.2088    0.2933    
## x2           0.62500    0.72386   0.8634    0.4366    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.3555, Adjusted R-squared:  0.0333 
## F-statistic: 1.103 on 2 and 4 DF,  p-value: 0.4153
## 
## Analysis of Variance Table
## 
## Response: Yield
##             Df Sum Sq Mean Sq F value  Pr(>F)
## FO(x1, x2)   2 4.6250  2.3125  1.1033 0.41534
## Residuals    4 8.3836  2.0959                
## Lack of fit  2 8.2969  4.1485 95.7335 0.01034
## Pure error   2 0.0867  0.0433                
## 
## Direction of steepest ascent (at radius 1):
##        x1        x2 
## 0.8137335 0.5812382 
## 
## Corresponding increment in original units:
##     Time     Temp 
## 4.068667 2.906191

Devido à falta de ajuste significativa (p-valor = 0,01), deve-se então utilizar um modelo de ordem maior. O modelo ajustado não apresenta caracteristicas que nos permita estimar valores para a variável resposta. Como não há termos quadráticos e de interação no ajuste do modelo, nota-se que as variáveis x1 e x2 não são significativas.

Incluem-se, então, interações:

CR.rsm1.5 <- update(CR.rsm1,.~.+TWI(x1,x2))
summary(CR.rsm1.5)
## 
## Call:
## rsm(formula = Yield ~ FO(x1, x2) + TWI(x1, x2), data = CR, subset = (Block == 
##     "B1"))
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 82.81429    0.62948 131.5604 9.683e-07 ***
## x1           0.87500    0.83272   1.0508    0.3705    
## x2           0.62500    0.83272   0.7506    0.5074    
## x1:x2        0.12500    0.83272   0.1501    0.8902    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.3603, Adjusted R-squared:  -0.2793 
## F-statistic: 0.5633 on 3 and 3 DF,  p-value: 0.6755
## 
## Analysis of Variance Table
## 
## Response: Yield
##             Df Sum Sq Mean Sq  F value   Pr(>F)
## FO(x1, x2)   2 4.6250  2.3125   0.8337 0.515302
## TWI(x1, x2)  1 0.0625  0.0625   0.0225 0.890202
## Residuals    3 8.3211  2.7737                  
## Lack of fit  1 8.2344  8.2344 190.0247 0.005221
## Pure error   2 0.0867  0.0433                  
## 
## Stationary point of response surface:
## x1 x2 
## -5 -7 
## 
## Stationary point in original units:
## Time Temp 
##   60  140 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1]  0.0625 -0.0625
## 
## $vectors
##         [,1]       [,2]
## x1 0.7071068 -0.7071068
## x2 0.7071068  0.7071068

Verifica-se novamente para este caso que a falta de ajuste é significativa, com p-valor = 0.005. Os dados do bloco 2 são acrescentados a fim de montar um modelo de segunda ordem. Isso pode ser feito utilizando “SO(x1,x2)”“, que inclui termos quadráticos e interação:

CR.rsm2 <- rsm(Yield ~ Block + SO(x1, x2), data = CR)
summary(CR.rsm2)
## 
## Call:
## rsm(formula = Yield ~ Block + SO(x1, x2), data = CR)
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 84.095427   0.079631 1056.067 < 2.2e-16 ***
## BlockB2     -4.457530   0.087226  -51.103 2.877e-10 ***
## x1           0.932541   0.057699   16.162 8.444e-07 ***
## x2           0.577712   0.057699   10.012 2.122e-05 ***
## x1:x2        0.125000   0.081592    1.532    0.1694    
## x1^2        -1.308555   0.060064  -21.786 1.083e-07 ***
## x2^2        -0.933442   0.060064  -15.541 1.104e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9981, Adjusted R-squared:  0.9964 
## F-statistic: 607.2 on 6 and 7 DF,  p-value: 3.811e-09
## 
## Analysis of Variance Table
## 
## Response: Yield
##             Df Sum Sq Mean Sq   F value    Pr(>F)
## Block        1 69.531  69.531 2611.0950 2.879e-10
## FO(x1, x2)   2  9.626   4.813  180.7341 9.450e-07
## TWI(x1, x2)  1  0.063   0.063    2.3470    0.1694
## PQ(x1, x2)   2 17.791   8.896  334.0539 1.135e-07
## Residuals    7  0.186   0.027                    
## Lack of fit  3  0.053   0.018    0.5307    0.6851
## Pure error   4  0.133   0.033                    
## 
## Stationary point of response surface:
##        x1        x2 
## 0.3722954 0.3343802 
## 
## Stationary point in original units:
##      Time      Temp 
##  86.86148 176.67190 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.9233027 -1.3186949
## 
## $vectors
##          [,1]       [,2]
## x1 -0.1601375 -0.9870947
## x2 -0.9870947  0.1601375

Obtém-se agora falta de ajuste não significativa (p-valor = 0,69), ou seja, o modelo de segunda ordem se ajustou bem aos dados. Observa-se também que o ponto de estacionariedade do modelo ajustado está em (0,37; 0,33), que é ponto de máximo.

Para construir as curvas de níveis e a superfície de resposta, pode ser empregada a função lm (linear model) ou a função rsm (response surface methodology).

Gráfico de contorno

fit.model <- rsm(Yield~SO(Time,Temp))
contour(fit.model, ~Time+Temp, main="Gráfico de contornos", col = rainbow(40))

Percebe-se que quando os valores de x1(TIME) se aproxima de 85 min e os valores de x2(TEMP) se aproxima de 175ºF os valores da variável resposta é maximizado. Como foi dito anteriormente podemos localizar de forma mais precisa o ponto de estacionariedade obtido através da ANOVA, que foi x1 = 86.86148 e x2 = 176.67190.

e por fim o gráfico de Superfície de Resposta

persp(fit.model, ~Time+Temp, main="Superfície de resposta",
       col=rainbow(40))